text:
#=##############################################################################
# DESCRIPTION
    Methods for statistical analysis on geometric datasets
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2021
  * License   : MIT License
=###############################################################################

# ----------- CARTESIAN OPERATIONS ---------------------------------------------
_agglomerate_mean(name, din, dout, args...) = din[name]
function _agglomerate_variance(name, din, dout, args...; mean_suf="-mean")
    if length(size(din[name]))==1
        return ( (din[name][i]-dout[name*mean_suf][i])^2
                        for i in 1:length(din[name]))
    else
        return ( (din[name][i,j]-dout[name*mean_suf][i,j])^2
                       for i in 1:size(din[name],1), j in 1:size(din[name], 2) )
    end
end

"Correlation of k-th direction with all other directions"
function _agglomerate_correlation(k::Int, name, din, dout, args...; mean_suf="-mean")
    return (
                (din[name][k, j]-dout[name*mean_suf][k, j])*(din[name][(i-1)%3+1, j]-dout[name*mean_suf][(i-1)%3+1, j])
                        for i in 1:size(din[name], 1), j in 1:size(din[name], 2)
           )
end
_agglomerate_correlation1(args...; optargs...) = _agglomerate_correlation(1, args...; optargs...)
_agglomerate_correlation2(args...; optargs...) = _agglomerate_correlation(2, args...; optargs...)
_agglomerate_correlation3(args...; optargs...) = _agglomerate_correlation(3, args...; optargs...)

# ----------- CYLINDRICAL OPERATIONS -------------------------------------------
"Calculate radial and azimuthal direction at each centroid"
function _get_rad_azim_dirs(axial_dir, din)

    keyi = findfirst(x->contains(x, "centroid"), collect(keys(din)))
    if keyi == nothing
        error("Centroids not found."*
        " Make sure that centroids are calculated before calling this function")
    end

    # Fetch position of each element
    centroid_key = collect(keys(din))[keyi]
    centroids = din[centroid_key]
    nelems = size(centroids, 2)

    # Determine radial directions
    radial_dirs_aux = ( ((centroids[i, j] - dot(view(centroids, :, j), axial_dir)*axial_dir[i])
                                            for i in 1:3) for j in 1:nelems )
    radial_dirs = ( collect(R)/norm(R) for R in radial_dirs_aux )

    # Determine azimuthal directions
    azim_dirs = ( cross(axial_dir, R) for R in radial_dirs )

    return radial_dirs, azim_dirs
end

function _agglomerate_cyl_mean(axial_dir::Array{<:Real, 1}, name, din, args...)

    radial_dirs, azim_dirs = _get_rad_azim_dirs(axial_dir, din)

    if length(size(din[name]))==1
        # error("Cylindrical transformation requested on a scalar field!")
        return _agglomerate_mean(name, din, args...)
    end

    data = din[name]
    ndims = size(data, 1)
    nelems = size(data, 2)
    axial_dirs = (axial_dir for i in 1:nelems)

    return (
                dot( view(data, :, j[1]), j[i+1] )
                    for i in 1:ndims,
                        j in zip(1:nelems, radial_dirs, azim_dirs, axial_dirs)
            )
end

function _agglomerate_cyl_variance(axial_dir::Array{<:Real, 1}, name, din, dout, args...; mean_suf="-cyl-mean")

    radial_dirs, azim_dirs = _get_rad_azim_dirs(axial_dir, din)

    if length(size(din[name]))==1
        # error("Cylindrical transformation requested on a scalar field!")
        return _agglomerate_variance(name, din, dout, args...; mean_suf=mean_suf)
    end

    Din = din[name]
    Dout = dout[name*mean_suf]
    ndims = size(Din, 1)
    nelems = size(Din, 2)
    axial_dirs = (axial_dir for i in 1:nelems)

    return (
                (dot( view(Din, :, j[1]), j[i+1] ) - Dout[i,j[1]])^2
                    for i in 1:ndims,
                        j in zip(1:nelems, radial_dirs, azim_dirs, axial_dirs)
            )
end

"Correlation of k-th direction with all other directions"
function _agglomerate_cyl_correlation(k::Int, axial_dir::Array{<:Real, 1}, name, din, dout, args...; mean_suf="-cyl-mean")

    radial_dirs, azim_dirs = _get_rad_azim_dirs(axial_dir, din)

    if length(size(din[name]))==1
        # error("Cylindrical transformation requested on a scalar field!")
        return _agglomerate_correlation(k, name, din, dout, args...; mean_suf=mean_suf)
    end

    Din = din[name]
    Dout = dout[name*mean_suf]
    ndims = size(Din, 1)
    nelems = size(Din, 2)
    axial_dirs = (axial_dir for i in 1:nelems)

    return (
                (dot( view(Din, :, j[1]), j[k+1] ) - Dout[k,j[1]]) * (dot( view(Din, :, j[1]), j[(i-1)%3+1+1] ) - Dout[(i-1)%3+1, j[1]])
                    for i in 1:ndims,
                        j in zip(1:nelems, radial_dirs, azim_dirs, axial_dirs)
            )
end

_agglomerate_cyl_correlation1(args...; optargs...) = _agglomerate_cyl_correlation(1, args...; optargs...)
_agglomerate_cyl_correlation2(args...; optargs...) = _agglomerate_cyl_correlation(2, args...; optargs...)
_agglomerate_cyl_correlation3(args...; optargs...) = _agglomerate_cyl_correlation(3, args...; optargs...)

# ----------- STATISTICAL OPERATIONS ON VTKS -----------------------------------
function calculate_statistics_vtk(vtkfiles;                  # List of files to process
                                    read_path="",            # Folder where to read files from
                                    ites=2,                  # Number of processing iterations
                                                             # Operations to perform; (operation name, operation function, iteration to act on)
                                    operations=[("mean", _agglomerate_mean, 1), ("variance", _agglomerate_variance, 2)],
                                    add_centroid=true,       # Whether to calculate and operate on centroid
                                    verbose=true, v_lvl=0)

    data_out = Dict()
    ncells = nothing
    centroid_cells = nothing

    nfiles = length(vtkfiles)

    for ite in 1:ites

        if verbose
           println("\t"^(v_lvl)*"Iteration $(ite) out of $(ites)")
        end

        for (fi, filename) in enumerate(vtkfiles)

            if verbose && fi%(ceil(Int, nfiles/10))==0
               println("\t"^(v_lvl+1)*"Reading file $(fi) out of $(nfiles)")
            end

            # Read VTK
            points, cells, cell_types, data = read_vtk(filename; path=read_path)

            # Pre-allocate and initialize output data with zeroes
            if ite==1 && fi==1
                ncells = length(cells)
                centroid_cells = zeros(3, ncells)
            end

            # Add centroid to original dataset
            if add_centroid

                # Point centroid
                if !("POINT_DATA" in keys(data))
                    data["POINT_DATA"] = Dict()
                end
                data["POINT_DATA"]["centroid-point"] = points

                # Cell centroid
                if length(cells) != 0

                    centroid_cells .= 0
                    for (celli, cell) in enumerate(cells)
                        for pi in cell
                            centroid_cells[:, celli] .+= view(points, :, pi+1)
                        end
                        centroid_cells[:, celli] ./= ncells
                    end

                    if !("CELL_DATA" in keys(data))
                        data["CELL_DATA"] = Dict()
                    end
                    data["CELL_DATA"]["centroid-cell"] = centroid_cells
                end
            end

            # Pre-allocate and initialize output data with zeroes
            if ite==1 && fi==1

                for (dataparent, subdata) in data
                    data_out[dataparent] = Dict()

                    for (name, this_data) in subdata
                        for (op_name, op) in operations
                            data_out[dataparent][name*"-"*op_name] = zero(this_data)
                        end
                    end
                end

            end

            # Iterate over data sets
            for (dataparent, subdata) in data
                for (name, this_data) in subdata

                    # Apply operation to agglomerate each dataset
                    for (op_name, op, op_ite) in operations
                        if op_ite==ite
                            try
                                data_out[dataparent][name*"-"*op_name] .+= op(name, subdata, data_out[dataparent], points, cells)
                            catch e
                                println(op)
                                println(size(data_out[dataparent][name*"-"*op_name]))
                                println(size(collect(op(name, subdata, data_out[dataparent], points, cells))))
                                println()
                                throw(e)
                            end
                        end
                    end

                end
            end

            # Convert the agglomerates into averages
            if fi==nfiles
                for (dataparent, subdata) in data
                    for (name, this_data) in subdata
                        for (op_name, op, op_ite) in operations
                            if op_ite==ite
                                data_out[dataparent][name*"-"*op_name] ./= nfiles
                            end
                        end
                    end
                end
            end

        end

    end

    return data_out
end

function calculate_statistics_vtk(vtkfiles, save_path;
                                    read_path="",
                                    out_filename="statistics",
                                    verbose=true, v_lvl=0,
                                    optargs...)

    data_out = calculate_statistics_vtk(vtkfiles;
                                        read_path=read_path,
                                        verbose=verbose, v_lvl=v_lvl,
                                        optargs...)
    str = nothing

    if save_path != nothing

        if verbose
           println("\t"^(v_lvl)*"Saving statistics under '$(out_filename).vtk'")
        end

        # Read last VTK as default geometry
        points, cells, cell_types, data = read_vtk(vtkfiles[end]; path=read_path)

        # Re-format points
        points = [points[:, pi] for pi in 1:size(points, 2)]

        point_data, cell_data = nothing, nothing

        if "POINT_DATA" in keys(data_out)
            point_data = [
                            Dict(
                                "field_name" => name,
                                "field_type" => length(size(this_data))==1 ? "scalar" : "vector",
                                "field_data" => length(size(this_data))==1 ? this_data : [this_data[:, pi] for pi in 1:size(this_data, 2)]
                            ) for (name, this_data) in data_out["POINT_DATA"]
                        ]
        end

        if "CELL_DATA" in keys(data_out)
            cell_data = [
                            Dict(
                                "field_name" => name,
                                "field_type" => length(size(this_data))==1 ? "scalar" : "vector",
                                "field_data" => length(size(this_data))==1 ? this_data : [this_data[:, celli] for celli in 1:size(this_data, 2)]
                            ) for (name, this_data) in data_out["CELL_DATA"]
                        ]
        end

        # Generate VTK file
        str = generateVTK(out_filename, points;
                                cells=cells,
                                point_data=point_data,
                                cell_data=cell_data,
                                path=save_path)
    end

    return data_out, str
end


text:
using DataFrames, XLSX

# Definimos la estructura de TrainInfo
struct TrainInfo
    train_id::String
    train_info::DataFrame
    team_id::String
    train_team::DataFrame
    train_player::DataFrame
    train_physical::Vector{DataFrame}
    train_bpm::Vector{DataFrame}
end

# Funcion de lectura del fichero TrainInfo
function leer_TrainInfo(path::String)
    # Leemos fichero
    xf = XLSX.readxlsx(path)
    sps_names = XLSX.sheetnames(xf)

    # Leemos primera hoja
    info_match = xf[sps_names[1]]
    # Convertimos a DataFrame
    df_train = DataFrame(info_match[:], :auto)

    # Obtenemos la tabla de información del entrenamiento
    column_names_train = Vector{String}(df_train[1, :])
    df_train_info = df_train[1:2, :]
    rename!(df_train_info, column_names_train)
    delete!(df_train_info, 1)
    train_id = df_train_info[1, 1]

    # Obtenemos la tabla de información del team
    column_names_team = Vector{String}(df_train[5, 1:2])
    df_train_team = df_train[5:6, 1:2]
    rename!(df_train_team, column_names_team)
    delete!(df_train_team, 1)
    team_id = df_train_team[1, 1]

    # Obtenemos la tabla de información de los jugadores
    column_names_player = Vector{String}(df_train[9, 1:7])
    df_train_player = df_train[10:end, 1:7]
    rename!(df_train_player, column_names_player)

    ## Obtenemos las tablas de datos del entrenamiento
    # Datos físicos
    dfs_physical = Vector{DataFrame}()
    # Datos de los latidos por minuto
    dfs_bpm = Vector{DataFrame}()
    # Recorremos las hojas que contienen estas tablas
    for i in 2:length(sps_names)
        # Leemos la hoja i
        info_players = xf[sps_names[i]]
        df_players = DataFrame(info_players[:], :auto)

        # Al haber dos tablas divididas en la misma hoja de tamaño
        # variable, calculamos la división dependiendo de los
        # valores missing.
        player_id = df_players[2:end, 1]
        replace!(player_id, missing => 0)
        d = findfirst(x -> x == 0, player_id)

        # Obtenemos la tabla de los datos físicos
        column_names_physical = Vector{String}(df_players[1, :])
        df_train_physical = df_players[2:d, :]
        rename!(df_train_physical, column_names_physical)
        push!(dfs_physical, df_train_physical)

        # Obtenemos la tabla de los datos de latidos por minuto
        column_names_bpm = Vector{String}(player_id[1:d-1])
        df_train_bpm = df_players[d+5:end, 1:d-1]
        rename!(df_train_bpm, column_names_bpm)
        push!(dfs_bpm, df_train_bpm)
    end
    # Devolvemos la información pertienente
    TrainInfo(train_id, df_train_info, team_id, df_train_team, df_train_player, dfs_physical, dfs_bpm)
end

ext: jl

max_stars_repo_name: vntr-CC/Forward-Football

max_issues_repo_name: vntr-CC/Forward-Football

max_forks_repo_name: vntr-CC/Forward-Football


text:
#

"""
    `as_poly(f)` convert something into a polynomial

Here `f` can be a `Polynomial` object; a vector of coefficents in the
form `[a_0, a_1, ..., a_n]`; or a callable object, such as a function,
that implements a polynomial function. The latter is determined by whether
it can be evaluated on the `Polynomial` monomial `x`.
"""
as_poly(f::Poly) = f
as_poly(T, f::Poly) = convert(Poly{T}, f)

as_poly(xs::Vector{T})  where {T}= Poly(xs)
as_poly(S::T, xs::Vector) where {T} = Poly(convert(Vector{S},xs))

## Try to convert a callable object into a polynomial, `Poly{T}`. `T` can be specified, or guessed from calling `f(0)`.
function as_poly(f)
    T = typeof(f(0))
    as_poly(T, f)
end

"""
as_poly{T}(T, f)

Convert `f` to a polynomial of type `Poly{T}`, or throw a `DomainError`.
"""
function as_poly(T, f)
    p = try
        x = variable(T)
        convert(Poly{T}, f(x))
    catch err
        throw(DomainError)
    end
    p
end

"""

Find coefficients of polynomial expressed as Poly, Callable object, or values [a0,a1, ..., an]

"""
poly_coeffs(ps::Vector{T}) where {T} = ps
poly_coeffs(p::Poly{T}) where {T} = Polynomials.coeffs(p)
poly_coeffs(f) = poly_coeffs(as_poly(f))
poly_coeffs(T, f) = convert(Vector{T}, poly_coeffs(f))

" Type of polynomial "
e_type(p::Poly{T}) where {T} = T
e_type(ps::Vector{T}) where {T} = T
e_type(p) =  eltype(p(0))

"""
reverse coefficients of a polynomial
"""
rcoeffs(p::Poly) = reverse(coeffs(p))

"""
monic
"""
monic(p::Poly) = p[Polynomials.degree(p)] != 0 ? Poly(p.a * inv(p[Polynomials.degree(p)]), p.var) : p

# a robust bisection method
bisection(f, a::Real, b::Real) = Roots.bisection64(f, promote(float(a), b)...)


text:
abstract type AbstractRecursiveVector{T,L} <: StaticArrays.StaticArray{Tuple{L}, T, 1} end
abstract type MutableRecursiveVector{T,L} <: AbstractRecursiveVector{T,L} end

struct StaticRecursiveVector{T,L} <: AbstractRecursiveVector{T,L}
    data::NTuple{L,T}
    StaticRecursiveVector{T,L}(data::NTuple{L,T}) where {T,L} = new(data)
end
mutable struct RecursiveVector{T,L} <: MutableRecursiveVector{T,L}
    data::NTuple{L,T}
    RecursiveVector{T,L}(data::NTuple{L,T}) where {T,L} = new(data)
    RecursiveVector{T,L}() where {T,L} = new()
    function RecursiveVector{T,L}(x::AbstractVector{T}) where {T,L}
        out = RecursiveVector{T,L}()
        @inbounds for l ∈ 1:L
            out[l] = x[l]
        end
        out
    end
end
struct PointerRecursiveVector{T,M,N,L} <: MutableRecursiveVector{T,L}
    data::Ptr{StaticRecursiveVector{T,cutoff}}
    # PointerRecursiveVector(ptr::Ptr{StaticRecursiveVector{T,cutoff,cutoff,cutoff2}}, ::Val{M}, ::Val{N}, ::Val{L}) where {T,M,N,L} = new{T,M,N,L}(ptr)
end

randvec(n) = randvec(Val(n))
function randvec(::Val{L}) where L
    out = RecursiveVector{Float64,L}()
    @inbounds for i ∈ 1:L
        out[i] = randn()
    end
    out
end

@inline float_point(A::RecursiveVector{T}) where T = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
@inline float_point(A::PointerRecursiveVector{T}) where T = Base.unsafe_convert(Ptr{T}, A.data)

Base.similar(::RecursiveVector{T,L}) where {T,L} = RecursiveVector{T,L}()
Base.similar(::RecursiveVector{T,L}, ::Type{Tnew}) where {T,L,Tnew} = RecursiveVector{Tnew,L}()
@inline Base.size(::RecursiveVector{T,L}) where {T,L} = (L,)
@inline Base.length(::RecursiveVector{T,L}) where {T,L} = L
@inline function Base.getindex(x::RecursiveVector{T,L}, i::Int) where {T,L}
    @boundscheck i > L && throw(BoundsError())
    x.data[i]
end
function Base.getindex(x::RecursiveVector{T,L}, i, ::Val{R}) where {T,L,R}
    @boundscheck i > L && throw(BoundsError())
    # ntuple(j -> x.data[i+j-1].value, Val(R))
    RT = Ptr{StaticSIMD{R,T}}
    Base.unsafe_load(Base.unsafe_convert(RT, pointer_from_objref(x)), i)
end
@inline function Base.setindex!(x::RecursiveVector{T,L}, val, i::Int) where {T,L}
    @boundscheck i > L && throw(BoundsError())
    unsafe_store!(float_point(x), convert(T, val), i )
    val
end

function Base.copyto!(x::RecursiveVector{T,L}, y::RecursiveVector{T,L}) where {T,L}
    @inbounds for l ∈ 1:L
        x[l] = y[l]
    end
    x
end
Base.copy(x::RecursiveVector{T,L}) where {T,L} = copyto!(RecursiveVector{T,L}(), x)

function LinearAlgebra.scale!(x::RecursiveVector{T,L}, α::Number) where {T,L}
    @inbounds @simd for l ∈ 1:L
        x[l] *= α
    end
    x
end

Base.vec(x::RecursiveVector) = x

struct StaticSIMD{N,T} <: AbstractVector{T}
    data::NTuple{N,Core.VecElement{T}}
end
Base.size(::StaticSIMD{N}) where N = (N,)
Base.length(::StaticSIMD{N}) where N = N
function Base.getindex(x::StaticSIMD{N}, i::Int) where N
    @boundscheck i > N && throw(BoundsError())
    x.data[i].value
end
@generated function Base.:*(x::StaticSIMD{N,T}, y::StaticSIMD{N,T}) where {N,T}
    quote
        Base.@_inline_meta
        StaticSIMD{$N,$T}( $(Expr(:tuple, [:(x.data[$i].value * y.data[$i].value) for i ∈ 1:N]...)) )
    end
end
@generated function Base.:+(x::StaticSIMD{N,T}, y::StaticSIMD{N,T}) where {N,T}
    quote
        Base.@_inline_meta
        StaticSIMD{$N,$T}( $(Expr(:tuple, [:(x.data[$i].value + y.data[$i].value) for i ∈ 1:N]...)) )
    end
end
@generated function Base.fma(x::StaticSIMD{N,T}, y::StaticSIMD{N,T}, z::StaticSIMD{N,T}) where {N,T}
    quote
        Base.@_inline_meta
        StaticSIMD{$N,$T}( $(Expr(:tuple, [:(x.data[$i].value * y.data[$i].value + z.data[$i].value) for i ∈ 1:N]...)) )
    end
end
@generated function Base.sum(x::StaticSIMD{N,T}) where {N,T}
    quote
        Base.@_inline_meta
        $(Expr(:call, :+, [:(x.data[$i].value) for i ∈ 1:N]...))
    end
end
@generated function LinearAlgebra.dot(x::RecursiveVector{T,L}, y::RecursiveVector{T,L}) where {T,L}
    vec_length = 8
    VVal = Val(vec_length)
    Lo4, r = divrem(L, vec_length)
    if L >= vec_length
        q = quote
            # @inbounds out_tup = x[1:4] .* y[1:4]
            @inbounds out_tup = x[1, $VVal] * y[1, $VVal]
            @fastmath for i ∈ 2:$Lo4
                @inbounds out_tup = fma(x[i, $VVal], y[i, $VVal], out_tup)
            end
            out = sum(out_tup)
            @inbounds for i ∈ $(L-r+1):$L
                out += x[i] * y[i]
            end
            out
        end
    elseif L > 0
        q, qa = create_quote()
        push!(qa, :(out = x[1] * y[1]))
        for i ∈ 2:L
            push!(qa, :(out += x[$i] * y[$i]))
        end
        push!(q.args, :(out))
    else
        q = :(NaN)
    end
    q
end

@generated function abs2norm(x::RecursiveVector{T,L}) where {T,L}
    vec_length = 8
    VVal = Val(vec_length)
    Lo4, r = divrem(L, vec_length)
    quote
        # @inbounds out_tup = x[1:4] .* y[1:4]
        simdv = x[1, $VVal]
        @inbounds out_tup = simdv * simdv
        @fastmath for i ∈ 2:$Lo4
            simdv = x[i, $VVal]
            @inbounds out_tup = fma(simdv, simdv, out_tup)
        end
        out = sum(out_tup)
        @inbounds for i ∈ $(L-r+1):$L
            out += abs2(x[i])
        end
        out
    end
end


text:
"""
    set_tokenizer(fun)

Call this to set the default tokenizer to invoke the passed in function `fun`
It will be used by `tokenize`.
Calling this will trigger recompilation of any functions that use `tokenize`.

Calling `set_tokenizer`  will give method overwritten warnings. They are expected, be worried if they do not occur
"""
function set_tokenizer(fun)
    @eval tokenize(str::AbstractString) = $(fun)(str)
end

"""
    set_sentence_splitter(fun)

Call this to set the default sentence splitter to invoke the passed in function `fun`
It will be used by `split_sentences`.
Calling this will trigger recompilation of any functions that use `split_sentences`.

Calling `set_sentence_splitter`  will give method overwritten warnings. They are expected, be worried if they do not occur
"""
function set_sentence_splitter(fun)
    @eval split_sentences(str::AbstractString) = $(fun)(str)
end

set_tokenizer(toktok_tokenize)
set_sentence_splitter(rulebased_split_sentences)


text:
"""
module Restarttools

Primary tools to do calculations.
Also contains all methods to restart calculations.
"""
module Restarttools

export calculate_adaptive_sample_inputs,
       calculate_energy_for_xyzfile,
       calculate_with_different_method,
       calculate_potential,
       continue_calculation,
       create_inputs,
       load_clusters_and_make_input,
       load_clusters_and_sample_input,
       load_data_file,
       load_restart_file,
       test_work,
       write_restart_file,
       write_save_file

using ..Calculators
using ..Clusters
using ..Fileaccess
using ..Sample
using Distributed
using ProgressMeter

"""
    write_save_file(fname, calculator, points, energy, cluster1, cluster2)

Saves final information for energy calculation.

# Arguments
- `fname` : name of restartfileusing PotentialCalculation
- `calculator::Calculator` : calculator used in calculations
- `points` : 2d array of point where to calculate energy
- `energy` : energy for points that have been caculeted
- `cluster1` : cluster1
- `cluster2` : cluster2
"""
function write_save_file(
    fname::AbstractString,
    calculator::Calculator,
    points,
    energy,
    cluster1,
    cluster2
)
    data = Dict("Energy" => energy, "Points" => points, "cluster1" => cluster1, "cluster2" => cluster2,
                "Basis" => calculator.basis, "Method" => calculator.method )
    flush(stdout)
    save_jld_data(fname, data)
end

"""
    write_restart_file(fname, calculator, points, restart_energy, cluster1, cluster2)

Saves restart information for energy calculation.

# Arguments
- `fname` : name of restartfile
- `calculator::Calculator` : calculator used in calculations
- `points` : 2d array of point where to calculate energy
- `energy` : energy for points that have been caculeted
- `cluster1` : cluster1
- `cluster2` : cluster2
"""
function write_restart_file(
    fname::AbstractString,
    calculator::Calculator,
    points,
    restart_energy,
    cluster1,
    cluster2
)
    @info "Writing restart file $(fname)"
    data = Dict("restart_energy" => restart_energy, "Points" => points, "cluster1" => cluster1, "cluster2" => cluster2,
                "Basis" => calculator.basis, "Method" => calculator.method )
    flush(stdout)
    save_jld_data(fname, data)
end

"""
    load_data_file(fname)

Loads saved data
"""
function load_data_file(fname::AbstractString)
    return load_jld_data(fname)
end

"""
    load_restart_file(fname)

Loads restart information from file `fname` and adds to it key `not_calculated`,
which holds information of which collumns of `Points` have not been calculated.
"""
function load_restart_file(fname::AbstractString)
    data = load_jld_data(fname)
    if haskey(data, "restart_energy")
        if ! ( typeof(data["restart_energy"]) <: Tuple)
            #NOTE Old style restart file
            not_calculated = Set(1:size(data["Points"])[2])
            setdiff!(not_calculated, x[1] for x in data["restart_energy"])
            push!(data, "not_calculated"=>not_calculated)
            iscalculated = trues(size(data["Points"]))
            for i in not_calculated
                iscalculated[:,i] .= false
            end
            energy = similar(data["Points"],Float64)
            for (i,e) in data["restart_energy"]
                energy[:,i] .= e
            end
            data["restart_energy"] = (energy, iscalculated)
        end
    else
        @warn "File is not a restart file" fname
    end
    return data
end

"""
continue_calculation(fname, calculator::Calculator; save_file="", restart_file="",
                              save_after=nworkers(),  pbar=true)

Restarts calculation from given file

# Arguments
- `fname` : file from which calculation is restarted
- `calculator::Calculator` : calculator used in calculations
- `save_file` : file where final results are saved, if given
- `restart_file` : file where new restart information is saved, if given
- `save_after=nworkers()` : make restart file when given ammount of points is calculated
- `pbar=true` : show progress bar
"""
function continue_calculation(
    fname::AbstractString,
    calculator::Calculator;
    save_file="",
    restart_file="",
    save_after=nworkers(),
    pbar=true
)
    data = load_restart_file(fname)
    @info "File $(fname) loaded - continuing calculation"
    calculator.basis = data["Basis"]
    calculator.method = data["Method"]

    energy,iscal = data["restart_energy"]

    lc1 = length(data["cluster1"])
    lc2 = length(data["cluster2"])
    rc1 =  1:lc1           #UnitRange for cluster1
    rc2 =  lc1+1:lc1+lc2   #UnitRange for cluster2

    lcal = length(data["Points"])-length(data["Points"][iscal])

    if pbar
        prog = Progress(lcal  ,dt=1, desc="Calculating points:")
    end

    jobs = RemoteChannel(()->Channel(2*nworkers()))
    results = RemoteChannel(()->Channel(2*nworkers()))

    iscalculated = deepcopy(iscal)
    @async for i in eachindex(data["Points"])
        if ! iscal[i]
            x = (calculator,data["Points"][i][rc1], data["Points"][i][rc2])
            put!(jobs, (i,x) )
        end
    end

    for p in workers()
        remote_do(parallel_work, p, jobs, results, _calculate_points)
    end

    n = 0
    for j in 1:lcal
        i,tmp = take!(results)
        energy[i] = tmp
        iscalculated[i] = true
        pbar && ProgressMeter.next!(prog)
        n += 1
        if restart_file != "" && n >= save_after && j < lcal
            write_restart_file(restart_file, calculator, data["Points"], (energy,iscalculated),
                               data["cluster1"], data["cluster2"])
            n = 0
        end
    end

    if save_file != ""
        write_save_file(save_file, calculator, data["Points"], energy,
                           data["cluster1"], data["cluster2"])
    end
    return Dict("Energy" => energy, "Points"=> data["Points"],
               "cluster1"=>data["cluster1"], "cluster2"=>data["cluster2"],
               "Method"=>calculator.method, "Basis"=>calculator.basis)

end

function continue_calculation(
    fname::AbstractString,
    calculator::AbstractCalculationProgram;
    save_file="",
    restart_file="",
    save_after=nworkers(),
    pbar=true
)
    return continue_calculation(
        fname,
        Calculator("","",calculator);
        save_file=save_file,
        restart_file=restart_file,
        save_after=save_after,
        pbar=pbar
    )
end

"""
    calculate_potential(fname::AbstractString, calculator::Calculator;
                        save_file="", restart_file="", pbar=true, save_after=nworkers()
                       )
With this function you can use already chosen points on to which to do energy calculation.

# Arguments
- `fname` : name of save/restart file where points are read
- `calculatror::Calculator` : calculator used for calculations

# Keywords
- `save_file=""` : save final results here, if not empty
- `restart_file=""` : save restarts here, if not empty
- `pbar=true` : show progress bar
- `save_after=nworkers()` : make restart file when given ammount of points is calculated
"""
function calculate_potential(fname::AbstractString, calculator::Calculator;
                             save_file="", restart_file="", pbar=true, save_after=nworkers()
                            )
    data = load_jld_data(fname)
    @info "File $(fname) loaded - calculating with different method"
    flush(stdout)
    lc1 = length(data["cluster1"])
    c1_points = map( x-> x[1:lc1], data["Points"])
    c2_points = map( x-> x[lc1+1:end], data["Points"])

    if pbar
        prog = Progress(length(c1_points) ,dt=1, desc="Calculating points:")
    end

    jobs = RemoteChannel(()->Channel(2*nworkers()))
    results = RemoteChannel(()->Channel(2*nworkers()))

    @async for i in eachindex(c1_points)
        x = (calculator,c1_points[i], c2_points[i])
        put!(jobs, (i,x) )
    end

    for p in workers()
        remote_do(parallel_work, p, jobs, results, _calculate_points)
    end

    energy = similar(c1_points, Float64)
    iscalculated = falses(size(energy))
    n = 0
    lcal=length(c1_points)
    for j in 1:lcal
        i,tmp = take!(results)
        energy[i] = tmp
        iscalculated[i] = true
        pbar && ProgressMeter.next!(prog)
        n += 1
        if restart_file != "" && n >= save_after && j < lcal
            write_restart_file(restart_file, calculator, data["Points"], (energy,iscalculated),
                               data["cluster1"], data["cluster2"])
            n = 0
        end
    end

    if save_file != ""
        write_save_file(save_file, calculator, data["Points"], energy,
                           data["cluster1"], data["cluster2"])
    end
    return Dict("Energy" => energy, "Points"=> data["Points"],
               "cluster1"=>data["cluster1"], "cluster2"=>data["cluster2"],
               "Method"=>calculator.method, "Basis"=>calculator.basis)
end

function calculate_with_different_method(fname::AbstractString, calculator::Calculator;
                             save_file="", restart_file="", pbar=true, save_after=nworkers()
                            )
    @warn "calculate_with_different_method is deprecated use calculate_potential instead"
    return calculate_potential(fname, calculator; save_file=save_file, restart_file=restart_file,
                               pbar=pbar, save_after=save_after)
end

"""
load_clusters_and_make_input(cluster1::String, cluster2::Cluster, calculator, nsamples;
                                      max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)

Loads cluster1 from xyz-file and returns them all as an Array
that can then be used with [`calculate_adaptive_sample_inputs`](@ref) to calculate energy data.
Differs from [`load_clusters_and_sample_input`](@ref) by taking every point from file

# Arguments
- `cluster1::String` : file from where cluster1 is sampled
- `cluster2::Cluster` :  cluster2
- `max_e` : Point that is closest and has less energy than this will be starting point for a line
- `unit` : Unit in which `max_e` is given
- `npoints` : Number of points in potential
- `maxdis` : Maximum distance in potential calculation
- `sstep` : Search step size for adaptive search
- `startdistance` : Distance from which adaptive seach is started
"""
function load_clusters_and_make_input(cluster1::String, cluster2::Cluster, calculator;
                                      max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)
    @warn "load_clusters_and_make_input deprecated use create_inputs instead"
    cluster1 = read_xyz(cluster1)

    return  [InputAdaptiveSampler(calculator, c, cluster2,
                1, max_e, unit=unit, npoints=npoints, maxdis=maxdis,
                sstep=sstep, startdistance=startdistance) for c in cluster1 ]
end

"""
load_clusters_and_make_input(cluster1::Cluster, cluster2::Cluster, calculator;
                                      nlines=1, max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)

Loads cluster1 from xyz-file and returns them all as an Array
that can then be used with [`calculate_adaptive_sample_inputs`](@ref) to calculate energy data.
Differs from [`load_clusters_and_sample_input`](@ref) by taking every point from file

# Arguments
- `cluster1::Cluster` : cluster1
- `cluster2::Cluster` : cluster2
- `nlines` : number of lines to be sampled in calculation
- `max_e` : Point that is closest and has less energy than this will be starting point for a line
- `unit` : Unit in which `max_e` is given
- `npoints` : Number of points in potential
- `maxdis` : Maximum distance in potential calculation
- `sstep` : Search step size for adaptive search
- `startdistance` : Distance from which adaptive seach is started
"""
function load_clusters_and_make_input(cluster1::Cluster, cluster2::Cluster, calculator;
                                      nlines=1, max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)

    @warn "load_clusters_and_make_input deprecated use create_inputs instead"
    return  [InputAdaptiveSampler(calculator, cluster1, cluster2,
                1, max_e, unit=unit, npoints=npoints, maxdis=maxdis,
                sstep=sstep, startdistance=startdistance) for _ in 1:nlines ]
end

"""
load_clusters_and_sample_input(fname_cluster1, cluster2, calculator, nsamples;
                                      max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)

Loads cluster1 from xyz-file and takes `nsamples` samples of it and returns them as an Array
that can then be used with [`calculate_adaptive_sample_inputs`](@ref) to calculate energy data

# Arguments
- `cluster1` : file from where cluster1 is sampled
- `cluster2` : cluster2
- `nsamples` : number of samples taken from `fname_cluster1`
- `max_e` : Point that is closest and has less energy than this will be starting point for a line
- `unit` : Unit in which `max_e` is given
- `npoints` : Number of points in potential
- `maxdis` : Maximum distance in potential calculation
- `sstep` : Search step size for adaptive search
- `startdistance` : Distance from which adaptive seach is started
"""
function load_clusters_and_sample_input(cluster1::String, cluster2, calculator, nsamples;
                                      max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)
    @warn "load_clusters_and_sample_input deprecated use create_inputs instead"
    cluster1 = read_xyz(cluster1)

    return  [InputAdaptiveSampler(calculator, rand(cluster1), cluster2,
                1, max_e, unit=unit, npoints=npoints, maxdis=maxdis,
                sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
end

"""
function load_clusters_and_sample_input(cluster1::String, cluster2::String, calculator, nsamples;
                                      nlines=1, max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)

# Arguments
- `cluster1` : file from where cluster1 is sampled
- `cluster2` : file from where cluster2 is sampled
- `nsamples` : number of samples taken from `fname_cluster1`
- `nlines` : number of lines to be sampled in calculation
- `max_e` : Point that is closest and has less energy than this will be starting point for a line
- `unit` : Unit in which `max_e` is given
- `npoints` : Number of points in potential
- `maxdis` : Maximum distance in potential calculation
- `sstep` : Search step size for adaptive search
- `startdistance` : Distance from which adaptive seach is started
"""
function load_clusters_and_sample_input(cluster1::String, cluster2::String, calculator, nsamples;
                                      nlines=1, max_e=0, unit="cm-1", npoints=10,
                                      maxdis=9.0, sstep=0.1, startdistance=2.5)
    @warn "load_clusters_and_sample_input is deprecated use create_inputs instead"
    cluster1 = read_xyz(cluster1)
    cluster2 = read_xyz(cluster2)

    return  [InputAdaptiveSampler(calculator, rand(cluster1), rand(cluster2),
                nlines, max_e, unit=unit, npoints=npoints, maxdis=maxdis,
                sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
end

"""
    calculate_adaptive_sample_inputs(inputs; save_after="", save_after=nworkers(), pbar=true)

Uses [`adaptive_line_sampler`](@ref) to `inputs` in distributed fashion

# Arguments
- `inputs` : calculation inputs array (eg. from [`create_inputs`](@ref))
- `save_after` : file where data saved
- `save_after=nworkers()` : number of calculated items before data is saved
- `pbar=true` : show progress bar
"""
function calculate_potential(inputs::AbstractArray{InputAdaptiveSampler};
                             save_file="", save_after=nworkers(), pbar=true)
    @info "Starting adaptive sampling"
    flush(stdout)
    t= collect(1:save_after:length(inputs))
    i_range = [ t[i-1]:t[i]-1  for i in 2:length(t) ]
    if t[end] < length(inputs)
        push!(i_range, t[end]:length(inputs))
    end

    if pbar
        prog = Progress(length(inputs) ,dt=1, desc="Calculating points:")
    end

    jobs = RemoteChannel(()->Channel(2*nworkers()))
    results = RemoteChannel(()->Channel(2*nworkers()))

    @async for x in inputs
        put!(jobs, (1, (x,) ) )
    end

    for p in workers()
        remote_do(parallel_work, p, jobs, results, _sample_and_calculate)
    end

    _,tmp = take!(results)
    energy = tmp["Energy"]
    points = tmp["Points"]
    mindis = tmp["Mindis"]
    pbar && ProgressMeter.next!(prog)
    if save_file != ""
        @info "Saving data to file $(save_file)"
        write_save_file(save_file, inputs[1].cal, points, energy,
                        inputs[1].cl1, inputs[1].cl2)
    end
    if length(inputs) > 1
        for i in 2:length(inputs)
            _,tmp = take!(results)
            energy = hcat(energy, tmp["Energy"] )
            points = hcat(points, tmp["Points"] )
            mindis = vcat(mindis, tmp["Mindis"] )
            pbar && ProgressMeter.next!(prog)
            if save_file != ""
                @info "Saving data to file $(save_file)"
                write_save_file(save_file, inputs[1].cal, points, energy,
                                inputs[1].cl1, inputs[1].cl2)
            end
        end
    end
    #pbar && ProgressMeter.finish!(prog)
    return Dict("Energy" => energy, "Points" => points, "Mindis" => mindis,
                "cluster1"=>inputs[1].cl1, "cluster2"=>inputs[1].cl2,
                "Method"=>inputs[1].cal.method, "Basis"=>inputs[1].cal.basis)
end

function calculate_adaptive_sample_inputs(inputs::AbstractArray{InputAdaptiveSampler};
                                         save_file_name="", save_step=nworkers(), pbar=true)
    @warn "calculate_adaptive_sample_inputs is deprecated use calculate_potential instead"
    return calculate_potential(inputs, save_file=save_file_name, save_after=save_step, pbar=pbar)
 end

"""
    create_inputs(cluster1, cluster2, calculator::Calculator;  kwargs...)

Creates inputs for the calculations. Use [`calculate_adaptive_sample_inputs`](@ref) on
results to do the calculation itself.

Inputs for the molecules (named cluster1/2) can be loaded from xyz-file or given by as
[`Cluster`](@ref) structures. If given as file names, then `nsamples` points are randomly
picked from the trajectories for calculation. Total number of lines sampled is thus
`nlines` times `nsamples`.

Parallelisation is done over `nsamples`, thus it is recommended for it to be multiple of
number workers processes.

# Arguments
- `cluster1` : something that can be interpreted as [`Cluster`](@ref) or an Array of it
- `cluster2` : something that can be interpreted as [`Cluster`](@ref) or an Array of it
- `calculator::Calculator` : [`Calculator`](@ref) used in sampling

# Keywords
- max_e=15000 : energy treshold in sampling. Enegy is lower than this value.
- maxdis=9.0 : maximun distance in calculation
- nline=1  : number of lines sampled for given structure pair
- npoints=10 : number of points in each sampled line
- nsamples=2 : number of points picked randomly from trajectories for line sampling
- sstep=0.1 : step size used in sampling
- startdistance= 3.5 : starting distance in sampling
- unit="cm-1" : unit for `max_e`

"""
function create_inputs(cluster1::AbstractString,
                cluster2::AbstractString,
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    c1 = read_xyz(cluster1)
    c2 = read_xyz(cluster2)

    return  create_inputs(c1, c2, calculator;
                nlines=nlines,
                nsamples=nsamples,
                max_e=max_e,
                unit=unit,
                npoints=npoints,
                maxdis=maxdis,
                sstep=sstep,
                startdistance=startdistance
            )
end

function create_inputs(cluster1::AbstractString,
                cluster2,
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    c1 = read_xyz(cluster1)

    return  create_inputs(c1, cluster2, calculator;
                nlines=nlines,
                nsamples=nsamples,
                max_e=max_e,
                unit=unit,
                npoints=npoints,
                maxdis=maxdis,
                sstep=sstep,
                startdistance=startdistance
            )
end

function create_inputs(cluster1,
                cluster2::AbstractString,
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    c2 = read_xyz(cluster2)

    return  create_inputs(cluster1, c2, calculator;
                nlines=nlines,
                nsamples=nsamples,
                max_e=max_e,
                unit=unit,
                npoints=npoints,
                maxdis=maxdis,
                sstep=sstep,
                startdistance=startdistance
            )
end

function create_inputs(cluster1::AbstractCluster,
                cluster2::AbstractCluster,
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    return  [InputAdaptiveSampler(calculator, cluster1, cluster2, nlines, max_e;
                                  unit=unit, npoints=npoints, maxdis=maxdis,
                                  sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
end

function create_inputs(cluster1::AbstractCluster,
                cluster2::AbstractArray{<:AbstractCluster},
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    return  [InputAdaptiveSampler(calculator, cluster1, rand(cluster2), nlines, max_e;
                                  unit=unit, npoints=npoints, maxdis=maxdis,
                                  sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
end

function create_inputs(cluster1::AbstractArray{<:AbstractCluster},
                cluster2::AbstractCluster,
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    return  [InputAdaptiveSampler(calculator, rand(cluster1), cluster2, nlines, max_e;
                                  unit=unit, npoints=npoints, maxdis=maxdis,
                                  sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
end

function create_inputs(cluster1::AbstractArray{<:AbstractCluster},
                cluster2::AbstractArray{<:AbstractCluster},
                calculator::Calculator;
                nlines=1,
                nsamples=2,
                max_e=15000,
                unit="cm-1",
                npoints=10,
                maxdis=9.0,
                sstep=0.1,
                startdistance=3.0
               )

    return  [InputAdaptiveSampler(calculator, rand(cluster1), rand(cluster2), nlines, max_e;
                                  unit=unit, npoints=npoints, maxdis=maxdis,
                                  sstep=sstep, startdistance=startdistance) for _ in 1:nsamples ]
end

"""
    _sample_and_calculate(inputs; pchannel=undef )

Internal helper function to help implement distributed calculations
"""
function _sample_and_calculate(inputs::InputAdaptiveSampler; pchannel=undef)
    start_dir = pwd()
    if typeof(inputs.cal.calculator) == Orca
        if inputs.cal.calculator.tmp_dir != start_dir
            cd(inputs.cal.calculator.tmp_dir)
        end
    end
    return sample_multiple_adaptive_lines(inputs, basename="base-$(myid())",
                                         id="Pid $(myid())", pchannel=pchannel)
end

"""
    _calculate_points(cal, c1_points, c2_points; pchannel=undef)

Internal helper function to ease use of distributed calculations.
Users should call "calculate_points" instead.
"""
function _calculate_points(cal, c1_points, c2_points; pchannel=undef)
    return bsse_corrected_energy(cal, c1_points, c2_points, basename="base-$(myid())",
                                 id="Pid $(myid())", pchannel=pchannel)
end

function _calculate_points(cal::Calculator{Orca}, c1_points, c2_points; pchannel=undef)
    start_dir = pwd()
    if cal.calculator.tmp_dir != start_dir
        cd(cal.calculator.tmp_dir)
    end
    return bsse_corrected_energy(cal, c1_points, c2_points, basename="base-$(myid())",
                                 id="Pid $(myid())", pchannel=pchannel)
end

"""
    calculate_energy_for_xyzfile(fname, cal; pbar=true)

Reads xyz-file and calculates energy of each point on it
"""
function calculate_energy_for_xyzfile(fname, cal; pbar=true)
    println("Calculating points for a file $(fname)")
    points = read_xyz(fname)
    out=@showprogress pmap(points) do x
        _calculate_energy(cal,x)
    end
    return out
end

function _calculate_energy(cal, points; pchannel=undef)
    return calculate_energy(cal, points, basename="base-$(myid())",
                                 id="Pid $(myid())", pchannel=pchannel)
end

function _calculate_energy(cal::Calculator{Orca}, points; pchannel=undef)
    start_dir = pwd()
    if cal.calculator.tmp_dir != start_dir
        cd(cal.calculator.tmp_dir)
    end
    return calculate_energy(cal, points, basename="base-$(myid())",
                                 id="Pid $(myid())", pchannel=pchannel)
end

function parallel_work(jobs, results, f)
    while true
        i,t = take!(jobs)
        x = f(t...)
        put!(results, (i,x))
    end
end

end  # module restarttools


text:
# Autogenerated wrapper script for SPIRV_Headers_jll for any
export spirv_xml

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `spirv_xml`
const spirv_xml_splitpath = ["include", "spirv", "spir-v.xml"]

# This will be filled out by __init__() for all products, as it must be done at runtime
spirv_xml_path = ""

# spirv_xml-specific global declaration
# This will be filled out by __init__()
spirv_xml = ""

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"SPIRV_Headers")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global spirv_xml_path = normpath(joinpath(artifact_dir, spirv_xml_splitpath...))

    global spirv_xml = spirv_xml_path
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
const stageForEnding = Dict{String, vk.VkShaderStageFlagBits}([
    ("vert.spv", vk.VK_SHADER_STAGE_VERTEX_BIT),
    ("vert", vk.VK_SHADER_STAGE_VERTEX_BIT),
    ("frag.spv", vk.VK_SHADER_STAGE_FRAGMENT_BIT),
    ("frag", vk.VK_SHADER_STAGE_FRAGMENT_BIT),
    ("geom.spv", vk.VK_SHADER_STAGE_GEOMETRY_BIT),
    ("geom", vk.VK_SHADER_STAGE_GEOMETRY_BIT),
    ("tcsh.spv", vk.VK_SHADER_STAGE_TESSELLATION_CONTROL_BIT),
    ("tcsh", vk.VK_SHADER_STAGE_TESSELLATION_CONTROL_BIT),
    ("tesc.spv", vk.VK_SHADER_STAGE_TESSELLATION_CONTROL_BIT),
    ("tesc", vk.VK_SHADER_STAGE_TESSELLATION_CONTROL_BIT),
    ("tesh.spv", vk.VK_SHADER_STAGE_TESSELLATION_EVALUATION_BIT),
    ("tesh", vk.VK_SHADER_STAGE_TESSELLATION_EVALUATION_BIT),
    ("tese.spv", vk.VK_SHADER_STAGE_TESSELLATION_EVALUATION_BIT),
    ("tese", vk.VK_SHADER_STAGE_TESSELLATION_EVALUATION_BIT),
    ("comp.spv", vk.VK_SHADER_STAGE_COMPUTE_BIT),
    ("comp", vk.VK_SHADER_STAGE_COMPUTE_BIT),
    ("rgen.spv", vk.VK_SHADER_STAGE_RAYGEN_BIT_NV),
    ("rgen", vk.VK_SHADER_STAGE_RAYGEN_BIT_NV),
    ("rint.spv", vk.VK_SHADER_STAGE_INTERSECTION_BIT_NV),
    ("rint", vk.VK_SHADER_STAGE_INTERSECTION_BIT_NV),
    ("rahit.spv", vk.VK_SHADER_STAGE_ANY_HIT_BIT_NV),
    ("rahit", vk.VK_SHADER_STAGE_ANY_HIT_BIT_NV),
    ("rchit.spv", vk.VK_SHADER_STAGE_CLOSEST_HIT_BIT_NV),
    ("rchit", vk.VK_SHADER_STAGE_CLOSEST_HIT_BIT_NV),
    ("rmiss.spv", vk.VK_SHADER_STAGE_MISS_BIT_NV),
    ("rmiss", vk.VK_SHADER_STAGE_MISS_BIT_NV),
    ("rcall.spv", vk.VK_SHADER_STAGE_CALLABLE_BIT_NV),
    ("rcall", vk.VK_SHADER_STAGE_CALLABLE_BIT_NV),
    ("task.spv", vk.VK_SHADER_STAGE_TASK_BIT_NV),
    ("task", vk.VK_SHADER_STAGE_TASK_BIT_NV),
    ("mesh.spv", vk.VK_SHADER_STAGE_MESH_BIT_NV),
    ("mesh", vk.VK_SHADER_STAGE_MESH_BIT_NV)
])

function identifyShader(filename::String)::vk.VkShaderStageFlagBits
    for s in stageForEnding
        if endswith(filename, s.first)
            return s.second
        end
    end
    return vk.VK_SHADER_STAGE_ALL
end


text:
#
# WIP: work in progress on deepvalidate. This code is not expected to work at this time
#

using Distributions
using DelimitedFiles
using SemanticModels
using Suppressor
using Test

include("varextract.jl")

@testset "TraceExtract" begin
    out_file = "traces.dat";
    mode = "w"
    traces = 300;

    function add(a, b)
        c = a + b
        return c
    end
    g(x) = begin
        y = add(x.*x, -x)
        z = 1
        v = y .- z
        s = sum(v)
        return s
    end
    h(x) = begin
        z = g(x)
        zed = sqrt(z)
        return zed
    end

    # Error conditions happen when our inputs are sufficiently small, so
    # Normal(0,2) gives us a good range of values to generate a reasonable
    # percentage of "bad" traces on which to train.

    if isfile(out_file)
        mode = "a"
    end

    seeds = rand(Normal(0,2),traces,3)

    for i=1:size(seeds,1)
        orig_stdout = stdout
        (read_in, write_out) = redirect_stdout()
        ctx = Extract.TraceCtx(pass=Extract.ExtractPass, metadata = (Any[], Any[]))

        try
            result = Extract.OverDub(ctx, h, seeds[i,:])
        catch DomainError
            dump(ctx.metadata)
        finally
            println()
        end

        open(out_file, mode) do f
            write(f, read(read_in))
        end

        redirect_stdout(orig_stdout)
        close(write_out)

        if i%1000 == 0
            @info string(i)
        end
    end
end


text:
using DashBootstrapComponents

list_group = dbc_listgroup([
    dbc_listgroupitem("The primary item", color="primary"),
    dbc_listgroupitem("A secondary item", color="secondary"),
    dbc_listgroupitem("A successful item", color="success"),
    dbc_listgroupitem("A warning item", color="warning"),
    dbc_listgroupitem("A dangerous item", color="danger"),
    dbc_listgroupitem("An informative item", color="info"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# Implements passing in Jacobians (and possibly other functions) via
# function overloading:
#
# - f(...) - objective function
# - f(Val{:jac}, ...) - Jacobian of objective function
# - the details of what `...` needs to be depends on the
#   AbstractODEProblem subtype

# Method_exists does not work:
#
# julia> f(::Val{:jac}, a, b, c) = 5
# f (generic function with 1 method)
#
# julia> hasmethod(f, Tuple{Val{:jac}, Vararg})
# false
#
# Thus hand-code it:
check_first_arg(f,T::Type) = check_first_arg(typeof(f),T)
function check_first_arg(::Type{F}, T::Type) where F
    typ = Tuple{Any, T, Vararg}
    typ2 = Tuple{Any, Type{T}, Vararg} # This one is required for overloaded types
    method_table = Base.MethodList(F.name.mt) # F.name.mt gets the method-table
    for m in method_table
        (m.sig<:typ || m.sig<:typ2) && return true
    end
    return false
end
# Standard
__has_jac(f) = check_first_arg(f, Val{:jac})
__has_tgrad(f) = check_first_arg(f, Val{:tgrad})

# Performance
__has_invW(f) = check_first_arg(f, Val{:invW})
__has_invW_t(f) = check_first_arg(f, Val{:invW_t})
has_invW(f::T) where {T} = istrait(HasInvW{T})
has_invW_t(f::T) where {T} = istrait(HasInvW_t{T})

# Parameter-Based
__has_paramderiv(f) = check_first_arg(f, Val{:deriv})
__has_paramjac(f) = check_first_arg(f, Val{:paramjac})

## Parameter Names Check
__has_syms(f) = :syms ∈ fieldnames(typeof(f))

## Analytical Solution Check
__has_analytic(f) = check_first_arg(f, Val{:analytic})


text:
"""
```julia
triangulate(p1, p2, P1, P2)
```

Triangulate point given its two projection coordinates and projection matrices.

# Arguments:

- `p1`: Pixel coordinates of a point in `(x, y)` format in the first view.
- `p2`: Pixel coordinates of a point in `(x, y)` format in the second view.
- `P1`: `3x4` or `4x4` projection matrix `K * P` for the first view.
- `P2`: `3x4` or `4x4` projection matrix `K * P` for the second view.

# Returns:

Triangulated point in `(x, y, z, w)` format.
To get actual coordinates, divide by `w`.
"""
function triangulate(p1, p2, P1, P2, cache::GEEV4x4Cache)
    p, P = (p1, p2), (P1, P2)
    @inbounds A = SMatrix{4, 4, Float64, 16}(p[i][j] * P[i][3, k] - P[i][j, k] for j ∈ 1:2, i ∈ 1:2, k ∈ 1:4)
    copy!(cache.A, A' * A)
    __geev_4x4!(cache)
    @inbounds SVector{4, Float64}(@view(cache.VL[:, argmin(cache.WR)]))
end

function triangulate(p1, p2, P1, P2)
    p, P = (p1, p2), (P1, P2)
    @inbounds A = SMatrix{4, 4, Float64, 16}(p[i][j] * P[i][3, k] - P[i][j, k] for j ∈ 1:2, i ∈ 1:2, k ∈ 1:4)
    cache = GEEV4x4Cache()
    copy!(cache.A, A' * A)
    __geev_4x4!(cache)
    @inbounds SVector{4, Float64}(@view(cache.VL[:, argmin(cache.WR)]))
end

chirality_test!(inliers, points1, points2, P1, P2, K1, K2; max_repr_error = 1.0) =
    chirality_test!(inliers, points1, points2, P1, P2, K1, K2, GEEV4x4Cache(); max_repr_error)

function chirality_test!(
    inliers, pixels1, pixels2, P1, P2, K1, K2, cache::GEEV4x4Cache;
    max_repr_error = 1.0,
)
    n_points = length(pixels1)
    Pr1, Pr2 = K1 * P1, K2 * P2

    repr_error = 0.0
    n_inliers = 0
    xy_ids = SVector{2, UInt8}(1, 2)

    @inbounds for i ∈ 1:n_points
        pixel1, pixel2 = pixels1[i], pixels2[i]
        pt3d = triangulate(pixel1, pixel2, Pr1, Pr2, cache)
        if (pt3d[3] < 0 && pt3d[4] > 0) || (pt3d[3] > 0 && pt3d[4] < 0)
            inliers[i] = false
            continue
        end

        pt3d *= 1.0 / pt3d[4]
        if pt3d[3] ≥ 50
            inliers[i] = false
            continue
        end

        x2 = P2 * pt3d
        if !(0 < x2[3] < 50)
            inliers[i] = false
            continue
        end

        projection = Pr2 * pt3d
        pixel = projection[xy_ids] .* (1.0 / projection[3])
        Δ = pixel2 .- pixel
        error = √(Δ[1]^2 + Δ[2]^2)
        if error > max_repr_error
            inliers[i] = false
            continue
        end

        repr_error += error
        n_inliers += 1
    end
    repr_error /= n_inliers
    n_inliers, repr_error
end

function compute_essential_error!(inliers, p1, p2, E, threshold)
    threshold *= threshold
    n_inliers = 0
    avg_error = 0.0
    Et = E'

    @inbounds for i ∈ 1:length(p1)
        p1i = SVector{3, Float64}(p1[i]..., 1.0)
        p2i = SVector{3, Float64}(p2[i]..., 1.0)
        Ep1 = E * p1i
        Ep2 = Et * p2i

        error = p2i ⋅ Ep1
        error = (error * error) * (
            (1.0 / Ep1[1] * Ep1[1] + Ep1[2] * Ep1[2]) +
            (1.0 / Ep2[1] * Ep2[1] + Ep2[2] * Ep2[2]))
        if error < threshold
            inliers[i] = true
            n_inliers += 1
            avg_error += error
        end
    end
    avg_error /= n_inliers
    n_inliers, avg_error
end

function compute_projections(E)
    W = SMatrix{3, 3, Float64, 9}(0, 1, 0, -1, 0, 0, 0, 0, 1)

    F = svd(E; full=true)
    U, Vt = F.U, F.Vt
    det(U) < 0 && (U *= -1;)
    det(Vt) < 0 && (Vt *= -1;)

    t = @inbounds U[:, 3]
    t2 = -t
    R1 = U * W * Vt
    R2 = U * W' * Vt

    P1 = get_transformation(R1, t)
    P2 = get_transformation(R1, t2)
    P3 = get_transformation(R2, t)
    P4 = get_transformation(R2, t2)
    P1, P2, P3, P4
end


text:
"""
    WDError <: Exception
A Web Driver Exception.
"""
struct WDError <: Exception
    status::Int
    error::String
    message::String
    function WDError(err::StatusError)
        json = JSON3.read(err.response.body).value
        new(err.status, json.error, json.message)
    end
end


text:
type ConvolutionLayer
    filters::Vector{Vol}
    biases::Vol
    stride::Int
    padding::Int
    inputActivation::Vol
    outputActivation::Vol

    function ConvolutionLayer(numOutputs::Int, stride::Int, padding::Int, filterDimensions, inputDimensions)
        (filterX, filterY, filterDepth) = filterDimensions
        (inputX, inputY, inputDepth) = inputDimensions
        outputX = floor((inputX + 2 * padding - filterX) / stride + 1)
        outputY = floor((inputY + 2 * padding - filterY) / stride + 1)
        filters = [Vol(dimensions) for _ in 1:numOutputs]
        biases = Vol((1, 1, numOutputs))
        inputActivation = Vol(inputDimensions)
        outputActivation = Vol((outputX, outputY, numOutputs))
        new(filters, biases, stride, padding, inputActivation, outputActivation)
    end
end

function foward(self::ConvolutionLayer, inputActivation::Vol)
    (inputX, inputY, inputDepth) = size(l.inputActivation)
    (outputX, outputY, numOutputs) = size(l.outputActivation)
    self.inputActivation = inputActivation
    for d in 1:length(self.filters)
        filter = self.filters[d]
        filterX, filterY, filterDepth = size(filters)
        for ay in 1:outputY, ax in 1:outputX
            x, y = (stride - 1) * ax - padding, (stride - 1) * ay - padding
            a = 0.0
            for fd in 1:filterDepth, fy in 1:filterY, fx in 1:filterX
                ox = x + fx
                oy = y + fy
                if 1 <= ox <= inputX && 1 <= oy <= inputY && 1 <= fd <= inputDepth
                    a += filter.w[fx, fy, fd] * inputActivation.w[ox, oy, fd]
                end
            end
            a += self.biases.w[1, 1, d]
            self.outputActivation.w[ax, ay, d] = a
        end
    end
    return self.outputActivation
end

function backward(self::ConvolutionLayer)
    (inputX, inputY, inputDepth) = size(self.inputActivation)
    (outputX, outputY, numOutputs) = size(self.outputActivation)
    fill!(self.inputActivation.dw, 0.0)
    for d in 1:length(self.filters)
        filter = self.filters[d]
        filterX, filterY, filterDepth = size(filters)
        for ay in 1:outputY, ax in 1:outputX
            x, y = (stride - 1) * ax - padding, (stride - 1) * ay - padding
            chainGrad = this.outputActivation.dw[ax, ay, d]
            for fd in 1:filterDepth, fy in 1:filterY, fx in 1:filterX
                ox = x + fx
                oy = y + fy
                if 1 <= ox <= inputX && 1 <= oy <= inputY && 1 <= fd <= inputDepth
                    filter.dw[fx, fy, d] += self.inputActivation[ox, oy, d] * chainGrad
                    self.inputActivation.dw[ox, oy, d] += filter.w[fx, fy, d] * chainGrad
                end
            end
            self.biases.dw[1, 1, d] += chain_grad
        end
    end
end

function gradientStates(self::ConvolutionLayer)
    grads = [GradientState(filter.w, filter.dw,  1.0, 1.0) for filter in self.filters]
    append!(grads, [GradientState(self.biases.w, self.biases.dw, 0.0, 0.0)])
    return grads
end


text:
"""
    Buffer(xs, ...)

`Buffer` is an array-like type which is mutable when taking gradients. You can
construct a `Buffer` with the same syntax as `similar` (e.g. `Buffer(xs, 5)`)
and then use normal indexing. Finally, use `copy` to get back a normal array.

For example:

```julia
julia> function vstack(xs)
           buf = Buffer(xs, length(xs), 5)
           for i = 1:5
             buf[:, i] = xs
           end
           return copy(buf)
         end
vstack (generic function with 1 method)

julia> vstack([1, 2, 3])
3×5 Array{Int64,2}:
 1  1  1  1  1
 2  2  2  2  2
 3  3  3  3  3

julia> gradient(x -> sum(vstack(x)), [1, 2, 3])
([5.0, 5.0, 5.0],)
```

`Buffer` is not an `AbstractArray` and can't be used for linear algebra
operations like matrix multiplication. This prevents it from being captured by
pullbacks.

`copy` is a semantic copy, but does not allocate memory. Instead the `Buffer`
is made immutable after copying.
"""
mutable struct Buffer{T,A<:AbstractArray{T}}
  data::A
  freeze::Bool
end

Buffer(xs::AbstractArray, args...) = Buffer(similar(xs, args...), false)

bufferfrom(xs::AbstractArray) = Buffer(xs, false)

Base.getindex(b::Buffer, i...) = b.data[i...]

function Base.setindex!(b::Buffer, v, i...)
  b.freeze && error("Buffer is frozen")
  b.data[i...] = v
end

function Base.copyto!(b::Buffer, data)
  b.freeze && error("Buffer is frozen")
  copyto!(b.data, data)
end

function Base.push!(b::Buffer, data)
  b.freeze && error("Buffer is frozen")
  push!(b.data, data)
end

function Base.copy(b::Buffer)
  b.freeze = true
  return b.data
end

function Base.deleteat!(b::Buffer, i)
  b.freeze && error("Buffer is frozen")
  deleteat!(b.data, i)
  return b
end

@forward Buffer.data  Base.eltype, Base.length, Base.ndims, Base.size, Base.axes,
                      Base.eachindex, Base.stride, Base.strides, Base.findfirst,
                      Base.keys

Base.IteratorSize(::Type{<:Buffer{<:Any, A}}) where {A} = Base.IteratorSize(A)

# Buffer iteration mirrors iteration for AbstractArray
function Base.iterate(b::Buffer, state=(eachindex(b),))
  y = iterate(state...)
  y === nothing && return nothing
  b[y[1]], (state[1], Base.tail(y)...)
end

grad_mut(b::Buffer) = fill!(similar(b.data, Any), nothing)
grad_mut(b::Buffer{T}) where T<:Number = fill!(similar(b.data, float(T)), 0)

function ZygoteRules._pullback(cx::AContext, ::Type{T}, args...) where {T<:Buffer}
  T(args...), _ -> nothing
end

@adjoint function getindex(b::Buffer, i...)
  b[i...], function (Δ)
    grad = grad_mut(__context__, b)
    grad[i...] = accum(grad[i...], Δ)
    return
  end
end

@adjoint! function setindex!(b::Buffer, v, i...)
  setindex!(b, v, i...), function (_)
    grad = grad_mut(__context__, b)
    v̄ = grad[i...]
    zero = eltype(grad) <: Number ? 0 : nothing
    if i isa NTuple{N,Integer} where N
      grad[i...] = zero
    else
      grad[i...] .= zero
    end
    (nothing, v̄, map(_->nothing, i)...)
  end
end

@adjoint! function copyto!(b::Buffer, xs)
  copyto!(b, xs), function (_)
    grad = grad_mut(__context__, b)
    x̄s = copy(grad)
    grad .= eltype(grad) <: Number ? 0 : nothing
    return (nothing, x̄s)
  end
end

@adjoint! function push!(b::Buffer, x)
  push!(b, x), function (y)
    grad = grad_mut(__context__, b)
    return (nothing, pop!(grad))
  end
end

function ZygoteRules._pullback(cx::AContext,
                               ::typeof(Broadcast.materialize!),
                               b::Buffer,
                               x::AbstractArray)
  ZygoteRules._pullback(cx, copyto!, b, x)
end

@adjoint function copy(b::Buffer)
  copy(b), function (b̄)
    grad_mut(__context__, b)[:] = b̄
    return
  end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@test convert(Tuple, (1,2)) == (1,2)
@testset "indexing" begin
    @test length(()) == 0
    @test length((1,)) === 1
    @test length((1,2)) === 2

    @test_throws ArgumentError Base.front(())
    @test_throws ArgumentError first(())

    @test endof(()) === 0
    @test endof((1,)) === 1
    @test endof((1,2)) === 2

    @test size((), 1) === 0
    @test size((1,), 1) === 1
    @test size((1,2), 1) === 2
end
# @test_throws ArgumentError size((), 2)
# @test_throws ArgumentError size((1,), 2)
# @test_throws ArgumentError size((1,2), 2)

@test getindex((1,), 1) === 1
@test getindex((1,2), 2) === 2
@test_throws BoundsError getindex((), 1)
@test_throws BoundsError getindex((1,2), 0)
@test_throws BoundsError getindex((1,2), -1)

@test getindex((1,), 1.0) === 1
@test getindex((1,2), 2.0) === 2
@test_throws BoundsError getindex((), 1.0)
@test_throws BoundsError getindex((1,2), 0.0)
@test_throws BoundsError getindex((1,2), -1.0)

@test getindex((5,6,7,8), [1,2,3]) === (5,6,7)
@test_throws BoundsError getindex((1,2), [3,4])

@test getindex((5,6,7,8), [true, false, false, true]) === (5,8)
@test_throws BoundsError getindex((5,6,7,8), [true, false, false, true, true])

@test getindex((5,6,7,8), []) === ()

## filling to specified length
@test @inferred(Base.fill_to_length((1,2,3), -1, Val{5})) == (1,2,3,-1,-1)
@test_throws ArgumentError Base.fill_to_length((1,2,3), -1, Val{2})

## iterating ##
@test start((1,2,3)) === 1

@test done((), 1)
@test !done((1,2,3), 3)
@test done((1,2,3), 4)

@test next((5,6,7), 1) === (5, 2)
@test next((5,6,7), 3) === (7, 4)
@test_throws BoundsError next((5,6,7), 0)
@test_throws BoundsError next((), 1)

@test collect(eachindex((2,5,"foo"))) == collect(1:3)
@test collect(eachindex((2,5,"foo"), (1,2,5,7))) == collect(1:4)

## eltype ##
@test eltype((1,2,3)) === Int
@test eltype((1.0,2.0,3.0)) <: AbstractFloat
@test eltype((true, false)) === Bool
@test eltype((1, 2.0, false)) === typejoin(Int, Float64, Bool)
@test eltype(()) === Union{}
@test eltype(Tuple{Int, Float64, Vararg{Bool}}) === typejoin(Int, Float64, Bool)
@test eltype(Tuple{Int, T, Vararg{Bool}} where T <: AbstractFloat) ===
    typejoin(Int, AbstractFloat, Bool)
@test eltype(Tuple{Int, Bool, Vararg{T}} where T <: AbstractFloat) ===
    typejoin(Int, AbstractFloat, Bool)
@test eltype(Union{Tuple{Int, Float64}, Tuple{Vararg{Bool}}}) === typejoin(Int, Float64, Bool)

begin
    local foo
    ## mapping ##
    foo() = 2
    foo(x) = 2x
    foo(x, y) = x + y
    foo(x, y, z) = x + y + z
    longtuple = ntuple(identity, 20)

    # 1 argument
    @test map(foo, ()) === ()
    @test map(foo, (1,)) === (2,)
    @test map(foo, (1,2)) === (2,4)
    @test map(foo, (1,2,3,4)) === (2,4,6,8)
    @test map(foo, longtuple) === ntuple(i->2i,20)

    # 2 arguments
    @test map(foo, (), ()) === ()
    @test map(foo, (1,), (1,)) === (2,)
    @test map(foo, (1,2), (1,2)) === (2,4)
    @test map(foo, (1,2,3,4), (1,2,3,4)) === (2,4,6,8)
    @test map(foo, longtuple, longtuple) === ntuple(i->2i,20)
    @test_throws BoundsError map(foo, (), (1,))
    @test_throws BoundsError map(foo, (1,), ())

    # n arguments
    @test map(foo, (), (), ()) === ()
    @test map(foo, (1,), (1,), (1,)) === (3,)
    @test map(foo, (1,2), (1,2), (1,2)) === (3,6)
    @test map(foo, (1,2,3,4), (1,2,3,4), (1,2,3,4)) === (3,6,9,12)
    @test map(foo, longtuple, longtuple, longtuple) === ntuple(i->3i,20)
    @test_throws BoundsError map(foo, (), (1,), (1,))
    @test_throws BoundsError map(foo, (1,), (1,), ())
end

## comparison ##
@test isequal((), ())
@test isequal((1,2,3), (1,2,3))
@test !isequal((1,2,3), (1,2,4))
@test !isequal((1,2,3), (1,2))

@test ==((), ())
@test ==((1,2,3), (1,2,3))
@test !==((1,2,3), (1,2,4))
@test !==((1,2,3), (1,2))

@test isless((1,2), (1,3))
@test isless((1,), (1,2))
@test !isless((1,2), (1,2))
@test !isless((2,1), (1,2))

## functions ##
@test isempty(())
@test !isempty((1,))

@test reverse(()) === ()
@test reverse((1,2,3)) === (3,2,1)

## specialized reduction ##
@test sum((1,2,3)) === 6

@test prod(()) === 1
@test prod((1,2,3)) === 6

@test all(()) === true
@test all((false,)) === false
@test all((true,)) === true
@test all((true, true)) === true
@test all((true, false)) === false
@test all((false, false)) === false

@test any(()) === false
@test any((true,)) === true
@test any((false,)) === false
@test any((true, true)) === true
@test any((true, false)) === true
@test any((false, false)) === false
@test any((false,false,false)) === false
@test any((true,false,false)) === true
@test any((false,true,false)) === true
@test any((false,false,true)) === true
@test any((true,true,false)) === true
@test any((true,false,true)) === true
@test any((true,true,false)) === true
@test any((true,true,true)) === true

@test @inferred(ntuple(abs2, Val{0})) == ()
@test @inferred(ntuple(abs2, Val{2})) == (1, 4)
@test @inferred(ntuple(abs2, Val{3})) == (1, 4, 9)
@test @inferred(ntuple(abs2, Val{4})) == (1, 4, 9, 16)
@test @inferred(ntuple(abs2, Val{5})) == (1, 4, 9, 16, 25)
@test @inferred(ntuple(abs2, Val{6})) == (1, 4, 9, 16, 25, 36)

# issue #12854
@test_throws TypeError ntuple(identity, Val{1:2})

for n = 0:20
    t = ntuple(identity, n)
    @test length(t) == n
    for i = 1:n
        @test t[i] == i
    end
end

# issue #19719
@test_throws BoundsError (1,2,3)[falses(4)]
@test_throws BoundsError (1,2,3)[[false,false,true,true]]
@test_throws BoundsError (1,2,3)[trues(2)]
@test_throws BoundsError (1,2,3)[falses(2)]
@test_throws BoundsError ()[[false]]
@test_throws BoundsError ()[[true]]

struct BitPerm_19352
    p::NTuple{8,UInt8}
    function BitPerm(p::NTuple{8,UInt8})
        sort(collect(p)) != collect(0:7) && error("$p is not a permutation of 0:7")
        new(p)
    end
    BitPerm_19352(b0,b1,b2,b3,b4,b5,b6,b7) = BitPerm((UInt8(b0),UInt8(b1),UInt8(b2),UInt8(b3),
        UInt8(b4),UInt8(b5),UInt8(b6),UInt8(b7)))
end

@testset "side effect in tuple constructor #19352" begin
    @test BitPerm_19352(0,2,4,6,1,3,5,7).p[2] == 0x02
end

# issue #15703
let
    struct A_15703{N}
        keys::NTuple{N, Int}
    end

    struct B_15703
        x::A_15703
    end

    function bug_15703(xs...)
        [x for x in xs]
    end

    function test_15703()
        s = (1,)
        a = A_15703(s)
        ss = B_15703(a).x.keys
        @test ss === s
        bug_15703(ss...)
    end

    test_15703()
end

# PR #15516
@test Tuple{Char,Char}("za") === ('z','a')
@test_throws ArgumentError Tuple{Char,Char}("z")

@test NTuple{20,Int}(Iterators.countfrom(2)) === (2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21)
@test NTuple{20,Float64}(Iterators.countfrom(2)) === (2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,21.)
@test_throws ArgumentError NTuple{20,Int}([1,2])

@test Tuple{Vararg{Float32}}(Float64[1,2,3]) === (1.0f0, 2.0f0, 3.0f0)
@test Tuple{Int,Vararg{Float32}}(Float64[1,2,3]) === (1, 2.0f0, 3.0f0)
@test Tuple{Int,Vararg{Any}}(Float64[1,2,3]) === (1, 2.0, 3.0)
@test Tuple(ones(5)) === (1.0,1.0,1.0,1.0,1.0)
@test_throws MethodError convert(Tuple, ones(5))

@testset "Multidimensional indexing (issue #20453)" begin
    @test_throws MethodError (1,)[]
    @test_throws MethodError (1,1,1)[1,1]
end

@testset "ambiguity between tuple constructors #20990" begin
    Tuple16Int = Tuple{Int,Int,Int,Int,Int,Int,Int,Int,Int,Int,Int,Int,Int,Int,Int,Int}
    tuple16int = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
    @test Tuple16Int(tuple16int) isa Tuple16Int
end

# PR #21446
for n = 0:15
    @test ntuple(identity, Val{n}) == ntuple(identity, n)
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
using MPSKit,TensorKit

mpo = nonsym_ising_mpo();
state = InfiniteMPS([ℂ^2],[ℂ^10]);
(state,pars,_) = leading_boundary(state,mpo,Vumps(tol_galerkin=1e-10));

(state,pars,_) = leading_boundary(state,mpo,PowerMethod(tol_galerkin=1e-12, maxiter=400));


text:
"""
    doSom(train::Array, distMatrix, kernelFun, len, η,
            r, toroidal, rDecay, ηDecay)

Train a SOM.

This worker function is called by the high-level-API-functions
`som(), somHexagonal() and somSpherical()`.

# Arguments:
- `x`: training Data
- `dm`: distance matrix of all neurons of the SOM
- `kernelFun`: distance kernel function of type fun(x, r)
- `len`: number of training steps (*not* epochs)
- `η`: learning rate
- `r`: training radius
- `toroidal`: if true, the SOM is toroidal.
- `rDecay`: if true, r decays to 0.0 during the training.
- `ηDecay`: if true, learning rate η decays to 0.0 during the training.
"""
function doSom(x::Array{Float64}, codes::Array{Float64},
             dm::Array{Float64}, kernelFun::Function, len::Int, η::Float64,
             r::Number, toroidal::Bool, rDecay::Bool, ηDecay::Bool)

    # make Δs for linear decay of r and η:
    r = Float64(r)
    if rDecay
        if r < 1.5
            Δr = 0.0
        else
            Δr = (r-1.0) / len
        end
    else
        Δr = 0.0
    end

    if ηDecay
        Δη = η / len
    else
        Δη = 0.0
    end

    numDat = nrow(x)
    numCodes = nrow(codes)

    # Training:
    # 1) select random sample
    # 2) find winner
    # 3) train all neurons with gaussian kernel
    p = Progress(len, dt=1.0, desc="Training...", barglyphs=BarGlyphs("[=> ]"),
                 barlen=50, color=:yellow)
    @time for s in 1:len

        sampl = rowSample(x)
        winner = findWinner(codes, sampl)

        for i in 1:numCodes
            # v = view(codes, i, :)
            Δi = codes[i,:] .- sampl
            codes[i,:] -= Δi .* kernelFun(dm[winner,i], r) .* η
            # v -=  @. v - sampl * kernelFun(dm[winner,i], r) * η
        end

        η -= Δη
        r -= Δr
        next!(p)
    end

    return codes
end

"""
    function visual(codes, x)

Return the index of the winner neuron for each training pattern
in x (row-wise).
"""
function visual(codes, x)

    vis = zeros(Int, nrow(x))
    for i in 1:nrow(x)
        vis[i] = findWinner(codes, [x[i, col] for col in 1:size(x, 2)])
    end

    return(vis)
end

"""
    makePopulation(nCodes, vis)

Return a vector of neuron populations.
"""
function makePopulation(nCodes, vis)

    population = zeros(Int, nCodes)
    for i in 1:nrow(vis)
        population[vis[i]] += 1
    end

    return population
end

"""
    makeClassFreqs(som, vis, classes)

Return a DataFrame with class frequencies for all neurons.
"""
function makeClassFreqs(som, vis, classes)

    # count classes and construct DataFrame:
    #
    classLabels = sort(unique(classes))
    classNum = nrow(classLabels)

    cfs = DataFrame(index = 1:som.nCodes)
    cfs[!,:X] = som.indices.X
    cfs[!,:Y] = som.indices.Y

    cfs.Population = zeros(Int, som.nCodes)

    for class in classLabels
        cfs[!,Symbol(class)] = zeros(Float64, som.nCodes)
    end

    # loop vis and count:
    #
    for i in 1:nrow(vis)

        cfs[vis[i], :Population] += 1
        class = Symbol(classes[i])
        cfs[vis[i], class] += 1
    end

    # make frequencies from counts:
    #
    for i in 1:nrow(cfs)

        counts = [cfs[i, col] for col in 5:size(cfs, 2)]
        total = cfs[i,:Population]
        if total == 0
            freqs = counts * 0.0
        else
            freqs = counts ./ total
        end

        for c in 1:classNum
            class = Symbol(classLabels[c])
            cfs[i,class] = freqs[c]
        end
    end

    return cfs
end

"""
    function trainAll(om::SOM, rain::Array{Float64},
                len, η, kernelFun, r,
                rDecay, ηDecay)

Connects the high-level-API functions with
the backend.

# If `init == true` a new som is initialised with randomly sampled
# samples from train.
# Otherwise the som is trained.
"""
function trainAll(som::Som, train::Array{Float64,2},
                len, η, kernelFun, r,
                rDecay, ηDecay)

    # normalise training data:
    #
    if som.norm != :none
        train = normTrainData(train, som.normParams)
    end

    # set default radius:
    if r < 0.1
        if som.topol != :spherical
            r = √(som.xdim^2 + som.ydim^2) / 2
        else
            r = π * som.ydim
        end
    end

    if som.topol == :spherical
        dm = distMatrixSphere(som.grid)
    else
        dm = distMatrix(som.grid, som.toroidal)
    end

    # println("$(show(IOContext(STDOUT, limit=true), "text/plain", train))")
    codes = doSom(train, som.codes, dm,
                  kernelFun, len, η, r,
                  som.toroidal, rDecay, ηDecay)

    # map training samles to SOM and calc. neuron population:
    vis = visual(codes, train)
    population = makePopulation(som.nCodes, vis)

    # create X,Y-indices for neurons:
    #
    if som.topol != :spherical
        x = [mod(i-1, som.xdim)+1 for i in 1:som.nCodes]
        y = [div(i-1, som.xdim)+1 for i in 1:som.nCodes]
    else
        x = y = collect(1:som.nCodes)
    end
    indices = DataFrame(X = x, Y = y)

    # update SOM object:
    somNew = deepcopy(som)
    somNew.codes[:,:] = codes[:,:]
    somNew.population[:] = population[:]
    return somNew
end

"""
    initAll( train::Array{Float64,2}, colNames::Array{String,1},
                norm::Symbol,
                xdim::Int, ydim::Int, nCodes::Int,
                topology::Symbol, toroidal::Bool)

Initialise a new Self-Organising Map.
"""
function initAll( train::Array{Float64,2}, colNames::Array{String,1},
                norm::Symbol,
                xdim::Int, ydim::Int, nCodes::Int,
                topology::Symbol, toroidal::Bool)

    # normalise training data:
    train, normParams = normDefine(train, norm)
    codes = initCodes(nCodes, train, colNames)

    if topology == :rectangular
        grid = gridRectangular(xdim, ydim)
    elseif topology == :hexagonal
        grid = gridHexagonal(xdim, ydim)
    elseif topology == :spherical
        grid = gridSpherical(nCodes)
    else
        error("Topology $topology is not supported!")
    end

    normParams = convert(DataFrame, normParams)
    names!(normParams, Symbol.(colNames))

    # create X,y-indices for neurons:
    #
    if topology != :spherical
        x = [mod(i-1, xdim)+1 for i in 1:nCodes]
        y = [div(i-1, xdim)+1 for i in 1:nCodes]
    else
        x = y = collect(1:nCodes)
    end
    indices = DataFrame(X = x, Y = y)

    # make SOM object:
    som = Som(codes = codes, colNames = colNames,
              normParams = normParams, norm = norm,
              xdim = xdim, ydim = ydim,
              nCodes = nCodes,
              grid = grid, indices = indices,
              topol = topology,
              toroidal = toroidal,
              population = zeros(Int, nCodes))
    return som
end


text:
export Poisson, poisson_aux
export Gauss, gauss_aux
export ScaledGauss, scaled_gauss_aux
export Anscombe, anscombe_aux

"""
    poisson_aux(μ, meas, storage=similar(μ))

Calculates the Poisson loss for `μ` and `meas`.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
"""
function poisson_aux(μ, meas, storage=similar(μ))
    # due to numerical errors, μ can be negative or 0
    if minimum(μ) <= 0
        μ .= μ .+ eps(maximum(μ)) .+ abs.(minimum(μ))
    end
    storage .= μ .- meas .* log.(μ)
    return sum(storage)
end

 # define custom gradient for speed-up
 # ChainRulesCore offers the possibility to define a backward AD rule
 # which can be used by several different AD systems
function ChainRulesCore.rrule(::typeof(poisson_aux), μ, meas, storage)
    Y = poisson_aux(μ, meas, storage)

    function poisson_aux_pullback(xbar)
        storage .= xbar .* (one(eltype(μ)) .- meas ./ μ)
        return NoTangent(), storage, (ChainRulesCore.@not_implemented "Save computation"), (ChainRulesCore.@not_implemented "Save computation")
    end

    return Y, poisson_aux_pullback
end

"""
    Poisson()

Returns a function to calculate Poisson loss
Check the help of `poisson_aux`.
"""
function Poisson()
      return poisson_aux
end

"""
    gauss_aux(μ, meas, storage=similar(μ))

Calculates the Gauss loss for `μ` and `meas`.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
"""
function gauss_aux(μ, meas, storage=similar(μ))
    storage .= abs2.(μ - meas)
    return sum(storage)
end

 # define custom gradient for speed-up
function ChainRulesCore.rrule(::typeof(gauss_aux), μ, meas, storage)
    Y = gauss_aux(μ, meas)
    function gauss_aux_pullback(xbar)
        return NoTangent(), 2 .* xbar .* (μ - meas), (ChainRulesCore.@not_implemented "Save computation"), (ChainRulesCore.@not_implemented "Save computation")
    end
    return Y, gauss_aux_pullback
end

"""
    Gauss()

Returns a function to calculate Gauss loss.
Check the help of `gauss_aux`.
"""
function Gauss()
    return gauss_aux
end

"""
    scaled_gauss_aux(μ, meas, storage=similar(μ); read_var=0)
Calculates the scaled Gauss loss for `μ` and `meas`.
`read_var=0` is the readout noise variance of the sensor.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
"""
function scaled_gauss_aux(μ, meas, storage=similar(μ); read_var=0)
    μ[μ .<= 1f-8] .= 1f-8
    storage .= log.(μ .+ read_var) .+ (meas .- μ).^2 ./ ((μ .+ read_var))
    return sum(storage)
end

 # define custom gradient for speed-up
function ChainRulesCore.rrule(::typeof(scaled_gauss_aux), μ, meas, storage; read_var=0)
    Y = scaled_gauss_aux(μ, meas)
    function scaled_gauss_aux_pullback(xbar)
        ∇ = xbar .* (μ.^2 .- meas.^2 .+ μ .+ read_var.*(1 .- 2 .* (meas .- µ)))./((μ .+read_var).^2)
        ∇[μ .<= 1f-8] .= 0
        return NoTangent(), ∇, (ChainRulesCore.@not_implemented "Save computation"), (ChainRulesCore.@not_implemented "Save computation"), (ChainRulesCore.@not_implemented "Save computation")
    end
    return Y, scaled_gauss_aux_pullback
end

"""
    ScaledGauss()

Returns a function to calculate scaled Gauss loss.
Check the help of `scaled_gauss_aux`.
"""
function ScaledGauss(read_var=0)
    return (µ, meas, storage=similar(µ)) -> scaled_gauss_aux(µ, meas, storage, read_var=read_var)
end

"""
    anscombe_aux(μ, meas, storage=similar(μ); b=0)

Calculates the Poisson loss using the Anscombe-based norm for `μ` and `meas`.
`μ` can be of larger size than `meas`. In that case
we extract a centered region from `μ` of the same size as `meas`.
`b=0` is the optional parameter under the `√`.
Note that the data will be normalized to the mean of the data, which means that you have to
devide this parameter also by the mean of the data, i.e. b=3.0/8.0/mean(measured).
"""
function anscombe_aux(μ, meas, storage=similar(μ); b=0)
    # we cannot devide b here by meas, since meas is already normalized
    # due to numerical errors, μ can be negative or 0
    mm = minimum(μ)
    if mm <= 0
        μ .= μ .+ eps(maximum(μ)) .+ abs.(mm)
    end
    storage .= abs2.(sqrt.(meas .+ b) .- sqrt.(μ .+ b))
    return sum(storage)
end

 # define custom gradient for speed-up
 # ChainRulesCore offers the possibility to define a backward AD rule
 # which can be used by several different AD systems
function ChainRulesCore.rrule(::typeof(anscombe_aux), μ, meas, storage; b=1)
    Y = anscombe_aux(μ, meas, storage, b=b)
    function anscombe_aux_pullback(xbar)
            storage .= xbar .* (one(eltype(μ)) .- sqrt.((meas .+ b) ./ (μ.+b)))
        return NoTangent(), storage, (ChainRulesCore.@not_implemented "Save computation"), (ChainRulesCore.@not_implemented "Save computation"), (ChainRulesCore.@not_implemented "Save computation")
    end

    return Y, anscombe_aux_pullback
end

"""
    Anscombe(b=0)
Returns a function to calculate Poisson loss using the Anscombe transform
Check the help of `anscombe_aux`.
"""

function Anscombe(b=0)
    (μ, meas, storage=similar(μ)) -> anscombe_aux(μ, meas, storage, b=b)
end


text:
using MIToS.MSA
using MIToS.Information

const msa_long = read("../data/PF00089_aligned.fasta", FASTA)
const msa_wide = read("../data/PF16957_aligned.fasta", FASTA)

mip(msa) = APC!(mapcolpairfreq!(mutual_information, msa, Counts{Float64,2,GappedAlphabet}(ContingencyTable(Float64,Val{2},GappedAlphabet()))))
mip(msa_long)
mip(msa_wide)

println("[BENCH] MIp PF00089: ", @elapsed mip(msa_long))
println("[BENCH] MIp PF16957: ", @elapsed mip(msa_wide))

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
using SMC
using Test

@testset "kalman" begin
    include("kalman.jl")
end

@testset "resampler" begin
    include("resampler.jl")
end

@testset "smc" begin
    include("smc.jl")
end

@testset "smc2_inference" begin
    include("smc2_inference.jl")
end


text:
using Caesar, RoME, IncrementalInference
using DelimitedFiles

# Visualization
using Colors, Gadfly
using RoMEPlotting, KernelDensityEstimatePlotting

#posFile = joinpath(dirname(pathof(Caesar)),"..","test","testdata","test_array_positions.csv");
posFile = joinpath(ENV["HOME"],"data","sas","test_array_positions.csv");
posData = readdlm(posFile,',',Float64,'\n')

posData[:,1] .+= 13.0
posData[:,2] .-= 40.0

# gps covariance
rtkCov = Matrix(Diagonal([0.1;0.1].^2))

## build the factor graph
N = 100
fg = initfg()

#landmark
addVariable!(fg, :l1, Point2)
# setVal!(fg, :l1, zeros(2,N))

global pp = nothing
for i in 1:5
    sym = Symbol("x$i")
    addVariable!(fg, sym, Point2)
    setVal!(fg, sym, zeros(2,N))
    pp = PriorPoint2(MvNormal(posData[i,1:2], rtkCov) )
    addFactor!(fg, [sym;], pp)
    # manual init for now
    #setValKDE!(getVariable(fg,sym), kde!(getSample(pp, N)[1]))
    #getData(fg,sym).initialized = true
end

#dataFile = joinpath(dirname(pathof(Caesar)),"..","test","testdata","test_array_waveforms.csv");
dataFile = joinpath(ENV["HOME"],"data","sas","test_array_waveforms.csv");

csvWaveData = readdlm(dataFile,',',Float64,'\n')
csvWaveData = csvWaveData[:,1:5]
#
cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml");
cfgd=loadConfigFile(cfgFile)
chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");
sas2d = prepareSAS2DFactor(5, csvWaveData, rangemodel=:Correlator, cfgd=cfgd, chirpFile=chirpFile)

#
addFactor!(fg, [:l1;:x1;:x2;:x3;:x4;:x5], sas2d, autoinit=false)
#, threadmodel=SingleThreaded )

writeGraphPdf(fg, engine="fdp")

# sas2d = getfnctype(getVert(fg, :l1x1x2x3x4x5f1, nt=:fnc)) # deepcopy somewhere in IIF

# writeGraphPdf(fg)
# ls(fg, :l1)
# getVal(fg, :l1)

getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = true

tree = solveTree!(fg)

drawTree(tree, show=true)

# forward SAS case
# for ti in 1:Threads.nthreads() reset!(sas2d.threadreuse[ti].dbg) end

@time predL1 = IncrementalInference.approxConv(fg, ls(fg, :l1)[1], :l1, N=N)
# setVal!(fg, :l1, zeros(2,N))

pllim = 800
X1 = getVal(fg, :x1)
pl = Gadfly.plot(
  Gadfly.layer(x=predL1[1,:][-pllim .< predL1[1,:] .< pllim],y=predL1[2,:][-pllim .< predL1[2,:] .< pllim],   Gadfly.Geom.hexbin(xbincount=100, ybincount=100)),
  Gadfly.layer(x=X1[1,:],y=X1[2,:], Gadfly.Geom.hexbin(xbincount=100, ybincount=100))
) #; # for skipping Juno display

pl |> PDF("/tmp/test.pdf"); @async run(`evince /tmp/test.pdf`)
pl |> SVG("/tmp/test.svg", 35cm, 25cm); @async run(`firefox /tmp/test.svg`)

## backwards case

# backwards leave one out case
setVal!(fg, :l1, 1.0*randn(2,N))
setVal!(fg, :x2, rand(MvNormal(posData[2,1:2], 1.0*Matrix{Float64}(Diagonal([1.0;1.0]))),N))
# for ti in 1:Threads.nthreads() reset!(sas2d.threadreuse[ti].dbg) end

@time predX2 = IncrementalInference.approxConv(fg, ls(fg, :l1)[1], :x2, N=N)

pl = Gadfly.plot(
layer(x=predX2[1,:], y=predX2[2,:], Geom.histogram2d(xbincount=500,ybincount=500)),
layer(x=posData[1:5,1],y=posData[1:5,2], Geom.point)
)

pl.coord = Coord.Cartesian(xmin=0,xmax=20,ymin=-50,ymax=-35)
pl |> SVG("/tmp/test.svg", 35cm, 25cm)

pl_wait = plotKDE(kde!(predX2))
pl.layers = [pl.layers; pl_wait.layers]

pl |> SVG("/tmp/test.svg", 35cm, 25cm)

0


text:
# Mostly taken wholesale from https://juliacomputing.com/blog/2017/02/21/finding-ioctls-with-cxx.html
using Cxx
CCompiler = Cxx.new_clang_instance(
    false #= don't julia definitions =#,
    true #= C mode (as opposed to C++) =#)

Cxx.cxxinclude(CCompiler, "/home/olie/Dropbox/working_dir/julia/offline-laptop/mykavosh/pynauty-0.6.0/nauty26r7/nauty.h") # Gives us an error but we can probably ignore it

PP = icxx"&$(Cxx.active_instances[2].CI)->getPreprocessor();"

P  = Cxx.active_instances[2].Parser

#= icxx"$PP->dumpMacroInfo($PP->getIdentifierInfo(\"WORDSIZE\"));" =#

getIdentifierInfo(PP, name) = icxx"$PP->getIdentifierInfo($name);"
getMacroInfo(PP, II::pcpp"clang::IdentifierInfo") = icxx"$PP->getMacroInfo($II);"
getMacroInfo(PP, name::String) = getMacroInfo(PP, getIdentifierInfo(PP, name))
tokens(MI::pcpp"clang::MacroInfo") = icxx"$MI->tokens();"

# Convert Tokens that are identifiers to strings, we'll use these later
tok_is_identifier(Tok) = icxx"$Tok.is(clang::tok::identifier);"
Base.String(II::pcpp"clang::IdentifierInfo") = unsafe_string(icxx"$II->getName().str();")
function Base.String(Tok::cxxt"clang::Token")
    @assert tok_is_identifier(Tok)
    II = icxx"$Tok.getIdentifierInfo();"
    @assert II != C_NULL
    String(II)
end
getSpelling(PP, Tok) = unsafe_string(icxx"$PP->getSpelling($Tok);")
function Base.show(io::IO, Tok::Union{cxxt"clang::Token",cxxt"clang::Token&"})
    print(io, unsafe_string(icxx"clang::tok::getTokenName($Tok.getKind());"))
    print(io, " '", getSpelling(PP, Tok), "'")
end

# Iteration for ArrayRef
import Base: start, next, length, done
const ArrayRef = cxxt"llvm::ArrayRef<$T>" where T
start(AR::ArrayRef) = 0
function next(AR::cxxt"llvm::ArrayRef<$T>", i) where T
    (icxx"""
        // Force a copy, otherwise we'll retain reference semantics in julia
        // which is not what people expect.
        $T element = ($AR)[$i];
        return element;
    """, i+1)
end
length(AR::ArrayRef) = icxx"$AR.size();"
done(AR::ArrayRef, i) = i >= length(AR)

# Usage stringFromMacro("WORDSIZE") etc
function stringFromMacro(macroname)
    join([getSpelling(PP,x) for x in tokens(getMacroInfo(PP,macroname))])
end


text:
function gencenters(side, λ)
    n = pois_rand(λ * side^3)
    return reduce(hcat, (rand(1:side, 3) for i in 1:n))
end

function genballs(side, R, λ)
    sphere_side = 2R + 1
    padded_side = side + sphere_side
    spheres = zeros(Int8, (padded_side, padded_side, padded_side))
    sphere = zeros(Int8, (sphere_side, sphere_side, sphere_side))
    centers = gencenters(side, λ)
    for i in -R:R
        for j in -R:R
            for k in -R:R
                dist = i^2 + j^2 + k^2
                if dist < R^2
                    sphere[k+R+1, j+R+1, i+R+1] = 1
                end
            end
        end
    end

    for center in (centers[:,i] for i in 1:size(centers,2))
        x = center[1]
        y = center[2]
        z = center[3]
        spheres[x:x + 2R, y:y + 2R, z:z + 2R] .|= sphere
    end
    return spheres[R+1:end-R-1, R+1:end-R-1, R+1:end-R-1]
end

heaviside(x) = max(sign(x), 0)

function s2_theory(r, R, λ)
    tmp = r/R
    tmp2 = (r > 2R) ? 2 : 1 + 3/4*tmp - 1/16*tmp^3
    η = λ * 4/3 * π * R^3
    return exp(-η*tmp2)
end

function ss_theory(r, R, λ)
    η = λ * 4/3 * π * R^3
    tmp = 3η/R
    tmp2 = (tmp * (1 - (0.5 - r/(4R))*heaviside(2R-r)))^2 + tmp/(2r)*heaviside(2R-r)
    return s2_theory(r, R, λ) * tmp2
end

function sv_theory(r, R, λ)
    η = λ * 4/3 * π * R^3
    tmp = 3η/R
    return s2_theory(r, R, λ) * tmp * (1 - (0.5 - r/(4R))*heaviside(2R-r))
end

function l2_theory(r, R, λ)
    η = λ * 4/3 * π * R^3
    p = exp(-η)
    return p^(1 + 3r/(4R))
end

function pore_size_theory(r, R, λ)
    s(x) = 4π * x^2
    v(x) = 4/3 * π * x^3
    η  = v(R) * λ
    p  = exp(-η)
    s1 = s(r + R)
    v1 = v(r + R)
    return λ * s1 * exp(-λ * v1) / p
end

function chord_length_theory(r, R, λ)
    η = λ * 4/3 * π * R^3
    p = exp(-η)
    return 3η/(4R)*p^(3r/(4R))
end

mean_chord_length(R, λ) = 1/(λ * π *R^2)

@testset "S2 on random overlapping balls generated by Poisson process" begin
    # Volume = S^3, radius of a ball = R
    # Poisson parameter = λ
    S = 500; R = 20; λ = 3e-5
    balls = genballs(S, R, λ)

    # Calculate in one direction for speed
    calc = Directional.s2(balls, 0; periodic = true) |> mean
    theory = [s2_theory(r - 1, R, λ) for r in 0:length(calc) - 1]
    err = relerr.(calc, theory)

    @test maximum(err) < 0.25
end

@testset "SS on random overlapping balls generated by Poisson process" begin
    # Volume = S^3, radius of a ball = R
    # Poisson parameter = λ
    S = 500; R = 20; λ = 3e-5;
    balls = genballs(S, R, λ)

    # Calculate in one direction for speed
    calc = Directional.surfsurf(balls, 0; periodic = true) |> mean
    theory = [ss_theory(r - 1, R, λ) for r in 0:length(calc) - 1]
    err = relerr.(calc, theory)

    # Surface-surface function with r < 2R has a discontinuity and two
    # points where it goes to +Inf. Do not test it.
    @test maximum(err[2R + 10:end]) < 0.25
end

@testset "SV on random overlapping balls generated by Poisson process" begin
    # Volume = S^3, radius of a ball = R
    # Poisson parameter = λ
    S = 500; R = 20; λ = 3e-5;
    balls = genballs(S, R, λ)

    # Calculate in one direction for speed
    calc = Directional.surfvoid(balls, 0; periodic = true) |> mean
    theory = [sv_theory(r - 1, R, λ) for r in 0:length(calc) - 1]
    err = relerr.(calc, theory)

    @test maximum(err) < 0.25
end

@testset "L2 on random overlapping balls generated by Poisson process" begin
    # Volume = S^3, radius of a ball = R
    # Poisson parameter = λ
    # Compare l2(x) for x ∈ {0,1,...,N-1}
    S = 500; R = 20; λ = 3e-5; N = 100
    balls = genballs(S, R, λ)

    calc = Directional.l2(balls, 0; len = N, periodic = true) |> mean .|> log
    theory = [log(l2_theory(r - 1, R, λ)) for r in 0:N - 1]
    err = relerr.(calc, theory)

    @test maximum(err) < 0.25
end

@testset "Pore size on random overlapping balls generated by Poisson process" begin
    # Volume = S^3, radius of a ball = R
    # Poisson parameter = λ
    S = 500; R = 30; λ = 1e-5
    balls = genballs(S, R, λ)

    calc = pore_size(balls; nbins = 20)
    edges = calc.edges[1]
    s = step(edges)
    theory = [integrate(x -> pore_size_theory(x, R, λ), n:0.05:n+s)
              for n in 0:s:s*(length(edges) - 2)]

    # Compare cummulative distributions instead of probability
    # densities because density is almost zero for big radii.
    calc_cdf = scan(calc.weights)
    theory_cdf = scan(theory)

    err = relerr.(calc_cdf, theory_cdf)
    @test maximum(err) < 0.25
end

@testset "Chord length on random overlapping balls generated by Poisson process" begin
    # Volume = S^3, radius of a ball = R
    # Poisson parameter = λ
    S = 500; R = 30; λ = 1e-5
    balls = genballs(S, R, λ)

    calc, mc = Directional.chord_length(balls, 0; nbins = 30)
    edges = calc.edges[1]
    s = step(edges)
    theory = [integrate(x -> chord_length_theory(x, R, λ), n:0.05:n+s)
              for n in 0:s:s*(length(edges) - 2)]

    # Compare cummulative distributions instead of probability
    # densities because density is almost zero for big radii.
    calc_cdf = scan(calc.weights)
    theory_cdf = scan(theory)

    err = relerr.(calc_cdf, theory_cdf)
    @test maximum(err) < 0.25
    @test relerr(mc, mean_chord_length(R, λ)) < 0.25
end