text:
abstract FixedArray{T, NDim, SIZE}
abstract MutableFixedArray{T, NDim, SIZE} <: FixedArray{T, NDim, SIZE}

typealias MutableFixedVector{T, CARDINALITY} MutableFixedArray{T, 1, Tuple{CARDINALITY}}
typealias MutableFixedMatrix{T, M, N}        MutableFixedArray{T, 2, Tuple{M,N}}

typealias FixedVector{CARDINALITY, T}        FixedArray{T, 1, Tuple{CARDINALITY,}}
typealias FixedMatrix{Row, Column, T}        FixedArray{T, 2, Tuple{Row, Column}}

abstract FixedVectorNoTuple{CARDINALITY, T} <: FixedVector{CARDINALITY, T}
export FixedVectorNoTuple

_size{T <: Tuple}(::Type{T}) = (T.parameters...)
_size{N, N2}(::Type{Tuple{N, N2}}) = (N,N2)
_size{N}(::Type{Tuple{N}}) = (N,)

eltype{T <: FixedArray}(A::Type{T}) = eltype_or(T, Any)
eltype{T <: FixedArray,N,SZ}(A::FixedArray{T,N,SZ}) = T

length{T <: FixedArray}(A::Type{T}) = prod(size(T))
length{T <: FixedArray}(A::T) = length(T)

endof{T <: FixedArray}(A::Type{T}) = length(T)
endof{T <: FixedArray}(A::T) = endof(T)

@generated function ndims{T <: FixedArray}(A::Type{T})
    :($(fsa_abstract(T).parameters[2]))
end
ndims{T <: FixedArray}(A::T) = ndims(T)

size{T,N,SZ}(A::Type{FixedArray{T,N,SZ}}) = _size(SZ)
size{T <: FixedArray}(A::Type{T}) = size_or(T, ())
size{T <: FixedArray}(A::T) = size(T)

size{T <: FixedArray}(A::Type{T}, d::Integer) = size(T)[d]
size{T <: FixedArray}(A::T, d::Integer) = size(T, d)

@generated function fsa_abstract{FSA <: FixedArray}(::Type{FSA})
    ff = FSA
    while ff.name.name != :FixedArray
       ff = supertype(ff)
    end
    :($ff)
end
@generated function size_or{FSA <: FixedArray}(::Type{FSA}, OR)
    fsatype = fsa_abstract(FSA)
    sz = fsatype.parameters[3]
    isa(sz, TypeVar) && return :(OR)
    any(x->isa(x, TypeVar), sz.parameters) && return :(OR)
    :($(_size(sz)))
end
@generated function eltype_or{FSA <: FixedArray}(::Type{FSA}, OR)
    fsatype = fsa_abstract(FSA)
    T = fsatype.parameters[1]
    isa(T, TypeVar) && return :(OR)
    :($T)
end
@generated function ndims_or{FSA <: FixedArray}(::Type{FSA}, OR)
    fsatype = fsa_abstract(FSA)
    N = fsatype.parameters[2]
    isa(N, TypeVar) && return :(OR)
    :($N)
end

# Iterator
start(A::FixedArray) = 1
function next(A::FixedArray, state::Integer)
    @inbounds x = A[state]
    (x, state+1)
end
done(A::FixedArray, state::Integer) = length(A) < state

@generated function basetype{T<:FixedArray}(::Type{T})
    :($(T.name.primary))
end

"""
    similar_type(::Type{FSA}, [::Type{T}=eltype(FSA)], [sz=size(FSA)])

Given an array type `FSA`, element type `T` and size `sz`, return a `FixedArray`
subtype which is as similar as possible.  `similar_type` is used in the same
spirit as `Base.similar` to store the results of `map()` operations, etc.
(`similar` cannot work here, because the types are generally immutable.)

By default, `similar_type` introspects `FSA` to determine whether `T` and `sz`
can be used; if not a canonical FixedArray container is returned instead.
"""
@pure function similar_type{T}(::Type{FixedArray}, ::Type{T}, sz::Tuple)
    if length(sz) == 1
        return Vec{sz[1],T}
    elseif length(sz) == 2
        return Mat{sz[1],sz[2],T}
    else
        throw(ArgumentError("No built in FixedArray type is implemented for eltype $T and size $sz"))
    end
end

@pure function similar_type{FSA <: FixedArray, T}(::Type{FSA}, ::Type{T}, sz::Tuple)
    fsa_size = fsa_abstract(FSA).parameters[3].parameters
    if eltype(FSA) == T && fsa_size == sz
        return FSA # Common case optimization
    end

    # The default implementation for similar_type is follows.  It involves a
    # fair bit of crazy type introspection: We check whether the type `FSA` has
    # the necessary type parameters to replace with `T` and `sz`, and if so
    # figure out how to do the replacement.  It's complicated because users may
    # arbitrarily rearrange type parameters in their subtypes, and possibly
    # even add new type parameters which aren't related to the abstract
    # FixedArray but should be preserved.

    # Propagate the available type parameters of FSA down to the abstract base
    # FixedArray as `TypeVar`s.
    pritype = FSA.name.primary
    fsatype = fsa_abstract(pritype)
    T_parameter    = fsatype.parameters[1]
    ndim_parameter = fsatype.parameters[2]
    sz_parameter   = fsatype.parameters[3]
    sz_parameters  = fsatype.parameters[3].parameters

    # Figure out whether FSA can accommodate the new eltype `T` and size `sz`.
    # If not, delegate to the fallback by default.
    if !((eltype(FSA) == T          || isa(T_parameter,    TypeVar)) &&
         (ndims(FSA)  == length(sz) || isa(ndim_parameter, TypeVar)) &&
         (fsa_size    == sz         || all(i -> (sz[i] == fsa_size[i] || isa(sz_parameters[i],TypeVar)), 1:length(sz))))
        return similar_type(FixedArray, T, sz)
    end

    # Iterate type parameters, replacing as necessary with T and sz
    params = collect(FSA.parameters)
    priparams = pritype.parameters
    for i=1:length(params)
        if priparams[i] === T_parameter
            params[i] = T
        elseif priparams[i] === ndim_parameter
            params[i] = length(sz)
        elseif priparams[i] === sz_parameter
            params[i] = Tuple{sz...}
        else
            for j = 1:length(sz_parameters)
                if priparams[i] === sz_parameters[j]
                    params[i] = sz[j]
                end
            end
        end
    end
    pritype{params...}
end

# similar_type versions with defaulted eltype and size
@pure similar_type{FSA <: FixedArray, T}(::Type{FSA}, ::Type{T}) = similar_type(FSA, T, size(FSA))
@pure similar_type{FSA <: FixedArray}(::Type{FSA}, sz::Tuple) = similar_type(FSA, eltype(FSA), sz)

# Deprecated similar() -> similar_type()
function get_tuple(f::FixedArray)
    Base.depwarn("get_tuple(f::FixedArray) is deprecated, use Tuple(f) instead", :get_tuple)
    Tuple(f)
end
function similar{FSA <: FixedArray}(::Type{FSA}, args...)
    Base.depwarn("similar{FSA<:FixedArray}(::Type{FSA}, ...) is deprecated, use similar_type instead", :similar)
    similar_type(FSA, args...)
end
similar{FSA <: FixedArray}(::Type{FSA}, sz::Int...) = similar(FSA, eltype(FSA), sz)
similar{FSA <: FixedArray,T}(::Type{FSA}, ::Type{T}, sz::Int...) = similar(FSA, T, sz)

@compat @generated function (::Type{T}){T<:Tuple, N, T1}(f::FixedVectorNoTuple{N, T1})
    return Expr(:tuple, ntuple(i->:(f[$i]), N)...)
end
@compat function (::Type{T}){T<:Tuple}(f::FixedArray)
    getfield(f, 1)
end

# Infrastructure for construct_similar
"Compute promoted element type of the potentially nested Tuple type `ty`"
function promote_type_nested(ty)
    if ty.name.primary === Tuple
        promote_type([promote_type_nested(t) for t in ty.parameters]...)
    else
        ty
    end
end

"""
Construct tuple expression converting inner elements of the nested Tuple type
`ty` with name `varexpr` to the scalar `outtype`
"""
function convert_nested_tuple_expr(outtype, varexpr, ty)
    if ty.name.primary === Tuple
        Expr(:tuple, [convert_nested_tuple_expr(outtype, :($varexpr[$i]), t)
                      for (i,t) in enumerate(ty.parameters)]...)
    else
        :(convert($outtype, $varexpr))
    end
end

"""
Compute the N-dimensional array shape of a nested Tuple if it were used as
column-major storage for a FixedArray.
"""
function nested_Tuple_shape(ty)
    if ty.name.primary !== Tuple
        return 0
    end
    subshapes = [nested_Tuple_shape(t) for t in ty.parameters]
    if isempty(subshapes)
        return ()
    end
    if any(subshapes .!= subshapes[1])
        throw(DimensionMismatch("Nested tuples must have equal length to form a FixedSizeArray"))
    end
    if subshapes[1] == 0
        return (length(subshapes),) # Scalar elements
    end
    return (subshapes[1]..., length(subshapes))
end

"""
    construct_similar(::Type{FSA}, elements::Tuple)

Construct FixedArray as similar as possible to `FSA`, but with shape given by
the shape of the nested tuple `elements` and with `eltype` equal to the
promoted element type of the nested tuple `elements`.
"""
@generated function construct_similar{FSA <: FixedArray}(::Type{FSA}, elements::Tuple)
    etype = promote_type_nested(elements)
    shape = nested_Tuple_shape(elements)
    outtype = similar_type(FSA, etype, shape)
    converted_elements = convert_nested_tuple_expr(etype, :elements, elements)
    constructor_expr(outtype, converted_elements)
end


text:
using Documenter, AutoARIMA

makedocs(
    sitename="AutoARIMA.jl"
)

deploydocs(
    devbranch = "main",
    repo = "github.com/pierrenodet/AutoARIMA.jl.git",
)

ext: jl


text:
# NB: this test won't run automatically - you will have to include it manually
using JustJoshing
using Test
using CuArrays

using Statistics: mean

@testset "CPU/GPU consistency" begin
    let S_t=100, t=0.01, r=0.02, σ=0.05, T=1.0, trials=100_000_000
        payoff = S->max(S-100,0)
        kwargs = (S_0=S_t,t=t,T=T,r=r,σ=σ,trials=trials)

        cpu_vanilla = mc_pricer(payoff;kwargs...) |> first
        gpu_vanilla = JustJoshing.CUDA.mc_pricer(payoff;kwargs...) |> first

        @test cpu_vanilla ≈ gpu_vanilla rtol=1e-2

        # GPU seems to always be greater than CPU by a small amount
        @test_broken cpu_vanilla ≈ gpu_vanilla rtol=1e-3

        kwargs = (S_0=S_t,t=0:0.1:1,trials=trials)
        cpu_pathdep = mc_pricer_pathdep(S->max(mean(S)-100,0);kwargs...) |> first
        gpu_pathdep = JustJoshing.CUDA.mc_pricer_pathdep(S->max.(mean(S;dims=2).-100,0);kwargs...) |> first

        # These numbers are very different
        @test_broken cpu_pathdep ≈ gpu_pathdep rtol=1e-3

    end
end


text:
#=
    Figure 4
=#

using Inference
using Model

using JLD2
using Plots
using StatsPlots
using Random

include("../figure_defaults.jl")

##############################################################
## Load model and results
##############################################################

include("../../Results/MainResult_Setup.jl")
@load "Results/MainResult.jld2"

##############################################################
## Compute results
##############################################################

# Simulate model
Y = simulate_model(θ,T,param_dist,E;N=100_000,η)

##############################################################
## Figure 4(a,b)
##############################################################

fig4ab = plot_fit_univariate(X,Y,T)

##############################################################
## Figure 4(c,d)
##############################################################

fig4cd = plot_fit_bivariate(X,Y,T)

##############################################################
## Save...
##############################################################

savefig(fig4ab,"$(@__DIR__)/fig4ab.svg")
savefig(fig4cd,"$(@__DIR__)/fig4cd.svg")

max_stars_repo_name: ap-browning/internalisation

max_issues_repo_name: ap-browning/internalisation

max_forks_repo_name: ap-browning/internalisation


text:
"""
2D Predator-Prey Model

Taken from:
[14.3 pp2: A 2D Predator-Prey Model][auto07p/14.3]
of AUTO-07p manual.

[auto07p/14.3]: http://www.dam.brown.edu/people/sandsted/auto/auto07p.pdf#page=138
"""
module PredatorPrey

using DiffEqBase: ODEProblem
using StaticArrays: SVector
using Setfield: @lens
import Setfield

using ...Bifurcations: BifurcationProblem
using ...Codim2: cast_container, DiffEqCodim2Problem

f1(u, p) =
    p[2] * u[1] * (1 - u[1]) - u[1] * u[2] - p[1] * (1 - exp(- p[3] * u[1]))
f2(u, p) =
    - u[2] + p[4] * u[1] * u[2]

f(u::SVector, p, t) = SVector(f1(u, p), f2(u, p))

function f(du, u, p, t)
    du[1] = f1(u, p)
    du[2] = f2(u, p)
    nothing
end

const _f = f

function make_prob(
        p = (0.0, 3, 5, 3);
        f = _f,
        u0 = SVector(0.0, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis = (@lens _[1]),
        t_domain = (0.0, 1.0),
        kwargs...)
    x_min = SVector(-0.2, -Inf)
    x_max = SVector(1.0, Inf)
    x_min = cast_container(typeof(ode.u0), x_min)
    x_max = cast_container(typeof(ode.u0), x_max)
    return BifurcationProblem(
        ode, param_axis, t_domain;
        phase_space = (x_min, x_max),
        kwargs...)
end

prob = make_prob()
ode = prob.p.de_prob
p = ode.p
u0 = ode.u0
param_axis = prob.p.param_axis

function make_codim2_prob(
        p = ode.p;
        u0 = SVector(0.0, 0.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem{!(u0 isa SVector)}(f, u0, tspan, p),
        param_axis1 = (@lens _[1]),
        param_axis2 = (@lens _[3]),
        t_domain = ([-Inf, -Inf], [Inf, Inf]),
        kwargs...)
    return DiffEqCodim2Problem(
        ode,
        param_axis1,
        param_axis2,
        t_domain;
        kwargs...)
end

end  # module


text:
import Dates
import LibGit2
import Pkg
import UUIDs

if length(ARGS) != 4
    throw(
        ArgumentError(
            string(
                "Syntax: julia make.jl ",
                "\"REGISTRY_NAME\" \"REGISTRY_UUID\" ",
                "\"GIT_USER_NAME\" \"GIT_USER_EMAIL\"",
                )
            )
        )
end

const MY_REGISTRY_NAME = convert(String, strip(ARGS[1]))
const MY_REGISTRY_UUID = convert(String, strip(ARGS[2]))
const MY_REGISTRY_UUID_OBJECT = UUIDs.UUID(MY_REGISTRY_UUID)
const MY_GIT_USER_NAME = convert(String, strip(ARGS[3]))
const MY_GIT_USER_EMAIL = convert(String, strip(ARGS[4]))

original_directory = pwd()
project_root = joinpath(splitpath(@__DIR__)...)
cd(project_root)

function simple_retry_string(f::Function)::String
    need_to_continue::Bool = true
    f_result::String = ""
    while need_to_continue
        try
            f_result = f()
            need_to_continue = false
        catch exception
            need_to_continue = true
            @error("$(exception)", exception=exception,)
        end
    end
    return f_result
end

rm(
    joinpath(project_root, "build", "STARTED",);
    force = true,
    recursive = true,
    )
mkpath(joinpath(project_root, "build",))
open(joinpath(project_root, "build", "STARTED",), "w",) do f
    write(f, "$(repr(Dates.now()))\n",)
end

@info("loading utils.jl...")
include(joinpath(project_root,"utils.jl",))
@info("successfully loaded utils.jl")

@info("parsing `offline.toml`...")
configuration = Pkg.TOML.parsefile("offline.toml")
@info("successfully parsed `offline.toml`...")

@info("parsing `broken.packages.toml`...")
broken_packages = Pkg.TOML.parsefile("broken.packages.toml")
broken_packages_name = broken_packages["broken"]["packages"]["name"]
@info("successfully parsed `broken.packages.toml`...")

@info("downloading Project.toml and Manifest.toml files...")
projects_downloads = String[]
if haskey(configuration["include"], "projects")
    for project_toml_table in configuration["include"]["projects"]
        url = project_toml_table["url"]
        @info("downloading $(repr(url))")
        push!(
            projects_downloads,
            simple_retry_string(() -> Base.download(url)),
            )
    end
end
manifests_downloads = String[]
if haskey(configuration["include"], "manifests")
    for project_toml_table in configuration["include"]["manifests"]
        url = project_toml_table["url"]
        @info("downloading $(repr(url))")
        push!(
            manifests_downloads,
            simple_retry_string(() -> Base.download(url)),
            )
    end
end
@info("successfully downloaded all Project.toml and Manifest.toml files")

@info("cloning registries...")
registries_clones = String[]
packages_from_registries_to_download_all_packages = String[]
for registry_toml_table in configuration["include"]["registries"]
    url = registry_toml_table["url"]
    download_all_packages = get(
        registry_toml_table,
        "download_all_packages",
        false,
        )
    tmp = retry_function_until_success(
        () -> git_clone_registry(url);
        )
    push!(registries_clones, tmp,)
    registry_toml_file_path = joinpath(
        tmp,
        "Registry.toml",
        )
    registry_toml_file_parsed = Pkg.TOML.parsefile(
        registry_toml_file_path
        )
    registry_name = registry_toml_file_parsed["name"]
    registry_uuid = registry_toml_file_parsed["uuid"]
    if download_all_packages
        append!(
            packages_from_registries_to_download_all_packages,
            [x["name"] for x in
                collect(values(registry_toml_file_parsed["packages"]))],
            )
    end
    unique!(packages_from_registries_to_download_all_packages)
end
unique!(packages_from_registries_to_download_all_packages)
sort!(packages_from_registries_to_download_all_packages)
@info("successfully cloned all registries")

@info("processing the manifests of the given packages")
packages_to_manifest_process = String[]
names_from_combinations_to_build = String[]
pkgname_to_branchestobuild = Dict{String, Vector{String}}()
combinations_to_build = configuration["build"]["packages"]
for combination in combinations_to_build
    for package in combination
        if "name" in keys(package.values)
            name = convert(String, strip(package["name"]))
            push!(names_from_combinations_to_build, name)
            if "branch" in keys(package.values)
                branch = convert(String, strip(package["branch"]))
                if !haskey(pkgname_to_branchestobuild, name)
                    pkgname_to_branchestobuild[name] = String[]
                end
                push!(pkgname_to_branchestobuild[name], branch)
                unique!(pkgname_to_branchestobuild[name])
            end
        end
    end
end
unique!(names_from_combinations_to_build)
append!(
    packages_to_manifest_process,
    names_from_combinations_to_build,
    )
append!(
    packages_to_manifest_process,
    strip.(configuration["warmup"]["packages"]),
    )
for file in manifests_downloads
    append!(
        packages_to_manifest_process,
        sort(
            unique(strip.(collect(keys(Pkg.TOML.parsefile(file)))))
            ),
        )
end
for file in projects_downloads
    append!(
        packages_to_manifest_process,
        sort(
            unique(strip.(collect(keys(Pkg.TOML.parsefile(file)["deps"]))))
            ),
        )
end
my_own_project_toml_file = joinpath(
    project_root,
    "Project.toml",
    )
append!(
    packages_to_manifest_process,
    sort(
        unique(
            strip.(
                collect(
                    keys(
                        Pkg.TOML.parsefile(my_own_project_toml_file)["deps"]
                        )
                    )
                )
            )
        ),
    )
unique!(packages_to_manifest_process)
sort!(packages_to_manifest_process)
original_depot_path = [x for x in Base.DEPOT_PATH]
my_depot = joinpath(project_root, "depot",)
my_environment = joinpath(project_root, "environments", "temporary",)
rm(my_depot; force = true, recursive = true,)
rm(my_environment;force = true,recursive = true,)
empty!(Base.DEPOT_PATH)
pushfirst!(Base.DEPOT_PATH, my_depot,)
Pkg.activate(my_environment)
results_of_manifest_processing = String[]
for registry_toml_table in configuration["include"]["registries"]
    url = registry_toml_table["url"]
    Pkg.Registry.add(Pkg.RegistrySpec(url=url,))
end
Pkg.Registry.update()
for name in packages_to_manifest_process
    rm(
        joinpath(my_environment, "Project.toml",);
        force = true,
        recursive = true,
        )
    rm(
        joinpath(my_environment, "Manifest.toml",);
        force = true,
        recursive = true,
        )
    retry_function_until_success(
        () -> Pkg_add_name_ignore_julia_version_error(name);
        )
    environment_manifest_contents = Pkg.TOML.parsefile(
        joinpath(my_environment, "Manifest.toml",)
        )
    append!(
        results_of_manifest_processing,
        strip.(collect(keys(environment_manifest_contents))),
        )
end
rm(my_depot; force = true, recursive = true,)
rm(my_environment;force = true,recursive = true,)
unique!(results_of_manifest_processing)
sort!(results_of_manifest_processing)
@info("finished processing the manifests of the given packages")

@info("cloning package repositories...")
packages_to_clone = String[]
append!(
    packages_to_clone,
    strip.(packages_from_registries_to_download_all_packages),
    )
append!(
    packages_to_clone,
    strip.(packages_to_manifest_process),
    )
append!(
    packages_to_clone,
    strip.(results_of_manifest_processing),
    )
unique!(packages_to_clone)
sort!(packages_to_clone)
@debug(
    "Packages to clone: ",
    length(packages_to_clone),
    packages_to_clone,
    repr(packages_to_clone),
    )
rm(
    joinpath(project_root, "packages",);
    force = true,
    recursive = true,
    )
rm(
    joinpath(project_root, "repos",);
    force = true,
    recursive = true,
    )
registry_toml = Dict{String, Any}()
registry_toml["name"] = MY_REGISTRY_NAME
registry_toml["uuid"] = MY_REGISTRY_UUID
registry_toml["repo"] = project_root
registry_toml["description"] = """
This registry allows you to use Julia packages
 behind a firewall. For help, visit
 https://github.com/DilumAluthge/OfflineRegistry
"""
packages_section = get(registry_toml,"packages",Dict{String,Any}(),)
exclude = configuration["exclude"]["packages"]
append!(
    exclude,
    broken_packages["broken"]["packages"]["name"],
    )
n = length(registries_clones)
for i = 1:n
    registry_root = registries_clones[i]
    registry_file = joinpath(registry_root, "Registry.toml",)
    registry_config = Pkg.TOML.parsefile(registry_file)
    registry_packages = registry_config["packages"]
    registry_packages_uuids = collect(keys(registry_packages))
    p = length(registry_packages_uuids)
    for j = 1:p
        uuid = registry_packages_uuids[j]
        name = registry_packages[uuid]["name"]
        if !(name in exclude)
            old_package_path = registry_packages[uuid]["path"]
            new_package_path = joinpath(
                "packages",
                name[1:1],
                name,
                )
            new_repo_path = joinpath(
                "repos",
                name[1:1],
                string(name, ".jl.git",),
                )
            package_source = joinpath(registry_root, old_package_path,)
            package_destination = joinpath(project_root, new_package_path,)
            repo_destination = joinpath(project_root, new_repo_path,)
            mkpath(package_destination)
            package_toml = Pkg.TOML.parsefile(
                joinpath(package_source,"Package.toml",)
                )
            original_url = package_toml["repo"]
            package_toml["repo"] = repo_destination
            rm(
                joinpath(package_destination, "Package.toml",);
                force = true,
                recursive = true,
                )
            open(joinpath(package_destination, "Package.toml",), "w") do f
                Pkg.TOML.print(f, package_toml,)
            end
            cp(
                joinpath(package_source, "Compat.toml",),
                joinpath(package_destination, "Compat.toml",);
                force = true,
                )
            cp(
                joinpath(package_source, "Versions.toml",),
                joinpath(package_destination, "Versions.toml",);
                force = true,
                )
            try
                cp(
                    joinpath(package_source, "Deps.toml",),
                    joinpath(package_destination, "Deps.toml",);
                    force = true,
                    )
            catch e
                @debug(
                    "ignoring exception: ",
                    e,
                    joinpath(package_source, "Deps.toml",),
                    )
            end
            if name in packages_to_clone
                @debug("Cloning package $(j) of $(p)")
                mkpath(repo_destination)
                tmp = retry_function_until_success(
                    () -> git_clone_repo(original_url);
                    )
                cp(
                    tmp,
                    repo_destination;
                    force = true,
                    )
                if haskey(pkgname_to_branchestobuild, name)
                    branches = pkgname_to_branchestobuild[name]
                    repo = LibGit2.GitRepo(repo_destination)
                    try
                        all_remotes = LibGit2.remotes(repo)
                        if "origin" in all_remotes
                            remote_name = "origin"
                        else
                            remote_name = first(all_remotes)
                        end
                        remote = LibGit2.get(
                            LibGit2.GitRemote,
                            repo,
                            remote_name,
                            )
                        for branch in branches
                            branch_already_exists = isa(
                                LibGit2.lookup_branch(repo, branch),
                                LibGit2.GitReference,
                                )
                            if branch_already_exists
                                @debug(
                                    string(
                                        "Branch $(branch) already exists in ",
                                        "package $(name)",
                                        )
                                    )
                            else
                                @debug(
                                    string(
                                        "Checking out branch $(branch) for ",
                                        "package $(name)",
                                        )
                                    )
                                remote_tip = LibGit2.GitCommit(
                                    repo,
                                    "refs/remotes/$(remote_name)/$(branch)",
                                    )
                                LibGit2.branch!(
                                    repo,
                                    branch,
                                    string(LibGit2.GitHash(remote_tip));
                                    track = branch,
                                    )
                            end
                        end
                    catch e1
                        @warn("ignoring exception: ", e1,)
                    finally
                        try
                            close(remote)
                        catch e2
                            @warn("ignoring exception: ", e2,)
                        end
                    end
                end
            end
            packages_section[uuid] = Dict(
                "name" => name,
                "path" => new_package_path,
                )
        end
    end
end
@info("successfully cloned all package repositories")

@info("generating Registry.toml file...")
registry_toml["packages"] = packages_section
registry_toml_path = joinpath(project_root, "Registry.toml",)
open(registry_toml_path, "w") do f
    Pkg.TOML.print(f, registry_toml,)
end
@info("successfully wrote Registry.toml file to $(repr(registry_toml_path))")

@info("tracking files and staging changes...")
project_repo = LibGit2.GitRepo(project_root)
LibGit2.add!(project_repo, "Registry.toml",)
LibGit2.add!(project_repo, "packages",)
@info("successfully tracked all files and staged all changes")

@info("committing changes...")
commit_msg = "Automated commit made by https://github.com/DilumAluthge/OfflineRegistry on $(repr(Dates.now()))"
sig = LibGit2.Signature(
    MY_GIT_USER_NAME,
    MY_GIT_USER_EMAIL,
    )
LibGit2.commit(project_repo,commit_msg;author = sig,committer = sig,)
all_project_remotes = LibGit2.remotes(project_repo)
for project_remote in all_project_remotes
    LibGit2.remote_delete(project_repo, project_remote,)
end
@info("successfully committed all changes")

@info("adding packages to depot and building packages...")
rm(my_depot; force = true, recursive = true,)
rm(my_environment;force = true,recursive = true,)
empty!(Base.DEPOT_PATH)
pushfirst!(Base.DEPOT_PATH, my_depot,)
Pkg.activate(my_environment)
Pkg.Registry.add(Pkg.RegistrySpec(path=project_root,))
packages_to_warmup = sort(
    unique(
        strip.(configuration["warmup"]["packages"])
        )
    )
n = length(packages_to_clone)
for i = 1:n
    @debug("Building package $(i) of $(n)")
    name = packages_to_clone[i]
    rm(
        joinpath(my_environment, "Project.toml",);
        force = true,
        recursive = true,
        )
    rm(
        joinpath(my_environment, "Manifest.toml",);
        force = true,
        recursive = true,
        )
    retry_function_until_success(
        () -> Pkg_add_name_ignore_julia_version_error(name);
        )
    try
        Pkg.build(name; verbose = true,)
    catch e1
        @warn("ignoring exception: ", e1,)
    end
    for package_to_warmup in intersect(
            packages_to_warmup,
            keys(
                Pkg.TOML.parsefile(
                    joinpath(my_environment, "Manifest.toml",)
                    )
                ),
            )
        Pkg.add(package_to_warmup)
        try
            Pkg.build(package_to_warmup; verbose = true,)
        catch e2
            @warn("ignoring exception: ", e2,)
        end
        Base.eval(
            Main,
            Base.Meta.parse("import $(package_to_warmup)"),
            )
    end
end

combinations_to_build = configuration["build"]["packages"]
n = length(combinations_to_build)
for i = 1:n
    combination_i = combinations_to_build[i]
    num_packages_in_combination_i = length(combination_i)
    list_of_packagespecs = Vector{Pkg.Types.PackageSpec}(
        undef,
        num_packages_in_combination_i,
        )
    for j = 1:num_packages_in_combination_i
        combination_i_package_j = combination_i[j]
        my_kwargs_dict = Dict{Symbol, String}()
        if "name" in keys(combination_i_package_j.values)
            my_kwargs_dict[:name] = strip(
                combination_i_package_j["name"]
                )
        end
        if "branch" in keys(combination_i_package_j.values)
            my_kwargs_dict[:rev] = strip(
                combination_i_package_j["branch"]
                )
        end
        if "version" in keys(combination_i_package_j.values)
            my_kwargs_dict[:version] = strip(
                combination_i_package_j["version"]
                )
        end
        if length(my_kwargs_dict) > 0
            packagespec = Pkg.PackageSpec(; my_kwargs_dict...)
        end
        list_of_packagespecs[j] = packagespec
    end
    @debug(
        "Building combination $(i) of $(n)",
        combination,
        list_of_packagespecs,
        repr(combination),
        repr(list_of_packagespecs),
        )
    rm(
        joinpath(my_environment, "Project.toml",);
        force = true,
        recursive = true,
        )
    rm(
        joinpath(my_environment, "Manifest.toml",);
        force = true,
        recursive = true,
        )
    retry_function_until_success(
        () -> Pkg.add(list_of_packagespecs);
        )
    for ps in list_of_packagespecs
        try
            Pkg.build(ps; verbose = true,)
        catch e1
            @warn("ignoring exception: ", e1,)
        end
    end
    for package_to_warmup in intersect(
            packages_to_warmup,
            keys(
                Pkg.TOML.parsefile(
                    joinpath(my_environment, "Manifest.toml",)
                    )
                ),
            )
        Pkg.add(package_to_warmup)
        try
            Pkg.build(package_to_warmup; verbose = true,)
        catch e2
            @warn("ignoring exception: ", e2,)
        end
        Base.eval(
            Main,
            Base.Meta.parse("import $(package_to_warmup)"),
            )
    end
end
rm(my_environment;force = true,recursive = true,)

moved_out_of_depot_dir = joinpath(
    project_root,
    "movedoutofdepot",
    )
moved_out_of_depot_packages_dir = joinpath(
    moved_out_of_depot_dir,
    "packages",
    )
move_out_of_depot_list = configuration["move_out_of_depot"]["packages"]
unique!(move_out_of_depot_list)
sort!(move_out_of_depot_list)
n = length(move_out_of_depot_list)
for i = 1:n
    name = move_out_of_depot_list[i]
    @debug("Moving \"$(name)\" ($(i) of $(n)) out of depot")
    old_package_path = joinpath(
        my_depot,
        "packages",
        name,
        )
    new_package_path = joinpath(
        moved_out_of_depot_packages_dir,
        name,
        )
    rm(
        new_package_path;
        force = true,
        recursive = true,
        )
    mkpath(moved_out_of_depot_packages_dir)
    mv(
        old_package_path,
        new_package_path;
        force = true,
        )
end

rm(joinpath(my_depot, "compiled",);force = true,recursive = true,)
rm(my_environment;force = true,recursive = true,)

empty!(Base.DEPOT_PATH)
for x in original_depot_path
    push!(Base.DEPOT_PATH, x,)
end
unique!(Base.DEPOT_PATH)

@info("successfully added packages to depot and built packages")

rm(
    joinpath(project_root, "build", "FINISHED",);
    force = true,
    recursive = true,
    )
mkpath(joinpath(project_root, "build",))
open(joinpath(project_root, "build", "FINISHED",), "w",) do f
    write(f, "$(repr(Dates.now()))\n",)
end

cd(original_directory)

max_forks_repo_name: DilumAluthge/PredictMDRegistryOffline


text:
using VectorizationBase: Vec, MM, mask, pick_vector_width_shift, pick_vector_width, extractelement

using Base.Cartesian: @nexprs, @ntuple
using Base: IEEEFloat

import LinearAlgebra: lmul!, rmul!, dot, norm, axpy!

@generated function rmul!(A::StructArray{MultiFloat{T,N}, 2}, G::Rotation2{MultiFloat{T,N},MultiFloat{T,N}}) where {T<:IEEEFloat,N}

    W, Wshift = pick_vector_width_shift(T)

    quote
        @nexprs $N k -> p_k = stridedpointer(reinterpret(T, getfield(A, 1)[k]))

        m = size(A, 1)

        i = 1

        for step = Base.OneTo(m >> $Wshift)
            @nexprs $N k -> v_k_1 = vload(p_k, (MM{$W}(i), G.i))
            @nexprs $N k -> v_k_2 = vload(p_k, (MM{$W}(i), G.i+1))

            a₁ = MultiFloat(@ntuple $N k -> v_k_1)
            a₂ = MultiFloat(@ntuple $N k -> v_k_2)

            a₁′ = a₁ * Vec{$W,$T}(G.c) + a₂ * Vec{$W,$T}(G.s)
            a₂′ = a₁ * -Vec{$W,$T}(G.s) + a₂ * Vec{$W,$T}(G.c)

            @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (MM{$W}(i), G.i))
            @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (MM{$W}(i), G.i+1))

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return A
        remainder_mask = mask($T, remainder)

        # Finally do the masked load, multiply and store.

        @nexprs $N k -> v_k_1 = vload(p_k, (MM{$W}(i), G.i), remainder_mask)
        @nexprs $N k -> v_k_2 = vload(p_k, (MM{$W}(i), G.i+1), remainder_mask)

        a₁ = MultiFloat(@ntuple $N k -> v_k_1)
        a₂ = MultiFloat(@ntuple $N k -> v_k_2)

        a₁′ = a₁ * Vec{$W,$T}(G.c) + a₂ * Vec{$W,$T}(G.s)
        a₂′ = a₁ * -Vec{$W,$T}(G.s) + a₂ * Vec{$W,$T}(G.c)

        @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (MM{$W}(i), G.i), remainder_mask)
        @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (MM{$W}(i), G.i+1), remainder_mask)

        A
    end
end

@generated function lmul!(G::Rotation2{MultiFloat{T,N},MultiFloat{T,N}}, A::StructArray{MultiFloat{T,N}, 2}) where {T<:IEEEFloat,N}

    W, Wshift = pick_vector_width_shift(T)

    quote
        @nexprs $N k -> p_k = stridedpointer(reinterpret(T, getfield(A, 1)[k]))

        m = size(A, 2)

        i = 1

        for step = Base.OneTo(m >> $Wshift)
            @nexprs $N k -> v_k_1 = vload(p_k, (G.i, MM{$W}(i)))
            @nexprs $N k -> v_k_2 = vload(p_k, (G.i+1, MM{$W}(i)))

            a₁ = MultiFloat(@ntuple $N k -> v_k_1)
            a₂ = MultiFloat(@ntuple $N k -> v_k_2)

            a₁′ = Vec{$W,$T}(G.c) * a₁ + Vec{$W,$T}(G.s) * a₂
            a₂′ = -Vec{$W,$T}(G.s) * a₁ + Vec{$W,$T}(G.c) * a₂

            @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (G.i, MM{$W}(i)))
            @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (G.i+1, MM{$W}(i)))

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return A
        remainder_mask = mask($T, remainder)

        # Finally do the masked load, multiply and store.
        @nexprs $N k -> v_k_1 = vload(p_k, (G.i, MM{$W}(i)), remainder_mask)
        @nexprs $N k -> v_k_2 = vload(p_k, (G.i+1, MM{$W}(i)), remainder_mask)

        a₁ = MultiFloat(@ntuple $N k -> v_k_1)
        a₂ = MultiFloat(@ntuple $N k -> v_k_2)

        a₁′ = Vec{$W,$T}(G.c) * a₁ + Vec{$W,$T}(G.s) * a₂
        a₂′ = -Vec{$W,$T}(G.s) * a₁ + Vec{$W,$T}(G.c) * a₂

        @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (G.i, MM{$W}(i)), remainder_mask)
        @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (G.i+1, MM{$W}(i)), remainder_mask)

        A
    end
end

@generated function MultiFloatOfVec(fs::NTuple{M,MultiFloat{T,N}}) where {T,M,N}
    exprs = [:(Vec($([:(fs[$j]._limbs[$i]) for j=1:M]...))) for i=1:N]

    return quote
        $(Expr(:meta, :inline))
        MultiFloat(tuple($(exprs...)))
    end
end

@generated function TupleOfMultiFloat(fs::MultiFloat{Vec{M,T},N}) where {T,M,N}
    exprs = [:(MultiFloat(tuple($([:(extractelement(fs._limbs[$j], $i)) for j=1:N]...)))) for i=0:M-1]
    return quote
        $(Expr(:meta, :inline))
        tuple($(exprs...))
    end
end

@generated function LinearAlgebra.dot(xs::AbstractArray{MultiFloat{T,N}}, ys::AbstractArray{MultiFloat{T,N}}) where {T,N}
    M = pick_vector_width(T)

    load_xs = ntuple(k -> :(xs[i + $(k - 1)]), M)
    load_ys = ntuple(k -> :(ys[i + $(k - 1)]), M)

    quote
        vec_total = zero(MultiFloat{Vec{$M,$T},$N})

        # iterate in steps of M
        step, i = 1, 1
        @inbounds while step ≤ length(xs) ÷ $M
            x = MultiFloatOfVec(tuple($(load_xs...)))
            y = MultiFloatOfVec(tuple($(load_ys...)))
            vec_total += x * y
            step += 1
            i += $M
        end

        # sum the remainder
        @inbounds total = +(TupleOfMultiFloat(vec_total)...)
        while i ≤ length(xs)
            @inbounds total += xs[i] * ys[i]
            i += 1
        end

        return total
    end
end

@generated function LinearAlgebra.axpy!(a::MultiFloat{T,N}, xs::AbstractArray{MultiFloat{T,N}}, ys::AbstractArray{MultiFloat{T,N}}) where {T,N}
    M = pick_vector_width(T)

    load_xs = ntuple(k -> :(xs[i + $(k - 1)]), M)
    load_ys = ntuple(k -> :(ys[i + $(k - 1)]), M)
    assign = [:($(load_ys[i]) = t[$i]) for i = 1:M]
    broadcast_a = ntuple(i -> :a, M)

    quote
        # iterate in steps of M
        step, i = 1, 1
        av = MultiFloatOfVec(tuple($(broadcast_a...)))
        @inbounds while step ≤ length(xs) ÷ $M
            x = MultiFloatOfVec(tuple($(load_xs...)))
            y = MultiFloatOfVec(tuple($(load_ys...)))
            result = y + x * av
            t = TupleOfMultiFloat(result)
            $(assign...)
            step += 1
            i += $M
        end

        # sum the remainder
        while i ≤ length(xs)
            @inbounds ys[i] += xs[i] * a
            i += 1
        end

        return ys
    end
end

LinearAlgebra.norm(xs::AbstractArray{<:MultiFloat}) = sqrt(dot(xs, xs))


text:
@testset "triangular cylinder 3D" begin
    vxl = (SVector(0,0,0), SVector(1,1,1))
    nout = SVector(1,1,0)

    @test_nowarn (r₀ = SVector(0.5,0,0); @inferred(volfrac(vxl, nout, r₀)))
    @test (r₀ = SVector(0.5,0,0); volfrac(vxl, nout, r₀) ≈ 0.125)
    @test (r₀ = SVector(0.5,0,0); volfrac(vxl, -nout, r₀) ≈ 0.875)
    @test (r₀ = SVector(1,0,0); volfrac(vxl, nout, r₀) ≈ 0.5)
    @test (r₀ = SVector(1,0,0); volfrac(vxl, -nout, r₀) ≈ 0.5)
    @test (r₀ = SVector(1,0,0); nout = SVector(1,2,0); volfrac(vxl, nout, r₀) ≈ 0.25)
    @test (r₀ = SVector(1,0,0); nout = SVector(1,2,0); volfrac(vxl, -nout, r₀) ≈ 0.75)
end  # @testset "triangular cylinder 3D"

@testset "triangular cylinder 2D" begin
    vxl = (SVector(0,0), SVector(1,1))
    nout = SVector(1,1)

    @test_nowarn (r₀ = SVector(0.5,0); @inferred(volfrac(vxl, nout, r₀)))
    @test (r₀ = SVector(0.5,0); volfrac(vxl, nout, r₀) ≈ 0.125)
    @test (r₀ = SVector(0.5,0); volfrac(vxl, -nout, r₀) ≈ 0.875)
    @test (r₀ = SVector(1,0); volfrac(vxl, nout, r₀) ≈ 0.5)
    @test (r₀ = SVector(1,0); volfrac(vxl, -nout, r₀) ≈ 0.5)
    @test (r₀ = SVector(1,0); nout = SVector(1,2); volfrac(vxl, nout, r₀) ≈ 0.25)
    @test (r₀ = SVector(1,0); nout = SVector(1,2); volfrac(vxl, -nout, r₀) ≈ 0.75)
end  # @testset "triangular cylinder 2D"

@testset "quadrangular cylinder 3D" begin
    @test_nowarn @inferred(volfrac((SVector(0,0,0),SVector(1,1,1)), SVector(1,2,0), SVector(0.5,0.5,0.5)))
    @test begin
        result = true
        for i = 1:100
            vxl = (-@SVector(rand(3)), @SVector(rand(3)))
            r₀ = mean(vxl)
            nout = randn(3)
            nout[rand(1:3)] = 0
            result &= volfrac(vxl, SVector{3}(nout), r₀)≈0.5
        end
        result
    end
end  # @testset "quadrangular cylinder 3D"

@testset "quadrangular cylinder 2D" begin
    @test_nowarn @inferred(volfrac((SVector(0,0),SVector(1,1)), SVector(1,2), SVector(0.5,0.5)))
    @test begin
        result = true
        for i = 1:100
            vxl = (-@SVector(rand(2)), @SVector(rand(2)))
            r₀ = mean(vxl)
            nout = @SVector randn(2)
            result &= volfrac(vxl, nout, r₀)≈0.5
        end
        result
    end
end  # @testset "quadrangular cylinder 2D"

@testset "general cases" begin
    # Test random cases.
    @test_nowarn @inferred(volfrac((-@SVector(rand(3)),@SVector(rand(3))), SVector(0,0,0), @SVector(randn(3))))
    @test begin
        result = true
        for i = 1:100
            vxl = (-@SVector(rand(3)), @SVector(rand(3)))
            r₀ = mean(vxl)
            nout = @SVector randn(3)
            result &= volfrac(vxl, nout, r₀)≈0.5
        end
        result
    end

    @test begin
        result = true
        for i = 1:100
            vxl = (-@SVector(rand(3)), @SVector(rand(3)))
            r₀ = @SVector randn(3)
            nout = @SVector randn(3)
            result &= (volfrac(vxl, nout, r₀) + volfrac(vxl, -nout, r₀) ≈ 1)
        end
        result
    end

    # Test boundary cases.
    vxl = (SVector(0,0,0), SVector(1,1,1))
    r₀ = SVector(0,0,0)
    @test (nout = SVector(1,0,0); volfrac(vxl, nout, r₀) ≈ 0)  # completely outside
    @test (nout = SVector(-1,0,0); volfrac(vxl, nout, r₀) ≈ 1)  # completely inside
    @test (nout = SVector(1,1,0); volfrac(vxl, nout, r₀) ≈ 0)  # completely outside
    @test (nout = SVector(-1,-1,0); volfrac(vxl, nout, r₀) ≈ 1)  # completely inside
    @test (nout = SVector(1,1,1); volfrac(vxl, nout, r₀) ≈ 0)  # completely outside
    @test (nout = SVector(-1,-1,-1); volfrac(vxl, nout, r₀) ≈ 1)  # completely inside
    @test (nout = SVector(-1,1,0); volfrac(vxl, nout, r₀) ≈ 0.5)  # rvol_tricyl()
    @test (nout = SVector(-1,-1,1); volfrac(vxl, nout, r₀) ≈ 5/6)  # rvol_gensect()
    @test (nout = SVector(1,-2,1); volfrac(vxl, nout, r₀) ≈ 0.5)  # rvol_quadsect()
    r₀ = SVector(0.5,0.5,0)
    @test (nout = SVector(-2,1,0); volfrac(vxl, nout, r₀) ≈ 0.5)  # rvol_quadcyl()

    # Test tolerance to floating-point arithmetic.
    vxl = (SVector(0,9.5,0), SVector(1,10.5,1))
    r₀ = SVector(0.5,10.0,0.5)
    nout = SVector(0.0,1/√2,1/√2)
    @test volfrac(vxl, nout, r₀) ≈ 0.5

    # Test rvol_quadsect() for nontrivial cases.
    vxl = (SVector(0,0,0), SVector(1,1,2))
    r₀ = SVector(0.5, 0.5, 0.5)
    @test begin
        result = true
        for i = 1:100
            nout = SVector(randn()/20, randn()/20, 1)
            result &= volfrac(vxl, nout, r₀)≈0.5/2
        end
        result
    end
end  # @testset "general cases"


text:
using XGBoost

##
#  this script demonstrate how to fit generalized linear model in xgboost
#  basically, we are using linear model, instead of tree for our boosters
##

dtrain = DMatrix("../data/agaricus.txt.train")
dtest = DMatrix("../data/agaricus.txt.test")

# change booster to gblinear, so that we are fitting a linear model
# alpha is the L1 regularizer
# lambda is the L2 regularizer
# you can also set lambda_bias which is L2 regularizer on the bias term

param = ["booster"=>"gblinear", "eta"=>1, "silent"=>0,
         "objective"=>"binary:logistic", "alpha"=>0.0001, "lambda"=>1]
# normally, you do not need to set eta (step_size)
# XGBoost uses a parallel coordinate descent algorithm (shotgun),
# there could be affection on convergence with parallelization on certain cases
# setting eta to be smaller value, e.g 0.5 can make the optimization more stable

##
# the rest of settings are the same
#

watchlist  = [(dtest,"eval"), (dtrain,"train")]
num_round = 4

bst = xgboost(dtrain, num_round, param=param, watchlist=watchlist)

preds = predict(bst, dtest)
labels = get_info(dtest, "label")

print("test-error=", sum((preds .> 0.5) .!= labels) / float(size(preds)[1]), "\n")


text:
using DiffEqPhysics

function generate_bodies_in_cell_nodes(n::Int, m::Real, v_dev::Real, L::Real)

    rng = MersenneTwister(n);
    velocities = v_dev * randn(rng, Float64, (3, n))
    bodies = MassBody[]

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x = dL/2:dL:L, y = dL/2:dL:L, z = dL/2:dL:L
        if count > n
            break
        end
        r = SVector(x, y, z)
        v = SVector{3}(velocities[:,count])
        body = MassBody(r, v, m)
        push!(bodies, body)
        count += 1
    end
    return bodies
end

function generate_bodies_in_line(n::Int, m::Real, v_dev::Real, L::Real)
    dL = L / (ceil(n^(1 / 3)))
    n_line = floor(Int, L / dL)
    rng = MersenneTwister(n);
    velocities = v_dev * randn(rng, Float64, (3, n_line))
    bodies = MassBody[]
    x = y = L / 2
    for i = 1:n_line
        r = SVector(x, y, i * dL)
        v = SVector{3}(velocities[:,i])
        body = MassBody(r, v, m)
        push!(bodies, body)
    end
    return bodies
end

function generate_random_directions(n::Int)
    theta = acos.(1 - 2 * rand(n));
    phi = 2 * pi * rand(n);
    directions = [@SVector [sin(theta[i]) .* cos(phi[i]), sin(theta[i]) .* sin(phi[i]), cos(theta[i])] for i = 1:n]
end

units = :real
units = :reduced

const T = 120.0 # °K
const kb = 1.38e-23 # J/K
const ϵ = T * kb
const σ = 3.4e-10 # m
const ρ = 1374 # kg/m^3
const m = 39.95 * 1.6747 * 1e-27 # kg
const N = 125#floor(Int, ρ * L^3 / m)
const L = (m*N/ρ)^(1/3)#10.229σ
const R = 2.25σ
const v_dev = sqrt(kb * T / m)
const τ = 1e-14 # σ/v
const t1 = 0τ
const t2 = 2000τ
#bodies = generate_bodies_randomly(N, m, v_dev, L)
bodies = generate_bodies_in_cell_nodes(N, m, v_dev, L)
#bodies = generate_bodies_in_line(N, m, v_dev, L)
jl_parameters = LennardJonesParameters(ϵ, σ, R)
pbc = CubicPeriodicBoundaryConditions(L)
thermostat = AndersenThermostat(0.02, T, kb)
lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => jl_parameters));
simulation = NBodySimulation(lj_system, (t1, t2), pbc, thermostat);
#result = run_simulation(simulation, Tsit5())
result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)

#=
using Plots
import GR
(rs, grf) = rdf(result)
(ts, dr2) = msd(result)
plot(rs/σ, grf, xlim=[0, 0.4999L/σ], label=["Radial distribution function"],ylabel="g(r)", xlabel="r/sigma")

using JLD
time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))
#save("D:/water $Nactual molecules $timesteps steps.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2)
save("D:/liquid argon $Nactual molecules $timesteps steps $time_now.jld", "rs", rs, "grf", grf, "ts", ts, "dr2", dr2, "e_tot", e_tot, "e_kin", e_kin, "e_pot", e_pot)
=#
#=
using Plots
import GR
time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))
@time animate(result, "D:/$Nactual liquid argon particles with $timesteps timesteps $time_now.gif")

#plot(simulation)
=#


text:
include("make_structure.jl")

Documenter.doctest(Learn)

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
if Pair != Base.Pair
import Base: Base, IOContext, string, join, sprint
IOContext(io::IO, KV::Pair) = IOContext(io, Base.Pair(KV[1], KV[2]))
length(s::String) = Base.length(s)
^(s::String, i::Int) = Base.:^(s, i)
end

function Base.show(io::IO, cfg::CFG)
    foreach(pairs(cfg.blocks)) do (idx, block)
        Base.println("$idx\t=>\t", join(block.succs, ", "))
    end
end

print_ssa(io::IO, val) = isa(val, SSAValue) ? Base.print(io, "%$(val.id)") :
                         isa(val, Argument) ? Base.print(io, "%%$(val.n)") : Base.print(io, val)
function print_node(io::IO, idx, stmt, used, maxsize; color = true, print_typ=true)
    if idx in used
        pad = " "^(maxsize-length(string(idx)))
        Base.print(io, "%$idx $pad= ")
    else
        Base.print(io, " "^(maxsize+4))
    end
    if isa(stmt, PhiNode)
        args = map(1:length(stmt.edges)) do i
            e = stmt.edges[i]
            v = !isassigned(stmt.values, i) ? "#undef" :
                sprint() do io′
                    print_ssa(io′, stmt.values[i])
                end
            "$e => $v"
        end
        Base.print(io, "φ ", '(', join(args, ", "), ')')
    elseif isa(stmt, PhiCNode)
        Base.print(io, "φᶜ ", '(', join(map(x->sprint(print_ssa, x), stmt.values), ", "), ')')
    elseif isa(stmt, PiNode)
        Base.print(io, "π (")
        print_ssa(io, stmt.val)
        Base.print(io, ", ")
        if color
            Base.printstyled(io, stmt.typ, color=:red)
        else
            Base.print(io, stmt.typ)
        end
        Base.print(io, ")")
    elseif isa(stmt, UpsilonNode)
        Base.print(io, "ϒ (")
        isdefined(stmt, :val) ?
            print_ssa(io, stmt.val) :
            Base.print(io, "#undef")
        Base.print(io, ")")
    elseif isa(stmt, ReturnNode)
        if !isdefined(stmt, :val)
            Base.print(io, "unreachable")
        else
            Base.print(io, "return ")
            print_ssa(io, stmt.val)
        end
    elseif isa(stmt, GotoIfNot)
        Base.print(io, "goto ", stmt.dest, " if not ")
        print_ssa(io, stmt.cond)
    elseif isexpr(stmt, :call)
        print_ssa(io, stmt.args[1])
        Base.print(io, "(")
        Base.print(io, join(map(arg->sprint(io->print_ssa(io, arg)), stmt.args[2:end]), ", "))
        Base.print(io, ")")
        if print_typ && stmt.typ !== Any
            Base.print(io, "::$(stmt.typ)")
        end
    elseif isexpr(stmt, :new)
        Base.print(io, "new(")
        Base.print(io, join(map(arg->sprint(io->print_ssa(io, arg)), stmt.args), ", "))
        Base.print(io, ")")
    else
        Base.print(io, stmt)
    end
end

function Base.show(io::IO, code::IRCode)
    io = IOContext(io, :color=>true)
    used = IdSet{Int}()
    Base.println(io, "Code")
    foreach(stmt->scan_ssa_use!(used, stmt), code.stmts)
    cfg = code.cfg
    max_bb_idx_size = length(string(length(cfg.blocks)))
    bb_idx = 1
    if any(i->!isassigned(code.new_nodes, i), 1:length(code.new_nodes))
        printstyled(io, :red, "ERROR: New node array has unset entry\n")
    end
    new_nodes = code.new_nodes[filter(i->isassigned(code.new_nodes, i), 1:length(code.new_nodes))]
    foreach(nn -> scan_ssa_use!(used, nn.node), new_nodes)
    perm = sortperm(new_nodes, by = x->x[1])
    new_nodes_perm = Iterators.Stateful(perm)

    if isempty(used)
        maxsize = 0
    else
        maxused = maximum(used)
        maxsize = length(string(maxused))
    end

    for idx in eachindex(code.stmts)
        if !isassigned(code.stmts, idx)
            # This is invalid, but do something useful rather
            # than erroring, to make debugging easier
            printstyled(io, :red, "UNDEF\n")
            continue
        end
        stmt = code.stmts[idx]
        bbrange = cfg.blocks[bb_idx].stmts
        bbrange = bbrange.first:bbrange.last
        bb_pad = max_bb_idx_size - length(string(bb_idx))
        bb_start_str = string("$(bb_idx) ",length(cfg.blocks[bb_idx].preds) <= 1 ? "─" : "┄",  "─"^(bb_pad)," ")
        if idx != last(bbrange)
            if idx == first(bbrange)
                Base.print(io, bb_start_str)
            else
                Base.print(io, "│  "," "^max_bb_idx_size)
            end
        end
        print_sep = false
        if idx == last(bbrange)
            print_sep = true
        end
        floop = true
        while !isempty(new_nodes_perm) && new_nodes[peek(new_nodes_perm)][1] == idx
            node_idx = popfirst!(new_nodes_perm)
            new_node = new_nodes[node_idx]
            node_idx += length(code.stmts)
            if print_sep
                if floop
                    Base.print(io, bb_start_str)
                else
                    Base.print(io, "│  "," "^max_bb_idx_size)
                end
            end
            print_sep = true
            floop = false
            print_ssa_typ = !isa(new_node.node, PiNode) && node_idx in used
            Base.with_output_color(:yellow, io) do io′
                print_node(io′, node_idx, new_node.node, used, maxsize; color = false,
                    print_typ=!print_ssa_typ || (isa(new_node.node, Expr) && typ != new_node.node.typ))
            end
            if print_ssa_typ
                Base.printstyled(io, "::$(new_node.typ)", color=:red)
            end
            Base.println(io)
        end
        if print_sep
            if idx == first(bbrange) && floop
                Base.print(io, bb_start_str)
            else
                Base.print(io, idx == last(bbrange) ? string("└", "─"^(1+max_bb_idx_size), " ") :
                    string("│  ", " "^max_bb_idx_size))
            end
        end
        if idx == last(bbrange)
            bb_idx += 1
        end
        if !isassigned(code.types, idx)
            # Again, this is an error, but can happen if passes don't update their type information
            printstyled(io, "::UNDEF", color=:red)
            println(io)
            continue
        end
        typ = code.types[idx]
        print_ssa_typ = !isa(stmt, PiNode) && idx in used
        print_node(io, idx, stmt, used, maxsize,
            print_typ=!print_ssa_typ || (isa(stmt, Expr) && typ != stmt.typ))
        if print_ssa_typ
            Base.printstyled(io, "::$(typ)", color=:red)
        end
        Base.println(io)
    end
end


text:
using LinearAlgebra

@testset "lnuclear regularizer" begin
    x = [
        5.0 0.0 0.0 0.0 0.0 0.0
        0.0 4.0 0.0 0.0 0.0 0.0
        0.0 0.0 2.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0
    ]
    reg = regularizer_lnuclear(0.5)

    @test g(reg, x) == 5.5

    y, Man = prox_αg(reg, x, 4)
    @test embed(Man, y) == [
        3.0 0.0 0.0 0.0 0.0 0.0
        0.0 2.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0
        0.0 0.0 0.0 0.0 0.0 0.0
    ]
    @test Man == FixedRankMatrices{5,6,2,ℝ}()

    comparison = check_regularizer_gradient_hessian(Man, reg)
    for (k, slope) in [(:frgrad, 2), (:frhess, 3)]
        regressiondata = remove_small_functionvals(comparison[k])
        @test length(regressiondata) != 0
        slope, residual = get_sloperesidual(regressiondata)
        @test isapprox(slope, slope; atol = 1e-1)
        @test isapprox(residual, 0; atol = 5e-2)
    end
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module GnuTLS

using Compat
using Base.Meta

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

import Base: isopen, write, read, readall, readavailable, close, show, nb_available, eof

const GNUTLS_MINSECURE_VER = @compat Dict(
	v"2.12"	=>	23,
	v"3.1"	=>	28,
	v"3.2"	=>	20,
	v"3.3"	=>	10,
	v"3.4"	=>	1
)

const gnutls_version = convert(VersionNumber,bytestring(ccall((:gnutls_check_version,gnutls),Ptr{UInt8},(Ptr{UInt8},),C_NULL)))

function versionisdeprecated(v::VersionNumber)
	majmin = convert(VersionNumber,"$(v.major).$(v.minor)")
	patch = v.patch
	if haskey(GNUTLS_MINSECURE_VER, majmin)
		return GNUTLS_MINSECURE_VER[majmin] > patch
	else
		return true
	end
end

macro gnutls_since(v,f)
	if isexpr(v,:macrocall) && v.args[1] == symbol("@v_str")
		v = convert(VersionNumber,v.args[2])
	end
	if gnutls_version < v
		msg = """This function is only supported in GnuTLS versions > $v.
				 To force GnuTLS.jl to build a more recent version, you may set
				 ENV[\"GNUTLS_VERSION\"] and run Pkg.fixup()"""
		body = quote
			error($msg)
		end
		if isexpr(f,:function) || isexpr(f,:(=))
			@assert isexpr(f.args[1],:call)
			return esc(Expr(f.head,f.args[1],body))
		else
			return nothing
		end
	end
	esc(f)
end

# GnuTLS Error handling

include("errormap.jl")

immutable GnuTLSException
	msg::ASCIIString
	code::Int32
end

GnuTLSException(code::Integer) = GnuTLSException("",code)

show(io::IO,err::GnuTLSException) = print(io,"GnuTLS Exception: ",err.msg,error_codes[@compat Int(err.code)][1],"(",string(err.code),"): ",error_codes[@compat Int(err.code)][2])

gnutls_error(msg::ASCIIString,err::Int32) = err < 0 ? throw(GnuTLSException(err)) : nothing
gnutls_error(err::Int32) = err < 0 ? throw(GnuTLSException(err)) : nothing

# Chrome does not Handle EOF properly, but instead just shuts down the transport. This function
# is a heuristic to determine whether that's what we're dealing with. Using this functions
# to ignore errors in HTTPS is allowed by RFC2818.
function is_premature_eof(code::Integer)
	if gnutls_version < v"3.0"
		# GNUTLS_E_UNEXPECTED_PACKAGE_LENGTH
		return code == -9
	else
		# GNUTLS_E_PREMATURE_TERMINATION
		return code == -110
	end
end

is_premature_eof(err::GnuTLSException) = is_premature_eof(err.code)

# GnuTLS Session support

const GNUTLS_SERVER = 1
const GNUTLS_CLIENT = 1<<1
const GNUTLS_DATAGRAM = 1<<2
const GNUTLS_NONBLOCK = 1<<3
const GNUTLS_NO_EXTENSIONS = 1<<4

type Session <: IO
	handle::Ptr{Void}
	open::Bool
	# For rooting purposes
	read::IO
	write::IO
	function Session(isserver::Bool = false)
		x = Array(Ptr{Void},1)
		gnutls_error(ccall((:gnutls_init,gnutls),Int32,(Ptr{Ptr{Void}}, UInt32),x,isserver?GNUTLS_SERVER:GNUTLS_CLIENT))
		ret = new(x[1],false)
		finalizer(ret,free_session)
		ret
	end
end

const GNUTLS_SHUT_RDWR = 0
const GNUTLS_SHUT_WR = 1

free_session(s::Session) = ccall((:gnutls_deinit,gnutls),Void,(Ptr{Void},),s.handle)
isopen(s::Session) = (s.open && (isopen(s.read) || isopen(s.write)))
function close(s::Session)
	ret::Int32 = 0
	try # The remote might very well simply shut the stream rather than acknowledge the closure
		ret = ccall((:gnutls_bye,gnutls), Int32, (Ptr{Void},Int32), s.handle, GNUTLS_SHUT_WR)
	end
	gnutls_error(ret)
	s.open=false
	close(s.read)
	close(s.write)
	nothing
end

const GNUTLS_PK_UNKNOWN = 0
const GNUTLS_PK_RSA = 1
const GNUTLS_PK_DSA = 2
const GNUTLS_PK_DH = 3
const GNUTLS_PK_EC = 4

const GNUTLS_SEC_PARAM_INSECURE = -20
const GNUTLS_SEC_PARAM_EXPORT = -15
const GNUTLS_SEC_PARAM_VERY_WEAK = -12
const GNUTLS_SEC_PARAM_WEAK = -10
const GNUTLS_SEC_PARAM_UNKNOWN = 0
const GNUTLS_SEC_PARAM_LOW = 1
const GNUTLS_SEC_PARAM_LEGACY = 2
const GNUTLS_SEC_PARAM_NORMAL = 3
const GNUTLS_SEC_PARAM_HIGH = 4
const GNUTLS_SEC_PARAM_ULTRA = 5

type DHParameters
	handle::Ptr{Void}
	function DHParameters(handle)
		ret = new(handle)
		finalizer(ret,free_dh_parameters)
		ret
	end
end

free_dh_parameters(dh::DHParameters) = ccall((:gnutls_dh_params_deinit,gnutls),Void,(Ptr{Void},),dh.handle)

function generate_dh_parameters(;sec_level=GNUTLS_SEC_PARAM_NORMAL)
	x = Array(Ptr{Void},1)
	ccall((:gnutls_dh_params_generate2,gnutls),Int32,(Ptr{Void},UInt32),x,
		ccall((:gnutls_sec_param_to_pk_bits,gnutls), UInt32, (Int32, UInt32),GNUTLS_PK_DH,GNUTLS_SEC_PARAM_NORMAL))
	DHParameters(x[1])
end

type CertificateStore
	handle::Ptr{Void}
	# for rooting purposes
	dh_parameters::DHParameters
	function CertificateStore()
		x = Array(Ptr{Void},1)
		gnutls_error(ccall((:gnutls_certificate_allocate_credentials,gnutls),Int32,(Ptr{Ptr{Void}},),x))
		ret = new(x[1])
		finalizer(ret,free_certificate_store)
		ret
	end
end

function set_dh_parameters(c::CertificateStore,dh::DHParameters)
	c.dh_parameters = dh
	ccall((:gnutls_certificate_set_dh_params,gnutls),Void,(Ptr{Void},Ptr{Void}),c.handle,dh.handle)
end

free_certificate_store(x::CertificateStore) = ccall((:gnutls_certificate_free_credentials,gnutls),Void,(Ptr{Void},),x.handle)
@gnutls_since v"3.0" function set_system_trust!(c::CertificateStore)
	ret = ccall((:gnutls_certificate_set_x509_system_trust,gnutls),Int32,(Ptr{Void},),c.handle)
	if ret == -1250
		return false
	end
	gnutls_error(ret)
	true
end

function add_trusted_ca(c::CertificateStore,file,isPEM=false)
	gnutls_error(ccall((:gnutls_certificate_set_x509_trust_file,gnutls),Int32,(Ptr{Void},Ptr{UInt8},Int32),c.handle,file,isPEM?1:0))
end

function load_certificate(c::CertificateStore,certfile::AbstractString,keyfile::AbstractString,isPEM=false)
	gnutls_error(ccall((:gnutls_certificate_set_x509_key_file,gnutls),Int32,(Ptr{Void},Ptr{UInt8},Ptr{UInt8},Int32),c.handle,certfile,keyfile,isPEM?1:0))
end

@gnutls_since v"3.0" begin
	const system_certificate_store = CertificateStore()
end

@gnutls_since v"3.0" system_trust() = system_certificate_store()

type Certificate
	handle::Ptr{Void}
	function Certificate()
		x = Array(Ptr{Void},1)
		gnutls_error(ccall((:gnutls_x509_crt_init,gnutls),Int32,(Ptr{Void},),x))
		ret = new(x[1])
		finalizer(ret,free_certificate)
		ret
	end
end
free_certificate(x::Certificate) = ccall((:gnutls_x509_crt_deinit,gnutls),Void,(Ptr{Void},),x.handle)

function import_certificate(p)
	c = Certificate()
	gnutls_error(ccall((:gnutls_x509_crt_import,gnutls),Int32,(Ptr{Void},Ptr{Void},Int32),c.handle,p,0))
	c
end

immutable DistinguishedName
	handle::Ptr{Void}
	cert::Certificate
end

immutable Datum
	data::Ptr{UInt8}
	size::UInt32
end

immutable AvaSt
	oid::Datum
	value::Datum
	value_tag::Culong
end

immutable Ava
	ava::AvaSt
	dn::DistinguishedName
end

function subject(c::Certificate)
	x = Array(Ptr{Void},1)
	gnutls_error(ccall((:gnutls_x509_crt_get_subject,gnutls),Int32,(Ptr{Void},Ptr{Ptr{Void}}),c.handle,x))
	DistinguishedName(x[1],c)
end

import Base: getindex

#(gnutls_x509_dn_t dn, Int irdn, Int iava, gnutls_x509_ava_st * ava)
function getindex(dn::DistinguishedName,irdn::Integer,iava::Integer)
	x = Array(AvaSt,1)
	gnutls_error(ccall((:gnutls_x509_dn_get_rdn_ava,gnutls),Int32,(Ptr{Void},Int32,Int32,Ptr{AvaSt}),dn.handle,irdn,iava,x))
	Ava(x[1],dn)
end

function subject_alt_name(c::Certificate,seq = 0)
	s = Array(Csize_t,1)
	err = ccall((:gnutls_x509_crt_get_subject_alt_name,gnutls),Int32,(Ptr{Void},Cuint,Ptr{UInt8},Ptr{Csize_t},Ptr{Cuint}),c.handle,seq,C_NULL,s,C_NULL)
	@assert error_codes[@compat Int(err)][1] == "GNUTLS_E_SHORT_MEMORY_BUFFER"
	a = Array(UInt8,s[1])
	gnutls_error(ccall((:gnutls_x509_crt_get_subject_alt_name,gnutls),Int32,(Ptr{Void},Cuint,Ptr{UInt8},Ptr{Csize_t},Ptr{Cuint}),c.handle,seq,a,s,C_NULL))
	bytestring(a[1:(end-1)])
end

gnutls_free(d::Datum) = ccall(:jl_gc_counted_free,Void,(Ptr{UInt8},),d.data)

const GNUTLS_CRT_PRINT_FULL 			= 0
const GNUTLS_CRT_PRINT_ONELINE 			= 1
const GNUTLS_CRT_PRINT_UNSIGNED_FULL	= 2
const GNUTLS_CRT_PRINT_COMPACT 			= 3
const GNUTLS_CRT_PRINT_FULL_NUMBERS 	= 4

function show(io::IO,c::Certificate)
	a = Array(Datum,1)
	gnutls_error(ccall((:gnutls_x509_crt_print,gnutls),Int32,(Ptr{Void},Int32,Ptr{Datum}),c.handle,GNUTLS_CRT_PRINT_FULL,a))
	print(io,bytestring(a[1].data,a[1].size))
	gnutls_free(a[1])
end

import Base: bytestring

function bytestring(d::Datum)
	bytestring(d.data,d.size)
end

function poll_readable{S<:IO}(io::S,ms::UInt32)
	c = Condition()
	@async begin
		sleep(ms/1000)
		notify(c,@compat Int32(0))
	end
	@async begin
		err = wait(readnotify(io))
		notify(c,err.code == 0 ? 1 : 0)
	end
	wait(c)::Int32
end

function read_ptr(strm_ref,ptr::Ptr{UInt8},size::Csize_t)
        io = unsafe_pointer_to_objref(strm_ref)
	Base.wait_readnb(io,1)
	n = min(nb_available(io.buffer),size)
	@compat read!(io,pointer_to_array(ptr,Int(n)))
	ret = convert(Cssize_t,n)
	ret::Cssize_t
end

function write_ptr(strm_ref,ptr::Ptr{UInt8},size::Csize_t)
        io = unsafe_pointer_to_objref(strm_ref)
        ret = convert(Cssize_t, Base.write(io, ptr, size))
        ret::Cssize_t
end

function associate_stream{S<:IO,T<:IO}(s::Session, read_strm::S, write_strm::T=read_strm)
	s.read = read_strm
	s.write = write_strm
	if write_strm == read_strm
		ccall((:gnutls_transport_set_ptr,gnutls),Void,(Ptr{Void},Any),s.handle,read_strm)
	else
		ccall((:gnutls_transport_set_ptr2,gnutls),Void,(Ptr{Void},Any,Any),s.handle,read_strm,write_strm)
	end

	@gnutls_since v"3.0" ccall((:gnutls_transport_set_pull_timeout_function,gnutls),Void,(Ptr{Void},Ptr{Void}),s.handle,cfunction(poll_readable,Int32,(S, UInt32)))
	ccall((:gnutls_transport_set_pull_function,gnutls),Void,(Ptr{Void},Ptr{Void}),s.handle,cfunction(read_ptr,Cssize_t,(Ptr{Void},Ptr{UInt8},Csize_t)))
	ccall((:gnutls_transport_set_push_function,gnutls),Void,(Ptr{Void},Ptr{Void}),s.handle,cfunction(write_ptr,Cssize_t,(Ptr{Void},Ptr{UInt8},Csize_t)))
end

handshake!(s::Session) = (gnutls_error(ccall((:gnutls_handshake,gnutls),Int32,(Ptr{Void},),s.handle));s.open = true; nothing)
function set_priority_string!(s::Session,priority::ASCIIString="NORMAL")
	x = Array(Ptr{UInt8},1)
	old_ptr = pointer(priority)
	ret = ccall((:gnutls_priority_set_direct,gnutls),Int32,(Ptr{Void},Ptr{UInt8},Ptr{Ptr{UInt8}}),s.handle,priority,x)
	offset = x[1] - old_ptr
	gnutls_error("At priority string offset $offset: ",ret)
end

function get_peer_certificate(s::Session)
	a = Array(UInt32,1)
	p = ccall((:gnutls_certificate_get_peers,gnutls),Ptr{Void},(Ptr{Void},Ptr{UInt32}),s.handle,a)
	p != C_NULL ? import_certificate(p) : nothing
end

const GNUTLS_CRD_CERTIFICATE = 1
const GNUTLS_CRD_ANON 		 = 2
const GNUTLS_CRD_SRP		 = 3
const GNUTLS_CRD_PSK 		 = 4
const GNUTLS_CRD_IA 		 = 5

set_credentials!(s::Session,c::CertificateStore) = gnutls_error(ccall((:gnutls_credentials_set,gnutls),Int32,(Ptr{Void},Int32,Ptr{Void}),s.handle,GNUTLS_CRD_CERTIFICATE,c.handle))

const GNUTLS_CERT_IGNORE = 0
const GNUTLS_CERT_REQUEST = 1
const GNUTLS_CERT_REQUIRE = 2

function set_prompt_client_certificate!(s::Session,required::Bool = true)
	ccall((:gnutls_certificate_server_set_request,gnutls),Void,(Ptr{Void},Int32),s.handle,required?GNUTLS_CERT_REQUIRE:GNUTLS_CERT_REQUEST)
end

function write(io::Session, data::Ptr{UInt8}, size::Integer)
	total = 0
	while total < length(data)
		ret = ccall((:gnutls_record_send,gnutls), Int, (Ptr{Void},Ptr{UInt8},Csize_t), io.handle, data+total, size-total)
		if ret < 0
			gnutls_error(@compat Int32(ret))
		end
		total += ret
	end
	total
end

function write(io::Session, data::Array{UInt8,1})
	total = 0
	while total < length(data)
		ret = ccall((:gnutls_record_send,gnutls), Int, (Ptr{Void},Ptr{UInt8},Csize_t), io.handle, pointer(data,total+1), length(data)-total)
		if ret < 0
			gnutls_error(@compat Int32(ret))
		end
		total += ret
	end
	total
end

function read(io::Session, data::Array{UInt8,1})
	total = 0
	while total < length(data)
		ret = ccall((:gnutls_record_recv,gnutls), Int, (Ptr{Void},Ptr{UInt8},Csize_t), io.handle, data, length(data))
		if ret < 0
			gnutls_error(@compat Int32(ret))
		elseif ret == 0
			throw(EOFError())
		end
		total += ret
	end
	data
end

checkpending(io::Session)  = convert(Bool, ccall((:gnutls_record_check_pending, gnutls), Int, (Ptr{Void},), io.handle))
recordmaxsize(io::Session)  = ccall((:gnutls_record_get_max_size, gnutls), Int, (Ptr{Void},), io.handle)
recordcheckcorked(io::Session)  = ccall((:gnutls_record_check_corked, gnutls), Int, (Ptr{Void},), io.handle)

function readtobuf(io::Session,buf::IOBuffer,nb)
	Base.ensureroom(buf,@compat Int(nb))
	ret = ccall((:gnutls_record_recv,gnutls), Int, (Ptr{Void},Ptr{UInt8},Csize_t), io.handle, pointer(buf.data,buf.size+1), nb)
	iseof = is_premature_eof(ret)
	if ret < 0 && !(iseof)
		gnutls_error(@compat Int32(ret))
	elseif ret == 0  || iseof
		close(io)
		return false
	else
		buf.size += ret
		return true
	end
end

nb_available(s::Session) = ccall((:gnutls_record_check_pending,gnutls), Csize_t, (Ptr{Void},), s.handle)

eof(s::Session) = (nb_available(s) == 0 && eof(s.read))

function readavailable(io::Session)
	buf = IOBuffer(Array(UInt8,0),true,true,true,true,typemax(Int))
	n = nb_available(io)
	if n == 0
		readtobuf(io,buf,1)
		n = nb_available(io)
	end
	while n>0
		if !readtobuf(io,buf,n)
			break
		end
		n = nb_available(io)
	end
	takebuf_array(buf)
end

function readall(io::Session)
	arr = UInt8[]
	while isopen(io)
		append!(arr, readavailable(io))
	end
	arr
end

function write{T}(s::Session, a::Array{T})
    if isbits(T)
        write(s,reinterpret(UInt8,a))
    else
        invoke(write, (IO, Array), s, a)
    end
end

# GnuTLS initialization

function logging_func(level::Int32,msg::Ptr{UInt8})
	println(bytestring(msg))
	nothing
end

function init()
	if versionisdeprecated(gnutls_version)
		msg = """This version of the GnuTLS library ($gnutls_version) is deprecated
			and contains known security vulnerabilities. Please upgrade to a
			more recent version."""
		warn(msg)
	end
	ccall((:gnutls_global_set_mem_functions,gnutls),Void,
		(Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),
		  cglobal(:jl_gc_counted_malloc),   # Malloc
		  cglobal(:jl_gc_counted_malloc),   # Secure Malloc
		  C_NULL,							# is_secure (may be NULL)
          cglobal(:realloc),  				# Realloc
          cglobal(:jl_gc_counted_free))		# Free
	ccall((:gnutls_global_set_log_function,gnutls),Void,
		(Ptr{Void},),cfunction(logging_func,Void,(Int32,Ptr{UInt8})))
	ccall((:gnutls_global_init,gnutls),Int32,())
    @gnutls_since v"3.0" begin
        global const has_system_trust = set_system_trust!(system_certificate_store)
    end
end
deinit() = ccall((:gnutls_global_deinit,gnutls),Void,())

read(io::Session, ::Type{UInt8}) = (x=Array(UInt8,1);read(io,x);x[1])

export  handshake!, associate_stream, set_priority_string!, set_credentials!, CertificateStore, Certificate

end

GnuTLS.init()


text:
# Test spatial weights matrix
@testset "Spatial Weights" begin

    # Create 3x3 polygons
    polygons = Polygon[]
    for i in 1:3
        for j in 1:3
            push!(polygons, Polygon([[[0.0 + i, 0.0 + j], [1.0 + i, 0.0 + j], [1.0 + i, 1.0 + j], [0.0 + i, 1.0 + j]]]) )
        end
    end

    # Rook contiguity
    @testset "Rook Contiguity" begin
        W = polyneigh(polygons, criterion = :Rook)
        @test isa(W, SpatialWeights)

        @test nobs(W) == 9
        @test mean(W) ≈ 2.6666666666666665 atol = 1e-5
        @test minimum(W) == 2
        @test maximum(W) == 4
        @test median(W) == 3.0
        @test nislands(W) == 0
        @test islands(W) == []

        @test cardinalities(W) == [2, 3, 2, 3, 4, 3, 2, 3, 2]

        @test neighbors(W, 1) == [2, 4]
        @test neighbors(W, 2) == [1, 3, 5]
        @test neighbors(W, 3) == [2, 6]
        @test neighbors(W, 4) == [1, 5, 7]
        @test neighbors(W, 5) == [2, 4, 6, 8]
        @test neighbors(W, 6) == [3, 5, 9]
        @test neighbors(W, 7) == [4, 8]
        @test neighbors(W, 8) == [5, 7, 9]
        @test neighbors(W, 9) == [6, 8]

        # Default is row-standardized weights
        @test wtransformation(W) == :row

        @test weights(W, 1) == [0.5, 0.5]
        @test weights(W, 2) == [1/3, 1/3, 1/3]
        @test weights(W, 3) == [0.5, 0.5]
        @test weights(W, 4) == [1/3, 1/3, 1/3]
        @test weights(W, 5) == [1/4, 1/4, 1/4, 1/4]
        @test weights(W, 6) == [1/3, 1/3, 1/3]
        @test weights(W, 7) == [0.5, 0.5]
        @test weights(W, 8) == [1/3, 1/3, 1/3]
        @test weights(W, 9) == [0.5, 0.5]

        # Transform to Binary weights
        wtransform!(W, :binary)

        @test weights(W, 1) == [1, 1]
        @test weights(W, 2) == [1, 1, 1]
        @test weights(W, 3) == [1, 1]
        @test weights(W, 4) == [1, 1, 1]
        @test weights(W, 5) == [1, 1, 1, 1]
        @test weights(W, 6) == [1, 1, 1]
        @test weights(W, 7) == [1, 1]
        @test weights(W, 8) == [1, 1, 1]
        @test weights(W, 9) == [1, 1]

        # Copy transformed to Row-standardization
        Wrow = wtransform(W, :row)

        @test wtransformation(Wrow) == :row

        @test weights(Wrow, 1) == [0.5, 0.5]
        @test weights(Wrow, 2) == [1/3, 1/3, 1/3]
        @test weights(Wrow, 3) == [0.5, 0.5]
        @test weights(Wrow, 4) == [1/3, 1/3, 1/3]
        @test weights(Wrow, 5) == [1/4, 1/4, 1/4, 1/4]
        @test weights(Wrow, 6) == [1/3, 1/3, 1/3]
        @test weights(Wrow, 7) == [0.5, 0.5]
        @test weights(Wrow, 8) == [1/3, 1/3, 1/3]
        @test weights(Wrow, 9) == [0.5, 0.5]

        # In-place row standardization
        wtransform!(W, :row)

        @test wtransformation(W) == :row

        @test weights(W, 1) == [0.5, 0.5]
        @test weights(W, 2) == [1/3, 1/3, 1/3]
        @test weights(W, 3) == [0.5, 0.5]
        @test weights(W, 4) == [1/3, 1/3, 1/3]
        @test weights(W, 5) == [1/4, 1/4, 1/4, 1/4]
        @test weights(W, 6) == [1/3, 1/3, 1/3]
        @test weights(W, 7) == [0.5, 0.5]
        @test weights(W, 8) == [1/3, 1/3, 1/3]
        @test weights(W, 9) == [0.5, 0.5]
    end

    # Queen contiguity
    @testset "Queen Contiguity" begin
        W = polyneigh(polygons, criterion = :Queen)
        @test isa(W, SpatialWeights)

        @test nobs(W) == 9
        @test mean(W) ≈ 4.444444444444445 atol = 1e-5
        @test minimum(W) == 3
        @test maximum(W) == 8
        @test median(W) == 5.0
        @test nislands(W) == 0
        @test islands(W) == []

        @test cardinalities(W) == [3, 5, 3, 5, 8, 5, 3, 5, 3]

        @test neighbors(W, 1) == [2, 4, 5]
        @test neighbors(W, 2) == [1, 3, 4, 5, 6]
        @test neighbors(W, 3) == [2, 5, 6]
        @test neighbors(W, 4) == [1, 2, 5, 7, 8]
        @test neighbors(W, 5) == [1, 2, 3, 4, 6, 7, 8, 9]
        @test neighbors(W, 6) == [2, 3, 5, 8, 9]
        @test neighbors(W, 7) == [4, 5, 8]
        @test neighbors(W, 8) == [4, 5, 6, 7, 9]
        @test neighbors(W, 9) == [5, 6, 8]

        # Row-standardized weights
        @test wtransformation(W) == :row

        @test weights(W, 1) == [1/3, 1/3, 1/3]
        @test weights(W, 2) == [1/5, 1/5, 1/5, 1/5, 1/5]
        @test weights(W, 3) == [1/3, 1/3, 1/3]
        @test weights(W, 4) == [1/5, 1/5, 1/5, 1/5, 1/5]
        @test weights(W, 5) == [1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8]
        @test weights(W, 6) == [1/5, 1/5, 1/5, 1/5, 1/5]
        @test weights(W, 7) == [1/3, 1/3, 1/3]
        @test weights(W, 8) == [1/5, 1/5, 1/5, 1/5, 1/5]
        @test weights(W, 9) == [1/3, 1/3, 1/3]

        # Binary weights
        wtransform!(W, :binary)
        @test wtransformation(W) == :binary

        @test weights(W, 1) == [1, 1, 1]
        @test weights(W, 2) == [1, 1, 1, 1, 1]
        @test weights(W, 3) == [1, 1, 1]
        @test weights(W, 4) == [1, 1, 1, 1, 1]
        @test weights(W, 5) == [1, 1, 1, 1, 1, 1, 1, 1]
        @test weights(W, 6) == [1, 1, 1, 1, 1]
        @test weights(W, 7) == [1, 1, 1]
        @test weights(W, 8) == [1, 1, 1, 1, 1]
        @test weights(W, 9) == [1, 1, 1]

        @test_nowarn show(IOBuffer(), W)
    end

    # Test Multipolygon with multipolygon 1 at the left and at the right
    @testset "Multipolygon" begin
        mpolygons = MultiPolygon[]
        push!(mpolygons,
            MultiPolygon([[[[1.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0]]],
                        [[[5.0, 0.0], [6.0, 0.0], [6.0, 1.0], [5.0, 1.0]]]])
            )

        for i in 2:4
            push!(mpolygons, MultiPolygon(Polygon([[[0.0 + i, 0.0], [1.0 + i, 0.0 ], [1.0 + i, 1.0], [0.0 + i, 1.0]]])) )
        end

        W = polyneigh(mpolygons)

        @test neighbors(W, 1) == [2, 4]
        @test neighbors(W, 2) == [1, 3]
        @test neighbors(W, 3) == [2, 4]
        @test neighbors(W, 4) == [1, 3]
    end

    # Test tolerance for edge detection
    @testset "Polygon tolerance" begin
        polygonstol = Polygon[]
        push!(polygonstol, Polygon([[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]]) )
        push!(polygonstol, Polygon([[[1.0 + 0.25, 0.0], [2.0 + 0.25, 0.0], [2.0 + 0.25, 1.0], [1.0 + 0.25, 1.0]]]) )
        push!(polygonstol, Polygon([[[2.0 + 0.50, 0.0], [3.0 + 0.50, 0.0], [3.0 + 0.50, 1.0], [2.0 + 0.50, 1.0]]]) )

        W = polyneigh(polygonstol)

        @test neighbors(W, 1) == []
        @test neighbors(W, 2) == []
        @test neighbors(W, 3) == []

        @test nislands(W) == 3
        @test islands(W) == [1, 2, 3]

        W = polyneigh(polygonstol, tol = 0.25)

        @test neighbors(W, 1) == [2]
        @test neighbors(W, 2) == [1, 3]
        @test neighbors(W, 3) == [2]

        @test nislands(W) == 0
        @test islands(W) == []
    end

    # Generate points for test
    points = Point[]
    for i in 1:3
        for j in 1:3
            push!(points, Point([0.0 + i, 0.0 + j]))
        end
    end

    # Points distance threshold
    @testset "Distance Threshold" begin
        W = dnearneigh(points, threshold = 1.5)
        @test isa(W, SpatialWeights)

        @test nobs(W) == 9
        @test mean(W) ≈ 4.444444444444445 atol = 1e-5
        @test minimum(W) == 3
        @test maximum(W) == 8
        @test median(W) == 5.0
        @test nislands(W) == 0
        @test islands(W) == []

        @test neighbors(W, 1) == [2, 4, 5]
        @test neighbors(W, 2) == [1, 3, 4, 5, 6]
        @test neighbors(W, 3) == [2, 5, 6]
        @test neighbors(W, 4) == [1, 2, 5, 7, 8]
        @test neighbors(W, 5) == [1, 2, 3, 4, 6, 7, 8, 9]
        @test neighbors(W, 6) == [2, 3, 5, 8, 9]
        @test neighbors(W, 7) == [4, 5, 8]
        @test neighbors(W, 8) == [4, 5, 6, 7, 9]
        @test neighbors(W, 9) == [5, 6, 8]
    end

    @testset "Polygon Missing" begin
        PM = Vector{Union{Missing,Polygon}}(undef, 10)
        PM[1:9] = polygons[1:9]
        PM[10] = missing

        @test typeof(PM) == Vector{Union{Missing,Polygon}}
        @test_throws ArgumentError polyneigh(PM) # Missin geometry error
    end

    # Points k near neighbors
    @testset "K near neighbors " begin
        W = knearneigh(points, k = 3)
        @test isa(W, SpatialWeights)

        @test nobs(W) == 9
        @test mean(W) ≈ 3 atol = 1e-5
        @test minimum(W) == 3
        @test maximum(W) == 3
        @test median(W) == 3.0
        @test nislands(W) == 0
        @test islands(W) == []

        @test neighbors(W, 1) == [2, 4, 5]
        @test neighbors(W, 2) == [1, 3, 5]
        @test neighbors(W, 3) == [2, 5, 6]
        @test neighbors(W, 4) == [1, 5, 7]
        @test neighbors(W, 5) == [2, 4, 8]
        @test neighbors(W, 6) == [3, 5, 9]
        @test neighbors(W, 7) == [4, 5, 8]
        @test neighbors(W, 8) == [5, 7, 9]
        @test neighbors(W, 9) == [5, 6, 8]
    end

    # Polygon mean center
    @testset "Polygon meancenter" begin
        x, y = meancenter(reggeomlattice(3, 3, direction = :rightup))
        @test x == [1.5, 2.5, 3.5, 1.5, 2.5, 3.5, 1.5, 2.5, 3.5]
        @test y == [1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 3.5, 3.5, 3.5]
    end

    # Polygon centroids
    @testset "Polygon centroid" begin
        x, y = centroid(reggeomlattice(3, 3, direction = :rightup))
        @test x == [1.5, 2.5, 3.5, 1.5, 2.5, 3.5, 1.5, 2.5, 3.5]
        @test y == [1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 3.5, 3.5, 3.5]
    end

    # Spatial Weights from Matrix
    @testset "From Matrix " begin
        W = SpatialWeights([0 1 0; 1 0 1; 0 1 0])
        @test isa(W, SpatialWeights)

        @test nobs(W) == 3
        @test mean(W) ≈ 4/3 atol = 1e-5
        @test minimum(W) == 1
        @test maximum(W) == 2
        @test median(W) == 1.0
        @test nislands(W) == 0
        @test islands(W) == []

        @test neighbors(W, 1) == [2]
        @test neighbors(W, 2) == [1, 3]
        @test neighbors(W, 3) == [2]

        # Row-standardize by default
        @test wtransformation(W) == :row
        @test weights(W, 1) == [1]
        @test weights(W, 2) == [0.5, 0.5]
        @test weights(W, 3) == [1]

        # Do not standardize
        W = SpatialWeights([0 2 0; 3 0 4; 0 5 0], standardize = false)
        @test wtransformation(W) == :original
        @test weights(W, 1) == [2]
        @test weights(W, 2) == [3, 4]
        @test weights(W, 3) == [5]

    end

    # Test conversions
    @testset "Matrix converion" begin
        W = polyneigh(polygons, criterion = :Rook)

        @test Matrix(W) ==
        [
            0.0   0.5   0.0   0.5   0.0   0.0   0.0   0.0   0.0
            1/3   0.0   1/3   0.0   1/3   0.0   0.0   0.0   0.0
            0.0   0.5   0.0   0.0   0.0   0.5   0.0   0.0   0.0
            1/3   0.0   0.0   0.0   1/3   0.0   1/3   0.0   0.0
            0.0   0.25  0.0   0.25  0.0   0.25  0.0   0.25  0.0
            0.0   0.0   1/3   0.0   1/3   0.0   0.0   0.0   1/3
            0.0   0.0   0.0   0.5   0.0   0.0   0.0   0.5   0.0
            0.0   0.0   0.0   0.0   1/3   0.0   1/3   0.0   1/3
            0.0   0.0   0.0   0.0   0.0   0.5   0.0   0.5   0.0
        ]

        # Binary weights
        wtransform!(W, :binary)
        @test Matrix(W) ==
        [
            0.0  1.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0;
            1.0  0.0  1.0  0.0  1.0  0.0  0.0  0.0  0.0;
            0.0  1.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0;
            1.0  0.0  0.0  0.0  1.0  0.0  1.0  0.0  0.0;
            0.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  0.0;
            0.0  0.0  1.0  0.0  1.0  0.0  0.0  0.0  1.0;
            0.0  0.0  0.0  1.0  0.0  0.0  0.0  1.0  0.0;
            0.0  0.0  0.0  0.0  1.0  0.0  1.0  0.0  1.0;
            0.0  0.0  0.0  0.0  0.0  1.0  0.0  1.0  0.0
        ]

        @test Matrix(sparse(W)) == Matrix(W)

    end

    # Regular lattice
    @testset "Regular lattice" begin
        regrightup = reggeomlattice(3, 3, direction = :rightup)
        regrightdown = reggeomlattice(3, 3, direction = :rightdown)
        regleftdown = reggeomlattice(3, 3, direction = :leftdown)
        regleftup = reggeomlattice(3, 3, direction = :leftup)

        @test coordinates(regrightup[1])   == [[[[1.0, 1.0], [2.0, 1.0], [2.0, 2.0], [1.0, 2.0]]]]
        @test coordinates(regrightdown[1]) == [[[[1.0, 3.0], [2.0, 3.0], [2.0, 4.0], [1.0, 4.0]]]]
        @test coordinates(regleftdown[1])  == [[[[3.0, 3.0], [4.0, 3.0], [4.0, 4.0], [3.0, 4.0]]]]
        @test coordinates(regleftup[1])    == [[[[3.0, 1.0], [4.0, 1.0], [4.0, 2.0], [3.0, 2.0]]]]

        WRook = polyneigh(regrightup, criterion = :Rook)
        @test minimum(WRook) == 2
        @test maximum(WRook) == 4
        @test median(WRook) == 3.0

        WhtQueen = polyneigh(regrightup, criterion = :Queen)
        @test minimum(WhtQueen) == 3
        @test maximum(WhtQueen) == 8
        @test median(WhtQueen) == 5.0

        regpoints = reggeomlattice(3, 3, gtype = :point)
        Wpoint = dnearneigh(regpoints, threshold = 1.0)
        @test minimum(Wpoint) == 2
        @test maximum(Wpoint) == 4
        @test median(WRook) == 3.0
    end

end


text:
"""
    abstract type AbstractArrayItem{N, T}

Abstract type for all [`Item`]s that wrap an `N`-dimensional
array with element type `T`.
"""
abstract type AbstractArrayItem{N, T} <: Item end

"""
    ArrayItem(a)

An item that contains an array.
"""
struct ArrayItem{N, T} <: AbstractArrayItem{N, T}
    data::AbstractArray{T, N}
end

Base.show(io::IO, item::ArrayItem{N, T}) where {N, T} =
    print(io, "ArrayItem{$N, $T}() of size $(size(itemdata(item)))")

"""
    MapElem(f)

Applies `f` to every element in an [`AbstractArrayItem`].
"""
struct MapElem <: Transform
    f
end

function apply(tfm::MapElem, item::AbstractArrayItem; randstate = nothing)
    return setdata(item, map(tfm.f, itemdata(item)))
end

function apply!(
        bufitem::I,
        tfm::MapElem,
        item::I;
        randstate = nothing) where I <: AbstractArrayItem
    map!(tfm.f, itemdata(bufitem), itemdata(item))
    return bufitem
end


text:
import LinearAlgebra: diag

abstract type AbstractPersistenceReduction end
mutable struct StandardReduction <: AbstractPersistenceReduction end
mutable struct TwistReduction <: AbstractPersistenceReduction end

struct Interval
    dim::Int
    generator::AbstractChain
    b::Number
    d::Number
end
Interval(dim::Int, b::Number, d::Number) = Interval(dim, EmptyChain(), b, d)
Interval(dim::Int, p::Pair) = Interval(dim, p.first, p.second)
Interval(p::Pair) = Interval(0, p)
Base.show(io::IO, intr::Interval) = print(io, "[$(intr.b),$(intr.d))")
Base.isless(i1::Interval, i2::Interval) = i1.b < i2.b ? true : ( i1.b == i2.b ? i1.d < i2.d : false )
birth(i::Interval) = i.b - i.d
death(i::Interval) = i.b + i.d
diag(i::Interval) = let c = death(i)/2.0; Interval(i.dim, i.generator, c, c) end

pair(i::Interval) = i.b => i.d
intervals(d::Int, ps::Pair...) = [Interval(d, p) for p in ps]

# Boundary matrix reduction algorithms
lastindex(col::AbstractSet) = length(col) == 0 ? -1 : maximum(col)

"""Standart reduction"""
function reduce!(::Type{StandardReduction}, ∂::Vector{<:AbstractSet})
    lowest_one_lookup = fill(-1, length(∂))
    for col in eachindex(∂)
        lowest_one = lastindex(∂[col])
        while lowest_one != -1 && lowest_one_lookup[lowest_one] != -1
            lowest = lowest_one_lookup[lowest_one]
            symdiff!(∂[col], ∂[lowest])
            lowest_one = lastindex(∂[col])
        end
        if lowest_one != -1
            lowest_one_lookup[lowest_one] = col
        end
    end
    return ∂
end

"""Twist reduction"""
function reduce!(::Type{TwistReduction}, ∂::Vector{<:AbstractSet})
    lowest_one_lookup = fill(-1, length(∂))

    for dim in maximum(map(length, ∂)):-1:1
        for col in eachindex(∂)
            if length(∂[col]) == dim
                lowest_one = lastindex(∂[col])
                while lowest_one != -1 && lowest_one_lookup[lowest_one] != -1
                    lowest = lowest_one_lookup[lowest_one]
                    symdiff!(∂[col], ∂[lowest])
                    lowest_one = lastindex(∂[col])
                end
                if lowest_one != -1
                    lowest_one_lookup[lowest_one] = col
                end
            end
        end
    end
    return ∂
end

Base.reduce(::Type{R}, ∂::Vector{<:AbstractSet}) where {R<:AbstractPersistenceReduction} = reduce!(R, deepcopy(∂))

function generate_pairs(∂::Vector{<:AbstractSet}; reduced = false)
    ridx = reduced ? 1 : 0
    births = eltype(∂)()
    ps = Pair[]
    for i in eachindex(∂)
        if length(∂[i]) > 0
            b = last(∂[i])
            d = i
            delete!(births, b)
            delete!(births, d)
            (d > b) && push!(ps, (b-ridx) => (d-ridx) )
        else
            push!(births, i)
        end
    end
    for i in births # no lowest, create semi-infinite interval
        push!(ps, (i-ridx)=>Inf)
    end
    return ps
end

"Compute raw persistence pairs (boundary matrix is reduced in a process)"
function pairs(::Type{R}, ∂::Vector{<:AbstractSet}; reduced = false) where {R <: AbstractPersistenceReduction}
    reduce!(R, ∂) # reduce  boundary matrix
    return generate_pairs(∂, reduced=reduced), ∂  # generate pairs
end
pairs(::Type{R}, flt::Filtration; reduced = false) where {R <: AbstractPersistenceReduction} =
    pairs(R, boundary(flt, reduced = reduced), reduced = reduced)

"""Return persistent diagram (birth-death pairs) per dimension."""
function intervals(flt::Filtration, R::Vector, length0=false, absolute=true)
    # resulting intervals
    intrs = Dict{Int,Vector{Interval}}()

    # compute intervals
    births = Set{Int}()
    ord = order(flt)
    for (i, (sdim, si, fv)) in enumerate(ord)
        col = R[i]
        if length(col) == 0
            push!(births, i)
        else
            b = last(col)
            d = i
            delete!(births, b)
            delete!(births, d)
            if d > b
                sdim = absolute ? ord[b][1] : ord[d][1]
                bv = ord[b][3]
                dv = ord[d][3]
                if dv > bv || length0
                    !haskey(intrs, sdim) && setindex!(intrs, Interval[], sdim)
                    push!(intrs[sdim], Interval(sdim, bv, dv))
                end
            end
        end
    end
    for i in births
        sdim, si, fv = ord[i]
        !haskey(intrs, sdim) && setindex!(intrs, Interval[], sdim)
        push!(intrs[sdim], Interval(sdim, fv, Inf))
    end

    return intrs
end

intervals(flt::Filtration; reduction=TwistReduction, length0=false, absolute=true) =
    intervals(flt, reduce!(reduction, boundary(flt)), length0, absolute)

"Calculate persistent Betti numbers for a filtration complex of dimension `dim`"
function betti(flt::Filtration, R::Vector, p::Int)
    maxdim =  dim(complex(flt))
    @assert maxdim >= p "Cannot calculate $p-dimensional Betti number for $maxdim-complex"

    # indexes of p-simplices
    p2i = findall(d->d[1]==p, order(flt))

    # the number of zero columns that correspond to p-simplices
    z = sum(map(length, R[p2i]) .== 0)

    # the number of lowest ones in rows that correspond to p-simplices
    l = length( findall(l->l ∈ p2i, map(lastindex, R)) )

    β = z - l
    return β < 0 ? 0 : β
end
betti(flt::Filtration, p::Int; reduction=TwistReduction) =
    betti(flt, reduce!(reduction, boundary(flt)), p)

"""
Persistent homology group iterator for a filtration
"""
mutable struct PersistentHomology <: AbstractHomology
    filtration::Filtration
    reduction::DataType
    R::Vector{<:AbstractSet}
end
function persistenthomology(::Type{R}, flt::Filtration;
                            reduced::Bool=false) where R<:AbstractPersistenceReduction
    RM = reduce!(R, boundary(flt, reduced = reduced))
    return PersistentHomology(flt, R, RM)
end
persistenthomology(flt::Filtration) = persistenthomology(TwistReduction, flt)

Base.show(io::IO, h::PersistentHomology) = print(io, "PersistentHomology[$(h.filtration) with $(h.reduction)]")

#
# Interface methods
#

group(h::PersistentHomology, p::Int) = betti(h.filtration, h.R, p)
intervals(h::PersistentHomology) = intervals(h.filtration, reduction=h.reduction)

#
# Iterator methods
#

Base.length(h::PersistentHomology) = dim(h.filtration.complex)+1

"""Return homology group type: dimension & Betti numbers."""
Base.eltype(h::PersistentHomology) = Tuple{Int, Int}

function Base.iterate(h::PersistentHomology, p=0)
    p > dim(h.filtration.complex) && return nothing
    βₚ = group(h, p)
    return (p, βₚ), p+1
end


text:
module GNSSDecoder

    using DocStringExtensions, GNSSSignals, Parameters, FixedPointNumbers, StaticArrays, LinearAlgebra
    using Unitful: Hz

    BUFFER_LENGTH = 1502 #Size of Buffer - All 5 Subframes + the last 2 Bits are stored
    SUBFRAME_LENGTH = 300 #Size of Subframes

    export decode,
        GNSSDecoderState,
        GPSData,
        GPSL1Constants

    include("gpsl1_structs.jl")

    """
        stores element in a bitbuffer

        $(SIGNATURES)

        # Buffer must be of Type BitArray{1}
    """
    function push_buffer(buffer::BitArray{1}, val::Bool)
        buffer = buffer >> 1
        buffer.chunks[1] = buffer.chunks[1] + val
        return buffer
    end

    """
        Decodes the bits inserted using a bit buffer
        $(SIGNATURES)

        # Details
        checks, if the buffer has the expected size of 1508 bits (length of subframe + preamble),
        then tries to find a preamble.

    """
    function decode(dc::GNSSDecoderState, data_bits::UInt64, num_bits)

        a = num_bits - dc.nb_prev
        dc.nb_prev = num_bits

        if a !== 0
            num_bits = a
        end

        for i = num_bits - 1:-1:0

            current_bit = data_bits & (1 << i)

            dc.buffer = push_buffer(dc.buffer, current_bit > 0)

            # * Count if Buffer is new to full stored buffer
            if dc.num_bits_buffered != BUFFER_LENGTH
                dc.num_bits_buffered += 1
            end

            # * Find first preamble
            if (!dc.preamble_found) && (dc.num_bits_buffered == BUFFER_LENGTH)
                dc.preamble_found = find_preamble(dc.buffer)
            end

            # * Begin decoding after preamble is confirmed
            if dc.preamble_found
                rev_buf = reverse(dc.buffer)
                words_in_subframes  = map( ##Splits Buffer in Subframes and Words
                        subfr_it->map(
                        word_it->rev_buf[(word_it * 30) + 3 + (SUBFRAME_LENGTH * subfr_it):((word_it + 1) * 30 ) + 2 + (SUBFRAME_LENGTH * subfr_it)],
                        0:9),
                        0:4) # gets the words of the subframe (word length = 30)
                 dc.prev_29 = rev_buf[1]
                 dc.prev_30 = rev_buf[2]
                 decode_words(dc, words_in_subframes)
                dc.preamble_found = false
            end
        end # end of for-loop
    end # end of decode()

    include("subframe_decoding.jl")
    include("find_preamble.jl")
    include("decode_words.jl")
    include("parity_check.jl")
    include("bin_2_dec.jl")
end# end of module


text:
using Nanoid
using Test

@testset "Nanoid.jl" begin
    @testset "Randomness" begin
        @test nanoid() != nanoid()
    end
    @testset "Length" begin
        @test length(nanoid()) == 21
        @test length(nanoid(1)) == 1
        @test length(nanoid(255)) == 255
        @test length(nanoid(25)) == 25
    end

    @testset "Errors" begin
        @test_throws ErrorException("size must be at least 1") nanoid(-1)
    end

    @testset "Custom alphabet" begin
        output = nanoid(;alphabet = "ABCDE")

        for chars in output
            @test chars in "ABCDE"
        end
    end
end


text:
function Base.show(io::IO, nd::NanoDate; sep::CharString="")
    str = string(nd; sep)
    print(io, str)
end

function Base.show(io::IO, nd::NanoDate,
                    df::DateFormat; sep::CharString="")
    str = Dates.format(nd, df; sep)
    print(io, str)
end

Base.show(io::IO, nd::NanoDate) =
    Base.show(io::IO, m::MIME"text/plain", nd::NanoDate)

Base.show(io::IO, m::MIME"text/plain", nd::NanoDate) =
    print(io, string(nd))

Base.show(io::IO, nd::NanoDate, sep::Tuple{Vararg{Char}}) =
    Base.show(io::IO, m::MIME"text/plain", nd::NanoDate, sep)

# to show values in the REPL and when using IJulia
Base.show(io::IO, m::MIME"text/plain",
          nd::NanoDate, sep::Tuple{Vararg{Char}}) =
    print(io, string(nd, sep))


text:
include("simple_diffusion.jl")

######
# Check on fully connected, prob = 1
######

g = CompleteGraph(5)
for i in 1:5
  add_vertex!(g)
end
neighbors(g, 1)

# Basic test. Watch connected vertices
result = diffusion_simulation(g,
                              1.0,
                              4,
                              to_watch=Set(1:5),
                              initial_at_risk=Set(1:5)
                              )
result
assert(result == [5, 5, 5, 5])

# Watching unconnected vertices
result = diffusion_simulation(g,
                              1.0,
                              4,
                              to_watch=Set(6:10),
                              initial_at_risk=Set(1:5)
                              )
assert( result == [0, 0, 0, 0])

# Watch subset
result = diffusion_simulation(g,
                              1.0,
                              4,
                              to_watch=Set(1:2),
                              initial_at_risk=Set(1:5)
                              )
                              result
assert( result == [2, 2, 2, 2] )

result = diffusion_simulation(g,
                              1.0,
                              4,
                              to_watch=Set(1:5),
                              initial_at_risk=Set(10)
                              )

assert( result == [0, 0, 0, 0] )

######
# Check along path graph
######

g2 = PathGraph(5)
result = diffusion_simulation(g2,
                              1.0,
                              4,
                              to_watch=Set(1:5),
                              initial_at_risk=Set(1)
                              )

assert( result == [2, 3, 4, 5] )

result = diffusion_simulation(g2,
                              1.0,
                              4,
                              to_watch=Set(1:5),
                              initial_at_risk=Set(3)
                              )

assert( result == [3, 5, 5, 5] )

# Check normalize
result = diffusion_simulation(g,
                              1.0,
                              4,
                              to_watch=Set(1:5),
                              initial_at_risk=Set(1:5),
                              normalize_p=true
                              )
result
assert(result != [5, 5, 5, 5])


text:
using Test

using FastFourierTransform

@testset "FFT Recursive Naive" begin
    for i in 1:5
        n = 2 ^ rand(0:10)
        t = complex(rand(n))
        f_direct = dft(t)
        @test idft(f_direct) ≈ t

        f_naive = fft_recursive_naive(t)
        f_recursive = fft_recursive(t)
        f_iterative = fft_iterative(t)
        @test f_naive ≈ f_direct
        @test f_recursive ≈ f_direct
        @test f_iterative ≈ f_direct

        t_naive = ifft_recursive_naive(f_naive)
        t_recursive = ifft_recursive(f_recursive)
        t_iterative = ifft_iterative(f_iterative)
        @test t_naive ≈ t
        @test t_recursive ≈ t
        @test t_iterative ≈ t

        # Polynomial multiplication
        p1 = rand(n)
        p2 = rand(n)
        x = rand(n)
        prod = @. eval_poly($Ref(p1), x) * eval_poly($Ref(p2), x)
        p = poly_mul(p1, p2)
        prod_fft = eval_poly.(Ref(p), x)
        @test prod ≈ prod_fft
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions