text:
"""
    AbstractEmbedding{D, T}

Abstract state space embedding.
"""
abstract type AbstractEmbedding{D, T} end

#####################################################
# Simple immutable type to store embedding lags.
#####################################################
struct EmbeddingLags
    lags::Vector{Int}
end

Base.getindex(rl::EmbeddingLags,
            i::Union{Int, Colon, UnitRange, Vector{Int}}) =
    Base.getindex(rl.lags, i)

#####################################################
# Simple immutable type to store embedding positions
#####################################################
"""
     EmbeddingPositions

Assume `r` is a `Embedding` instance.
Stores the information on which of the data series in
`r.which_dataseries`

"""
struct EmbeddingPositions
    positions::Vector{Int}
end

Base.getindex(rp::EmbeddingPositions,
            i::Union{Int, Colon, UnitRange, Vector{Int}}) =
    Base.getindex(rp.positions, i)

#####################################################
# Immutable type to store the data series that went
# into the embedding.
#####################################################
"""
    EmbeddingData{N, T}

The data and parameters that goes into a embedding.
Assume `rd` is a `EmbeddingData` instance. Then

* `rd.dataseries` gives a vector of dataseries (each itself being
a vector) that was used for the embedding.
* `rd.in_which_pos` indicates which of the dataseries appear in which
columns of the embedding. For example, `rd.in_which_pos[2] = 3` means
that the 3rd timeseries (which can accessed by `rd.dataseries[3]`)
appears in the second column of the embedding.
* `rd.at_what_lags[i]` gives the lags in the i-th column of the embedding.
* `rd.labels[i]` returns the label of the i-th dataseries
* `rd.descriptions[i]` returns the description of the i-th dataseries.
"""
struct EmbeddingData{N, T}
    dataseries::Vector{Vector{T}}
    labels::Vector{String}
    descriptions::Vector{String}
    in_which_pos::EmbeddingPositions
    at_what_lags::EmbeddingLags
end

#####################################################
# Embedding type.
#####################################################
"""
    Embedding{D, T}

An immutable type representing a state space embedding.
Assume `r` is a embedding, then

* `r.points` gives the embeded  state vectors
* `r.embeddingdata` returns the recon
* series that was used to make the embedding, ``

"""
struct Embedding{D, T} <: AbstractEmbedding{D, T}
    points::Dataset{D, T}
    embeddingdata::EmbeddingData
end

@inline Base.length(r::AbstractEmbedding) = length(r.points)

"""
    size(E::AbstractEmbedding) -> (Int, Int)

Get the size of the embedding, which is `(npoints, dim)`.
"""
@inline Base.size(r::AbstractEmbedding{D,T}) where {D,T} = (length(r.points), D)
@inline Base.size(r::AbstractEmbedding, i) = size(r.points)[i]

tps = Union{SVector{D, T} where {D, T}, Colon, UnitRange{Int}, AbstractVector{Int}}
@inline Base.getindex(r::AbstractEmbedding, i::Int) = r.points[i]
@inline Base.getindex(r::AbstractEmbedding, i::tps) = r.points[i]
@inline Base.getindex(r::AbstractEmbedding, i::Int, j::tps) = r.points[i, j]
@inline Base.getindex(r::AbstractEmbedding, i::tps, j::tps) = r.points[i, j]
@inline Base.getindex(r::AbstractEmbedding, i::Int, j::Colon) = r.points[i]
@inline Base.getindex(r::AbstractEmbedding, i::tps, j::Colon) = r.points[i]
@inline Base.getindex(r::AbstractEmbedding, i::Colon, j::Int) = r.points[i, j]
@inline Base.getindex(r::AbstractEmbedding, i::Colon, j::Colon) = r.points
@inline Base.getindex(r::AbstractEmbedding, i::Colon, j::tps) = r.points[i, j]

Base.unique(r::AbstractEmbedding) = Base.unique(r.points.data)
Base.unique(r::AbstractEmbedding, i::Int) = Base.unique(r.points.data, dims = i)

dimension(::AbstractEmbedding{D,T}) where {D,T} = D
@inline Base.eltype(r::AbstractEmbedding{D,T}) where {D,T} = Ti

import Base: ==
==(r₁::AbstractEmbedding, r₂::AbstractEmbedding) =
    r₁.points == r₂.points

"""
    npoints(E::AbstractEmbedding) -> Int

Get the dimension the number of points in the embedding.
"""
npoints(r::AbstractEmbedding) = length(r.points)

"""
    points(E::AbstractEmbedding) -> Int

Get the points of the embedding.
"""
points(r::AbstractEmbedding) = r.points

"""
    get_dataseries(E::AbstractEmbedding) -> Vector{Vector{T}} where {T<:Number}

Return the data series that was used for the embedding.
"""
get_dataseries(r::AbstractEmbedding) = r.embeddingdata.dataseries

"""
    get_dataseries(E::AbstractEmbedding) -> Vector{Vector{T}} where {T<:Number}

Return the i-th dataseries that was used for the embedding.
"""
get_dataseries(E::AbstractEmbedding, i::Int) = E.embeddingdata.dataseries[i]
series = get_dataseries

"""
    which_series(E::AbstractEmbedding, i::Int)

Return the index of the data series appear in the `i`-th column
of the embedding.

## Example
Assume `r` is a `Embedding` instance, then
`which_series(r, 2) = 1` indicates that the second column of
the embedding is populated by a lagged instance of the 1st
data series.
"""
which_series(E::AbstractEmbedding, i::Int) = E.embeddingdata.in_which_pos[i]

"""
    in_which_pos(E::AbstractEmbedding) -> EmbeddingPositions

Return the data series that went into the embedding.
"""
in_which_pos(E::AbstractEmbedding) = E.embeddingdata.in_which_pos

"""
    at_what_lags(E::AbstractEmbedding) -> Vector{Vector{T}} where {T<:Number}

Return the data series that went into the embedding.
"""
at_what_lags(E::AbstractEmbedding) = E.embeddingdata.at_what_lags

#####################################################
# Retrieve information about the embedding.
#####################################################
"""
    n(E::AbstractEmbedding) -> Int

How many points are there in the embedding?
"""
n(E::AbstractEmbedding) = length(E.points)

"""
    lags(E::AbstractEmbedding) -> EmbeddingLags

What is the embedding lag for each of the variables/columns
of the embedding?
"""
lags(E::AbstractEmbedding) = at_what_lags(E)

"""
    lag(E::AbstractEmbedding, i::Int) -> Int

What is the embedding lag for the i-th variable/column
of the embedding?
"""
lag(E::AbstractEmbedding, i::Int) = E.at_what_lags[i]

"""
    label(E::AbstractEmbedding, i::Int) -> String

Get the description label of the i-th time series that went
into the embedding.

"""
label(E::AbstractEmbedding, i::Int) = E.embeddingdata.labels[i]
"""
    labels(E::AbstractEmbedding) -> Vector{String}

Get the descriptions of the all time series that went
into the embedding.
"""
labels(E::AbstractEmbedding) = E.embeddingdata.labels

import Base.minimum, Base.maximum
import DelayEmbeddings.minima
import DelayEmbeddings.maxima

"""
    minimum(d::Dataset{D, T} where {D, T}) -> T

Finds the minimum value in a Dataset.
"""
minimum(d::Dataset{D, T}) where {D, T} = minimum(minima(d))

"""
    maximum(d::Dataset{D, T} where {D, T}) -> T

Finds the maximum value in a Dataset.
"""
maximum(d::Dataset{D, T}) where {D, T} = maximum(maxima(d))

"""
    minimum(r::Embedding{D, T}) where {D, T}) -> T

Finds the minimum value in a Dataset.
"""
minimum(r::AbstractEmbedding{D, T}) where {D, T} = minimum(r.points)

"""
    maximum(r::Dataset{D, T}) where {D, T}) -> T

Find the maximum value in a Dataset.
"""
maximum(r::AbstractEmbedding{D, T}) where {D, T} = maximum(r.points)

"""
    minima(r::Embedding) where {D, T}) -> `SVector{D, T}`

Returns an `SVector` containing the maxima of the embedding along
each axis.
"""
minima(r::AbstractEmbedding{D, T}) where {D, T} = minima(r.points)

"""
    maxima(r::Embedding{D, T}) where {D, T}) -> `SVector{D, T}`

Returns an `SVector` containing the maxima of the embedding along
each axis.
"""
maxima(r::AbstractEmbedding{D, T}) where {D, T} = maximum(r.points)

"""
    embed(ts::Vector{Vector{T}}, in_which_pos::Vector{Int},
                at_what_lags::Vector{Int}) where T <: Number ->
                Embedding

Perform a state space embedding of the vectors in `ts`.

## Arguments
1. `which_ts::Vector{Vector{T}} where T <: Number`. Contains the the time
    series to embed.
2. `in_which_pos::Vector{Int}``. The length of in_which_pos gives the dimension
    of the embedding. The value of the ith element of in_which_pos indicates
    which time series in the ith column of the embedding.
    - **Example 1**: if `which_ts = [ts1, ts2]`, then we index ts1 as 1 and
        ts2 as 2. Setting `in_which_pos = [2, 2, 1]` will result in a
        3-dimensional embedding where `ts2` will appear in columns 1 and 2,
        while `ts1` will appear in column 3.
    - **Example 2**: If `which_ts = [ts1, ts2, ts3]`, then
        `in_which_pos = [2, 1, 2, 3, 3]` results in a 5-dimensional embedding
        where `ts1`appears in column 2, `ts2` appears in columns 1 and 3, while
        `ts3`appears in columns 4 and 5.
3. `at_what_lags::Vector{Int}` sets the lag in each column. Must be the same
    length as `which_ts`.
    - **Example**: if `in_which_pos = [2, 2, 1]`, then
        `at_what_lags = [1, 0, -1]` means that the lag in column 1 is 1, the
        lag in the second column is 0 and the lag in the third column is -1.
"""
function embed(ts::Vector{Vector{T}},
               in_which_pos::Vector{Int},
               at_what_lags::Vector{Int};
                 labels::Vector{String} = ["" for x in 1:length(ts)]) where {T<:Number}
    dim = length(in_which_pos)
    minlag, maxlag = minimum(at_what_lags), maximum(at_what_lags)
    npts = length(ts[1]) - (maxlag + abs(minlag))
    E = zeros(T, dim, npts)

    for i in 1:length(in_which_pos)
        ts_ind = in_which_pos[i]
        TS = ts[ts_ind]
        lag = at_what_lags[i]

        if lag > 0
            E[i, :] = TS[((1 + abs(minlag)) + lag):(end - maxlag) + lag]
        elseif lag < 0
            E[i, :] = TS[((1 + abs(minlag)) - abs(lag)):(end - maxlag - abs(lag))]
        elseif lag == 0
            E[i, :] = TS[(1 + abs(minlag)):(end - maxlag)]
        end
    end
    N = length(ts)
    embeddingdata = EmbeddingData{N, T}(
        ts,  # the dataseries
        ["" for i = 1:length(ts)], # empty labels by default
        ["" for i = 1:length(ts)], # empty descriptions by default
        EmbeddingPositions(in_which_pos), # the positions in which the dataseries appear
        EmbeddingLags(at_what_lags) # the embedding lag for each column
        )
    Embedding(Dataset(transpose(E)), embeddingdata)
end

"""
    embed(A::AbstractArray)

Returns an embedding of a vector of vectors, treating each
vector as a dynamical variable. Zero lag is used for all the
columns.
"""
function embed(ts::Vector{Vector{T}}) where {T}
    D = length(ts)
    embed(ts, [i for i = 1:D], [0 for i in 1:D])
end

"""
    embed(A::AbstractArray)

Returns an embedding of an array, treating each
column as a dynamical variable. Zero lag is used
for all the columns.
"""
function embed(A::AbstractArray{T}) where {T}
    D = size(A, 2)
    embed([A[:, i] for i = 1:D],
            [i for i = 1:D],
            [0 for i in 1:D])
end

"""
    embed(d::AbstractArray{T, 2},
                    in_which_pos::Vector{Int},
                    at_what_lags::Vector{Int})

Embedding of data represented by an array. Each
column of the array must correspond to one data series.
"""
function embed(data::AbstractArray{T, 2},
                    in_which_pos::Vector{Int},
                    at_what_lags::Vector{Int}) where {T}
    D = size(data, 2)
    embed(
        [data[:, i] for i = 1:D],
        in_which_pos,
        at_what_lags
    )
end

Embedding(ts::Vector{Vector{T}},
               in_which_pos::Vector{Int},
               at_what_lags::Vector{Int}) where {T} =
               embed(ts, in_which_pos, at_what_lags)

function Embedding(v::Vector{Vector{T} where T})
    positions = [i for i in 1:length(v)]
    lags = [0 for i in 1:length(v)]
    embed(v, positions, lags)
end

"""
     embed(d::Dataset)

Returns a state space embedding of the column of `d`.
"""
function embed(d::Dataset)
    D = size(d, 2)
    embed([d[:, i] for i = 1:D])
end

"""
     embed(d::Dataset,
        in_which_pos::Vector{Int},
        at_what_lags::Vector{Int})

Returns a state space embedding of the column of `d`.

## Arguments
* `d::Dataset`: The columns of `d` contains the data series to use
    for the embedding.
* `in_which_pos::Vector{Int}``. The length of in_which_pos gives the dimension
    of the embedding. The value of the ith element of `in_which_pos`
    indicates which column of `d` goes in the i-th column of the embedding.
* `at_what_lags::Vector{Int}` sets the lag in each column of the reconstrution.
    Must be the same length as `dimension(d)`
    - **Example**: if `in_which_pos = [2, 2, 1]`, then
        `at_what_lags = [1, 0, -1]` means that the lag in column 1 is 1, the
        lag in the second column is 0 and the lag in the third column is -1.
"""
function embed(d::Dataset,
        in_which_pos::Vector{Int},
        at_what_lags::Vector{Int})
    embed(
        [d[:, i] for i = 1:DelayEmbeddings.dimension(d)],
        in_which_pos,
        at_what_lags
    )
end

function Embedding(d::Dataset,
        in_which_pos::Vector{Int},
        at_what_lags::Vector{Int})
    embed(
        [d[:, i] for i = 1:DelayEmbeddings.dimension(d)],
        in_which_pos,
        at_what_lags
    )
end

########################################
# SArrays and SVectors
########################################
function embed(ts::Vector{SVector{N, T}},
        in_which_pos::Vector{Int},
        at_what_lags::Vector{Int};
        labels::Vector{String} =  ["" for x in 1:length(ts)]) where {N, T}
    dim = length(in_which_pos)
    minlag, maxlag = minimum(at_what_lags), maximum(at_what_lags)
    npts = length(ts[1]) - (maxlag + abs(minlag))
    E = zeros(T, dim, npts)

    for i in 1:length(in_which_pos)
        ts_ind = in_which_pos[i]
        TS = ts[ts_ind]
        lag = at_what_lags[i]

        if lag > 0
            E[i, :] = TS[((1 + abs(minlag)) + lag):(end - maxlag) + lag]
        elseif lag < 0
            E[i, :] = TS[((1 + abs(minlag)) - abs(lag)):(end - maxlag - abs(lag))]
        elseif lag == 0
            E[i, :] = TS[(1 + abs(minlag)):(end - maxlag)]
        end
    end
    embeddingdata = EmbeddingData{N, T}(
        ts,  # the dataseries
        ["" for i = 1:length(ts)], # empty labels by default
        ["" for i = 1:length(ts)], # empty descriptions by default
        EmbeddingPositions(in_which_pos), # the positions in which the dataseries appear
        EmbeddingLags(at_what_lags) # the embedding lag for each column
        )
    Embedding(Dataset(Array(transpose(E))), embeddingdata)
end

embed(A::SArray) = embed([A[:, i] for i in 1:size(S, 2)],
                        [i for i in 1:size(A, 2)],
                        [0 for i in 1:size(A, 2)])

embed(A::SArray, in_which_pos::Vector{Int}, at_what_lags::Vector{Int}) =
    embed([A[:, i] for i in 1:size(A, 2)],
    in_which_pos, at_what_lags)

####################################
# Pretty printing.
####################################
function summarise(r::EmbeddingData)
    n  = length(r.dataseries)
    _type = typeof(r)
    """$_type consisting of $n data series."""
end

function description(r::EmbeddingData)
    S = String[]
    lags_s = string.(r.at_what_lags)
    positions_s = join(string.(r.in_which_pos))
    indexed_s = join(string.(1:length(r.dataseries)), " ")
    if all(r.labels .== "")

        push!(S, " The time series are both unnamed")
    else
        lbls = join([" The time series have labels", r.labels], ", ")
        push!(S, lbls)
    end
    push!(S, ", are indexed $indexed_s, and the data series appear")
    push!(S, " as columns of the embedding in $positions_s at $lags_s.")
    S = join([summarise(r); S], "")
    return S
end

using Base.Iterators: flatten

function matstring(d::AbstractDataset{D, T}) where {D, T}
    N = length(d)
    if N > 36
        mat = zeros(eltype(d), 40, D)
        for (i, a) in enumerate(flatten((1:20, N-19:N)))
            mat[i, :] .= d[a]
        end
    else
        mat = Matrix(d)
    end
    s = sprint(io -> show(IOContext(io, :limit=>true), MIME"text/plain"(), mat))
    s = join(split(s, '\n')[2:end], '\n')
    return s
end

function summarise(r::AbstractEmbedding)
    n_dataseries = length(r.embeddingdata.dataseries)
    embedding_type = typeof(r)
    npts = length(r.points)
    summary = "$embedding_type with $npts points\n"
    return join([summary, matstring(r.points)], "")
end

Base.show(io::IO, r::EmbeddingData) = println(io, description(r))
Base.show(io::IO, r::Embedding) = println(io, summarise(r))

export
AbstractEmbedding,
Embedding,
EmbeddingData,
EmbeddingLags,
EmbeddingPositions,
embed,
label, labels,
lag, lags,
n,
at_what_lags,
in_which_pos,
which_series,
get_dataseries,
points,
npoints,
dimension,
minimum, maximum,
minima, maxima


text:
using LLVM.Interop

@testset "interop" begin

@testset "base" begin

@test isa(JuliaContext(), Context)

@test isa(convert(LLVMType, Nothing), LLVM.VoidType)

@test_throws ErrorException convert(LLVMType, Ref)
convert(LLVMType, Ref, true)

@generated function add_one(i)
    T_int = convert(LLVMType, Int)

    f, ft = create_function(T_int, [T_int])

    Builder(JuliaContext()) do builder
        entry = BasicBlock(f, "entry", JuliaContext())
        position!(builder, entry)

        val = add!(builder, parameters(f)[1], ConstantInt(T_int, 1))

        ret!(builder, val)
    end

    call_function(f, Int, Tuple{Int}, :((i,)))
end

@test add_one(1) == 2

@eval struct GhostType end
@eval struct NonGhostType1
    x::Int
end
@eval mutable struct NonGhostType2 end

@test isghosttype(GhostType)
@test !isghosttype(NonGhostType1)
@test !isghosttype(NonGhostType2)
@test isboxed(NonGhostType2)

end

@testset "asmcall" begin

# only asm

a1() = @asmcall("nop")
@test a1() == nothing

a2() = @asmcall("nop", Nothing)
@test a2() == nothing

a3() = @asmcall("nop", Nothing, Tuple{})
@test a3() == nothing

# asm + constraints

b1() = @asmcall("nop", "")
@test b1() == nothing

b2() = @asmcall("nop", "", Nothing)
@test b2() == nothing

b3() = @asmcall("nop", "", Nothing, Tuple{})
@test b3() == nothing

# asm + constraints + side-effects

c1() = @asmcall("nop", "", false)
@test c1() == nothing

c2() = @asmcall("nop", "", false, Nothing)
@test c2() == nothing

c3() = @asmcall("nop", "", false, Nothing, Tuple{})
@test c3() == nothing

if Sys.ARCH == :x86 || Sys.ARCH == :x86_64

# arguments

d1(a) = @asmcall("bswap \$0", "=r,r", Int32, Tuple{Int32}, a)
@test d1(Int32(1)) == Int32(16777216)

# multiple output registers

e1() = @asmcall("mov \$\$1, \$0; mov \$\$2, \$1;", "=r,=r", Tuple{Int16,Int32})
@test e1() == (Int16(1), Int32(2))

# TODO: alternative test snippets for other platforms

end

end

if VERSION >= v"1.2.0-DEV.531"
@testset "passes" begin

Context() do ctx
LLVM.Module("SomeModule", ctx) do mod
ModulePassManager() do pm

alloc_opt!(pm)
barrier_noop!(pm)
gc_invariant_verifier!(pm)
gc_invariant_verifier!(pm, true)
lower_exc_handlers!(pm)
combine_mul_add!(pm)
multi_versioning!(pm)
propagate_julia_addrsp!(pm)
lower_ptls!(pm)
lower_ptls!(pm, true)
lower_simdloop!(pm)
late_lower_gc_frame!(pm)
final_lower_gc!(pm)

end
end
end

end
end

end

ext: jl


text:
@test begin
    l = Segment(Point(0, 0), Point(1, 1))
    L = Segment(Point(2, 0), Point(3, 1))
    !collides(l, L)
end

@test begin
    l = Segment(Point(0, 0), Point(1, 1))
    L = Segment(Point(1, 0), Point(2, 1))
    !collides(l, L)
end

@test begin
    l = Segment(Point(0, 0), Point(1, 1))
    L = Segment(Point(0.75, 0), Point(1, 0.25))
    !collides(l, L)
end

@test begin
    l = Segment(Point(0, 0), Point(1, 0))
    L = Segment(Point(0, 1), Point(1, 1))
    !collides(l, L)
end

@test begin
    l = Segment(Point(0, 0), Point(1, 0))
    L = Segment(Point(2, 0), Point(3, 0))
    !collides(l, L)
end

@test collides(Segment(Point(-1, 0), Point(0, 0)), Segment(Point(0, 0), Point(1, 0)))

@test begin
    l = Segment(Point(0, 0), Point(1, 0))
    L = Segment(Point(1, 0), Point(3, 0))
    collides(l, L)
end

@test begin
    l = Segment(Point(0, 0), Point(1, 1))
    L = Segment(Point(0, 1), Point(1, 0))
    collides(l, L)
end

@test begin
    l = Segment(Point(0, 0), Point(1, 1))
    collides(l, l)
end

# Box vs segment
@test begin
    b = Box(Point(0, 0), Point(1, 1))
    l = Segment(Point(2, 2), Point(3, 3))
    !collides(l, b)
end

@test begin
    b = Box(Point(0, 0), Point(1, 1))
    l = Segment(Point(1, 1), Point(3, 3))
    collides(l, b)
end

@test begin
    b = Box(Point(0, 0), Point(1, 1))
    l = Segment(Point(1, 1), Point(0, 1))
    collides(l, b)
end

@test begin
    b = Box(Point(0, 0), Point(1, 1))
    l = Segment(Point(1, 1), Point(0, 0))
    collides(l, b)
end

@test begin
    b = Box(Point(0, 0), Point(1, 1))
    l = Segment(Point(0.25, 0.25), Point(0.75, 0.75))
    collides(l, b)
end

# -------------------------------------------------------------------
# # Both in
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0.25, 0.25, 0.25), Point(0.75, 0.75, 0.75))
    collides(l, b)
end

# One in other out
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0.5, 0.5, 0.5), Point(3, 3, 3))
    collides(l, b)
end

# Both out but cross
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0.5, 0.5, -1), Point(0.5, 0.5, 3))
    collides(l, b)
end

# Both out but cross
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0.25, 0.25, -1), Point(0.75, 0.75, 3))
    collides(l, b)
end

# Both out but cross
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0.123, 0.1, -1), Point(0.567, 0.456, 3))
    collides(l, b)
end

# Point collision
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(1, 1, 1), Point(3, 3, 3))
    collides(l, b)
end

# Both out
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(2, 2, 2), Point(3, 3, 3))
    !collides(l, b)
end

# Both out
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0.5, 0.5, 0), Point(0.5, 0.5, -2))
    collides(l, b)
end

# Entire segment collision
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0, 0, 0), Point(1, 0, 0))
    collides(l, b)
end

# Intersect in plane
@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    l = Segment(Point(0, 0, 0), Point(1, 1, 0))
    collides(l, b)
end

max_stars_repo_name: MiroK/fractional-dim


text:
module TSeriesForecast

submodules = [
    "simpleexponentialsmoothing",
    "holts_trend_method",
    "seasonality_exponential_smoothing"
]

for sm ∈ submodules
    include("./$sm.jl")
end

end


text:
module BenchFindall

using BenchmarkTools
using Transducers

needle = 123.0
xs = rand(100_000)
xs[xs.>0.001] .= needle

suite = BenchmarkGroup()
suite["base"] = @benchmarkable findall(isequal($needle), $xs)
suite["xf-array"] = @benchmarkable collect(
    Enumerate() |> Filter(!isequal($needle) ∘ last) |> Map(first),
    $xs,
)
suite["xf-iter"] = @benchmarkable collect(
    Enumerate() |> Filter(isequal($needle) ∘ last) |> Map(first),
    Transducers.NoAdjoint(x for x in $xs),
)

end  # module
BenchFindall.suite


text:
include("header.jl")

@testset "iterate" begin
    o = (:delta => 0.01,)

    function itr1(w)
        total = 0.0
        for wi in w; total+=sum(wi); end
        return total
    end

    function itr2(w)
        total = 0.0
        for (i,wi) in enumerate(w); total+=sum(wi); end
        return total
    end

    function itr3(w)
        total = 0.0
        for (k,wi) in w; total+=sum(wi); end
        return total
    end

    function itr4(w)
        w1,w2 = w
        return sum(w1)+sum(w2)
    end

    for T in (Float32, Float64)
        warray  = [randn(T,2,3,5),randn(T,1,3,5)]
        wtuple  = (randn(T,2,3,5),randn(T,1,3,5))
        wdict   = Dict(:w1=>randn(T,2,3,5),:w2=>randn(T,1,3,5))

        @test gradcheck(itr1,warray; o...)
        @test gradcheck(itr2,warray; o...)
        @test gradcheck(itr1,warray[1]; o...)
        @test gradcheck(itr2,warray[1]; o...)
        @test gradcheck(itr3,wdict; o...)
        @test gradcheck(itr4,warray; o...)
        @test gradcheck(itr4,wtuple; o...)
    end
end


text:
using Documenter, TestImages, Images, OMETIFF, TOML

include("generate_imagelist.jl")

generate_imagelist(joinpath(@__DIR__, "src"))

format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true")

makedocs(;
    modules=[TestImages],
    format=format,
    pages=[
        "Home" => "index.md",
        "List of images" => "imagelist.md",
        "Function reference" => "functionreference.md",
        "Contributing" => "contributing.md",
    ],
    repo="https://github.com/JuliaImages/TestImages.jl/blob/{commit}{path}#L{line}",
    sitename="TestImages.jl",
)

deploydocs(;
    repo="github.com/JuliaImages/TestImages.jl",
    # A lot of thumbnails are generated, so compressing history in `gh-pages` avoids infinite
    # storage increasement.
    forcepush=true,
)


text:
# Command: Delete Cookie
"""
    delete!(session::Session, cookie::AbstractString)::Nothing

If cookie = "", delete all cookies. Otherwise, delete cookie matching the "cookie" name.
"""
function delete!(session::Session, cookie::AbstractString)::Nothing
    @unpack addr, id = session
    if isempty(cookie)
        response = HTTP.delete(
            "$addr/session/$id/cookie",
            [("Content-Type" => "application/json; charset=utf-8")],
        )
        @assert response.status == 200
    else
        response = HTTP.delete(
            "$addr/session/$id/cookie/$(escapeuri(cookie))",
            [("Content-Type" => "application/json; charset=utf-8")],
        )
        @assert response.status == 200
    end
    nothing
end


text:
@testset "Parameter" begin
    p = GenericPara(diagType = Ver4Diag, innerLoopNum = 1)
    q = GenericPara(diagType = Ver4Diag, innerLoopNum = 2)
    a = GenericPara(diagType = Ver4Diag, innerLoopNum = 2)

    @test p != q
    @test q == a

    aa = reconstruct(a, transferLoop = [0.0, 0.0, 0.0])
    @test a != aa

    # reconstruct with the same diagType but different interaction
    aaa = reconstruct(a, interaction = [])
    @test a != aaa

    # reconstruct with different diagram type leads to different parameter
    #reconstructed GenericPara uses the old parameters such as firstLoopIdx, totalLoopNum etc., which are different between diagTypes
    s = reconstruct(a, diagType = SigmaDiag)
    ss = GenericPara(diagType = SigmaDiag, innerLoopNum = 2)
    @test s != ss

end


text:
using AlphaStructures
using ViewerGL
GL =  ViewerGL

"""
Return points from file.
"""
function load_points(filename::String)::Matrix
    io = open(filename, "r")
    point = readlines(io)
    close(io)

    b = [tryparse.(Float64,split(point[i], " ")) for i in 1:length(point)]
    V = hcat(b...)
    return V
end

filename = "examples/examples2D/model/lar_text.txt";
V = load_points(filename);

# view points
GL.VIEW([
	GL.GLPoints(permutedims(V))
])

# points by columns
filtration = AlphaStructures.alphaFilter(V);
VV,EV,FV = AlphaStructures.alphaSimplex(V, filtration, 0.01)

# View triangle
GL.VIEW([
	GL.GLGrid(V, FV, GL.COLORS[1], 0.8)
])


text:
using  Literate
using  Plots

files =  filter( f -> startswith(f, "0"), readdir("src")) |> collect

for file in files
    Literate.notebook("src/$file", "notebooks",  execute=false)
end


text:
@with_kw struct SimpleLightDark <: POMDP{Int,Int,Float64}
    discount::Float64       = 0.95
    correct_r::Float64      = 100.0
    incorrect_r::Float64    = -100.0
    light_loc::Int          = 10
    radius::Int             = 60
end
discount(p::SimpleLightDark) = p.discount
isterminal(p::SimpleLightDark, s::Number) = !(s in -p.radius:p.radius)

const ACTIONS = [-10, -1, 0, 1, 10]
actions(p::SimpleLightDark) = ACTIONS
n_actions(p::SimpleLightDark) = length(actions(p))
const ACTION_INDS = Dict(a=>i for (i,a) in enumerate(actions(SimpleLightDark())))
actionindex(p::SimpleLightDark, a::Int) = ACTION_INDS[a]

states(p::SimpleLightDark) = -p.radius:p.radius + 1
n_states(p::SimpleLightDark) = length(states(p))
stateindex(p::SimpleLightDark, s::Int) = s+p.radius+1

function transition(p::SimpleLightDark, s::Int, a::Int)
    if a == 0
        return SparseCat(SVector(p.radius+1), SVector(1.0))
    else
        return SparseCat(SVector(clamp(s+a, -p.radius, p.radius)), SVector(1.0))
    end
end

observation(p::SimpleLightDark, sp) = Normal(sp, abs(sp - p.light_loc) + 0.0001)

function reward(p::SimpleLightDark, s, a)
    if a == 0
        return s == 0 ? p.correct_r : p.incorrect_r
    else
        return -1.0
    end
end

function initialstate(p::SimpleLightDark)
    ps = ones(2*div(p.radius,2)+1)
    ps /= length(ps)
    return SparseCat(div(-p.radius,2):div(p.radius,2), ps)
end

@with_kw struct DSimpleLightDark <: POMDP{Int, Int, Int}
    sld::SimpleLightDark       = SimpleLightDark()
    binsize::Float64           = 1.0
    _d_disc::Vector{SparseCat} = SparseCat[get_discrete_obs_dist(sld, sp, binsize) for sp in -sld.radius:sld.radius]
end

sld(s::SimpleLightDark) = s
sld(s::DSimpleLightDark) = s.sld

function get_discrete_obs_dist(p::SimpleLightDark, sp::Int, binsize::Float64)
    dist = observation(p, sp)
    # discretize observations
    interval_start = 0.0
    discrete_dist = Vector{Float64}(undef, floor(Int, div(p.radius*2, binsize)) + 2)
    @assert div(p.radius*2, binsize) + 2 == length(range(-p.radius, step=binsize, stop=p.radius)) + 1
    for (oi, o) in enumerate(range(-p.radius, step=binsize, stop=p.radius))
        interval_end = cdf(dist, o)
        discrete_dist[oi] = interval_end - interval_start
        interval_start = interval_end
    end
    discrete_dist[end] = 1.0 - interval_start
    return SparseCat(1:length(discrete_dist), discrete_dist)
end

observations(p::DSimpleLightDark) = 1:(div(p.sld.radius * 2, p.binsize) + 2)
n_observations(p::DSimpleLightDark) = length(observations(p))
obsindex(p::DSimpleLightDark, o::Int) = o

actions(p::DSimpleLightDark) = actions(p.sld)
n_actions(p::DSimpleLightDark) = n_actions(p.sld)
actionindex(p::DSimpleLightDark, a::Int) = actionindex(p.sld, a)

states(p::DSimpleLightDark) = states(p.sld)
n_states(p::DSimpleLightDark) = n_states(p.sld)
stateindex(p::DSimpleLightDark, s::Int) = stateindex(p.sld, s)

observation(p::DSimpleLightDark, sp::Int) = p._d_disc[stateindex(p, sp)]
transition(p::DSimpleLightDark, s::Int, a::Int) = transition(p.sld, s, a)
reward(p::DSimpleLightDark, s, a) = reward(p.sld, s, a)
initialstate(p::DSimpleLightDark) = initialstate(p.sld)

discount(p::DSimpleLightDark) = discount(p.sld)
isterminal(p::DSimpleLightDark, s::Number) = isterminal(p.sld, s)


text:
module Operations

using Base.Random: UUID
using Base: LibGit2
using Pkg3.TerminalMenus
using Pkg3.Types
import Pkg3: Pkg2, depots, BinaryProvider, USE_LIBGIT2_FOR_ALL_DOWNLOADS, NUM_CONCURRENT_DOWNLOADS

const SlugInt = UInt32 # max p = 4
const chars = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789"
const nchars = SlugInt(length(chars))
const max_p = floor(Int, log(nchars, typemax(SlugInt) >>> 8))

function slug(x::SlugInt, p::Int)
    1 ≤ p ≤ max_p || # otherwise previous steps are wrong
        error("invalid slug size: $p (need 1 ≤ p ≤ $max_p)")
    return sprint() do io
        for i = 1:p
            x, d = divrem(x, nchars)
            write(io, chars[1+d])
        end
    end
end
slug(x::Integer, p::Int) = slug(SlugInt(x), p)

function slug(bytes::Vector{UInt8}, p::Int)
    n = nchars^p
    x = zero(SlugInt)
    for (i, b) in enumerate(bytes)
        x = (x + b*powermod(2, 8(i-1), n)) % n
    end
    slug(x, p)
end

slug(uuid::UUID, p::Int=4) = slug(uuid.value % nchars^p, p)
slug(sha1::SHA1, p::Int=4) = slug(sha1.bytes, p)

version_slug(uuid::UUID, sha1::SHA1) = joinpath(slug(uuid), slug(sha1))

function find_installed(uuid::UUID, sha1::SHA1)
    slug = version_slug(uuid, sha1)
    for depot in depots()
        path = abspath(depot, "packages", slug)
        ispath(path) && return path
    end
    return abspath(depots()[1], "packages", slug)
end

function package_env_info(pkg::String, env::EnvCache = EnvCache(); verb::String = "choose")
    project = env.project
    manifest = env.manifest
    haskey(manifest, pkg) || return nothing
    infos = manifest[pkg]
    isempty(infos) && return nothing
    if haskey(project, "deps") && haskey(project["deps"], pkg)
        uuid = project["deps"][pkg]
        filter!(infos) do info
            haskey(info, "uuid") && info["uuid"] == uuid
        end
        length(infos) < 1 &&
            cmderror("manifest has no stanza for $pkg/$uuid")
        length(infos) > 1 &&
            cmderror("manifest has multiple stanzas for $pkg/$uuid")
        return first(infos)
    elseif length(infos) == 1
        return first(infos)
    else
        options = String[]
        paths = convert(Dict{String,Vector{String}}, find_registered(pkg))
        for info in infos
            uuid = info["uuid"]
            option = uuid
            if haskey(paths, uuid)
                for path in paths[uuid]
                    info′ = parse_toml(path, "package.toml")
                    option *= " – $(info′["repo"])"
                    break
                end
            else
                option *= " – (unregistred)"
            end
            push!(options, option)
        end
        menu = RadioMenu(options)
        choice = request("Which $pkg package do you want to use:", menu)
        choice == -1 && cmderror("Package load aborted")
        return infos[choice]
    end
end

get_or_make(::Type{T}, d::Dict{K}, k::K) where {T,K} =
    haskey(d, k) ? convert(T, d[k]) : T()

function load_versions(path::String)
    toml = parse_toml(path, "versions.toml")
    Dict(VersionNumber(ver) => SHA1(info["hash-sha1"]) for (ver, info) in toml)
end

function load_package_data(f::Base.Callable, path::String, versions)
    toml = parse_toml(path, fakeit=true)
    data = Dict{VersionNumber,Dict{String,Any}}()
    for ver in versions
        ver::VersionNumber
        for (v, d) in toml, (key, value) in d
            vr = VersionRange(v)
            ver in vr || continue
            dict = get!(data, ver, Dict{String,Any}())
            haskey(dict, key) && cmderror("$ver/$key is duplicated in $path")
            dict[key] = f(value)
        end
    end
    return data
end

load_package_data(f::Base.Callable, path::String, version::VersionNumber) =
    get(load_package_data(f, path, [version]), version, nothing)

function deps_graph(env::EnvCache, pkgs::Vector{PackageSpec})
    deps = Dict{UUID,Dict{VersionNumber,Tuple{SHA1,Dict{UUID,VersionSpec}}}}()
    uuids = [pkg.uuid for pkg in pkgs]
    seen = UUID[]
    while true
        unseen = setdiff(uuids, seen)
        isempty(unseen) && break
        for uuid in unseen
            push!(seen, uuid)
            deps[uuid] = valtype(deps)()
            for path in registered_paths(env, uuid)
                version_info = load_versions(path)
                versions = sort!(collect(keys(version_info)))
                dependencies = load_package_data(UUID, joinpath(path, "dependencies.toml"), versions)
                compatibility = load_package_data(VersionSpec, joinpath(path, "compatibility.toml"), versions)
                for (v, h) in version_info
                    d = get_or_make(Dict{String,UUID}, dependencies, v)
                    r = get_or_make(Dict{String,VersionSpec}, compatibility, v)
                    q = Dict(u => get_or_make(VersionSpec, r, p) for (p, u) in d)
                    # VERSION in get_or_make(VersionSpec, r, "julia") || continue
                    deps[uuid][v] = (h, q)
                    for (p, u) in d
                        u in uuids || push!(uuids, u)
                    end
                end
            end
        end
        find_registered!(env, uuids)
    end
    return deps
end

"Resolve a set of versions given package version specs"
function resolve_versions!(env::EnvCache, pkgs::Vector{PackageSpec})::Dict{UUID,VersionNumber}
    info("Resolving package versions")
    # anything not mentioned is fixed
    uuids = UUID[pkg.uuid for pkg in pkgs]
    uuid_to_name = Dict{String, String}()
    for (name::String, uuid::UUID) in env.project["deps"]
        uuid_to_name[string(uuid)] = name
        uuid in uuids && continue
        info = manifest_info(env, uuid)
        haskey(info, "version") || continue
        ver = VersionNumber(info["version"])
        push!(pkgs, PackageSpec(name, uuid, ver))
    end
    # construct data structures for resolver and call it
    reqs = Dict{String,Pkg2.Types.VersionSet}(string(pkg.uuid) => pkg.version for pkg in pkgs)
    deps = convert(Dict{String,Dict{VersionNumber,Pkg2.Types.Available}}, deps_graph(env, pkgs))
    for dep_uuid in keys(deps)
        info = manifest_info(env, UUID(dep_uuid))
        if info != nothing
            uuid_to_name[info["uuid"]] = info["name"]
        end
    end
    deps = Pkg2.Query.prune_dependencies(reqs, deps, uuid_to_name)
    vers = convert(Dict{UUID,VersionNumber}, Pkg2.Resolve.resolve(reqs, deps, uuid_to_name))
    find_registered!(env, collect(keys(vers)))
    # update vector of package versions
    for pkg in pkgs
        pkg.version = vers[pkg.uuid]
    end
    uuids = UUID[pkg.uuid for pkg in pkgs]
    for (uuid, ver) in vers
        uuid in uuids && continue
        name = registered_name(env, uuid)
        push!(pkgs, PackageSpec(name, uuid, ver))
    end
    return vers
end

"Find names, repos and hashes for each package UUID & version"
function version_data(env::EnvCache, pkgs::Vector{PackageSpec})
    names = Dict{UUID,String}()
    hashes = Dict{UUID,SHA1}()
    upstreams = Dict{UUID,Vector{String}}()
    for pkg in pkgs
        uuid = pkg.uuid
        ver = pkg.version::VersionNumber
        upstreams[uuid] = String[]
        for path in registered_paths(env, uuid)
            info = parse_toml(path, "package.toml")
            if haskey(names, uuid)
                names[uuid] == info["name"] ||
                    cmderror("$uuid: name mismatch between registries: ",
                             "$(names[uuid]) vs. $(info["name"])")
            else
                names[uuid] = info["name"]
            end
            repo = info["repo"]
            repo in upstreams[uuid] || push!(upstreams[uuid], repo)
            vers = load_versions(path)
            if haskey(vers, ver)
                h = vers[ver]
                if haskey(hashes, uuid)
                    h == hashes[uuid] ||
                        warn("$uuid: hash mismatch for version $ver!")
                else
                    hashes[uuid] = h
                end
            end
        end
        @assert haskey(hashes, uuid)
    end
    foreach(sort!, values(upstreams))
    return names, hashes, upstreams
end

const refspecs = ["+refs/*:refs/remotes/cache/*"]

function get_archive_url_for_version(url::String, version)
    if (m = match(r"https://github.com/(.*?)/(.*?).git", url)) != nothing
        return "https://github.com/$(m.captures[1])/$(m.captures[2])/archive/v$(version).tar.gz"
    end
    return nothing
end

function install(
    env::EnvCache,
    uuid::UUID,
    name::String,
    hash::SHA1,
    urls::Vector{String},
    version::Union{VersionNumber,Void} = nothing
)::Tuple{String,Bool}
    # returns path to version & if it's newly installed
    version_path = find_installed(uuid, hash)
    ispath(version_path) && return version_path, false
    http_download_successful = false
    env.preview[] && return version_path, true
    if !USE_LIBGIT2_FOR_ALL_DOWNLOADS && version != nothing
        for url in urls
            archive_url = get_archive_url_for_version(url, version)
            if archive_url != nothing
                path = joinpath(tempdir(), name * "_" * randstring(6) * ".tar.gz")
                url_success = true
                try
                    cmd = BinaryProvider.gen_download_cmd(archive_url, path);
                    run(cmd, (DevNull, DevNull, DevNull))
                catch e
                    e isa InterruptException && rethrow(e)
                    url_success = false
                end
                url_success || continue
                http_download_successful = true
                dir = joinpath(tempdir(), randstring(12))
                mkpath(dir)
                cmd = BinaryProvider.gen_unpack_cmd(path, dir);
                run(cmd, (DevNull, DevNull, DevNull))
                dirs = readdir(dir)
                # 7z on Win might create this spurious file
                filter!(x -> x != "pax_global_header", dirs)
                @assert length(dirs) == 1
                !isdir(version_path) && mkpath(version_path)
                mv(joinpath(dir, dirs[1]), version_path; remove_destination=true)
                Base.rm(path; force = true)
                break # object was found, we can stop
            end
        end
    end
    if !http_download_successful || USE_LIBGIT2_FOR_ALL_DOWNLOADS
        upstream_dir = joinpath(depots()[1], "upstream")
        ispath(upstream_dir) || mkpath(upstream_dir)
        repo_path = joinpath(upstream_dir, string(uuid))
        repo = ispath(repo_path) ? LibGit2.GitRepo(repo_path) : begin
            # info("Cloning [$uuid] $name")
            LibGit2.clone(urls[1], repo_path, isbare=true)
        end
        git_hash = LibGit2.GitHash(hash.bytes)
        for i = 2:length(urls)
            try LibGit2.GitObject(repo, git_hash)
                break # object was found, we can stop
            catch err
                err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
            end
            url = urls[i]
            LibGit2.fetch(repo, remoteurl=url, refspecs=refspecs)
        end
        tree = try
            LibGit2.GitObject(repo, git_hash)
        catch err
            err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow(err)
            error("$name: git object $(string(hash)) could not be found")
        end
        tree isa LibGit2.GitTree ||
            error("$name: git object $(string(hash)) should be a tree, not $(typeof(tree))")
        mkpath(version_path)
        opts = LibGit2.CheckoutOptions(
            checkout_strategy = LibGit2.Consts.CHECKOUT_FORCE,
            target_directory = Base.unsafe_convert(Cstring, version_path)
        )
        h = string(hash)[1:16]
        LibGit2.checkout_tree(repo, tree, options=opts)
    end
    return version_path, true
end

function update_manifest(env::EnvCache, uuid::UUID, name::String, hash::SHA1, version::VersionNumber)
    infos = get!(env.manifest, name, Dict{String,Any}[])
    info = nothing
    for i in infos
        UUID(i["uuid"]) == uuid || continue
        info = i
        break
    end
    if info == nothing
        info = Dict{String,Any}("uuid" => string(uuid))
        push!(infos, info)
    end
    info["version"] = string(version)
    info["hash-sha1"] = string(hash)
    delete!(info, "deps")
    for path in registered_paths(env, uuid)
        data = load_package_data(UUID, joinpath(path, "dependencies.toml"), version)
        data == nothing && continue
        info["deps"] = convert(Dict{String,String}, data)
        break
    end
    return info
end

function prune_manifest(env::EnvCache)
    keep = map(UUID, values(env.project["deps"]))
    while !isempty(keep)
        clean = true
        for (name, infos) in env.manifest, info in infos
            haskey(info, "uuid") && haskey(info, "deps") || continue
            UUID(info["uuid"]) ∈ keep || continue
            for dep::UUID in values(info["deps"])
                dep ∈ keep && continue
                push!(keep, dep)
                clean = false
            end
        end
        clean && break
    end
    if VERSION < v"0.7.0-DEV.1393"
        filter!(env.manifest) do _, infos
            filter!(infos) do info
                haskey(info, "uuid") && UUID(info["uuid"]) ∈ keep
            end
            !isempty(infos)
        end
    else
        filter!(env.manifest) do _infos # (_, info) doesn't parse on 0.6
            _, infos = _infos
            filter!(infos) do info
                haskey(info, "uuid") && UUID(info["uuid"]) ∈ keep
            end
            !isempty(infos)
        end
    end
end

function apply_versions(env::EnvCache, pkgs::Vector{PackageSpec})::Vector{UUID}
    names, hashes, urls = version_data(env, pkgs)
    # install & update manifest
    new_versions = UUID[]

    jobs = Channel(NUM_CONCURRENT_DOWNLOADS);
    results = Channel(NUM_CONCURRENT_DOWNLOADS);
    @schedule begin
        for pkg in pkgs
            put!(jobs, pkg)
        end
    end

    for i in 1:NUM_CONCURRENT_DOWNLOADS
        @schedule begin
            for pkg in jobs
                uuid = pkg.uuid
                version = pkg.version::VersionNumber
                name, hash = names[uuid], hashes[uuid]
                try
                    version_path, new = install(env, uuid, name, hash, urls[uuid], version)
                    put!(results, (pkg, version_path, version, hash, new))
                catch e
                    put!(results, e)
                end
            end
        end
    end

    textwidth = VERSION < v"0.7.0-DEV.1930" ? Base.strwidth : Base.textwidth
    max_name = maximum(textwidth(names[pkg.uuid]) for pkg in pkgs)

    for _ in 1:length(pkgs)
        r = take!(results)
        r isa Exception && cmderror("Error when installing packages:\n", sprint(Base.showerror, r))
        pkg, path, version, hash, new = r
        if new
            vstr = version != nothing ? "v$version" : "[$h]"
            new && info("Installed $(rpad(names[pkg.uuid] * " ", max_name + 2, "─")) $vstr")
        end
        uuid = pkg.uuid
        version = pkg.version::VersionNumber
        name, hash = names[uuid], hashes[uuid]
        update_manifest(env, uuid, name, hash, version)
        new && push!(new_versions, uuid)
    end
    prune_manifest(env)
    return new_versions
end

function dependency_order_uuids(env::EnvCache, uuids::Vector{UUID})::Dict{UUID,Int}
    order = Dict{UUID,Int}()
    seen = UUID[]
    k::Int = 0
    function visit(uuid::UUID)
        uuid in seen &&
            return warn("Dependency graph not a DAG, linearizing anyway")
        haskey(order, uuid) && return
        push!(seen, uuid)
        info = manifest_info(env, uuid)
        haskey(info, "deps") &&
            foreach(visit, values(info["deps"]))
        pop!(seen)
        order[uuid] = k += 1
    end
    visit(uuid::String) = visit(UUID(uuid))
    foreach(visit, uuids)
    return order
end

function build_versions(env::EnvCache, uuids::Vector{UUID})
    # collect builds for UUIDs with `deps/build.jl` files
    env.preview[] && (info("Skipping building in preview mode"); return)
    builds = Tuple{UUID,String,SHA1,String}[]
    for uuid in uuids
        info = manifest_info(env, uuid)
        name = info["name"]
        # TODO: handle development packages?
        haskey(info, "hash-sha1") || continue
        hash = SHA1(info["hash-sha1"])
        path = find_installed(uuid, hash)
        ispath(path) || error("Build path for $name does not exist: $path")
        build_file = joinpath(path, "deps", "build.jl")
        ispath(build_file) && push!(builds, (uuid, name, hash, build_file))
    end
    # toposort builds by dependencies
    order = dependency_order_uuids(env, map(first, builds))
    sort!(builds, by = build -> order[first(build)])
    # build each package verions in a child process
    withenv("JULIA_ENV" => env.project_file) do
        LOAD_PATH = filter(x -> x isa AbstractString, Base.LOAD_PATH)
        for (uuid, name, hash, build_file) in builds
            log_file = splitext(build_file)[1] * ".log"
            Base.info("Building $name [$(string(hash)[1:16])]...")
            Base.info(" → $log_file")
            code = """
                empty!(Base.LOAD_PATH)
                append!(Base.LOAD_PATH, $(repr(LOAD_PATH)))
                empty!(Base.LOAD_CACHE_PATH)
                append!(Base.LOAD_CACHE_PATH, $(repr(Base.LOAD_CACHE_PATH)))
                empty!(Base.DL_LOAD_PATH)
                append!(Base.DL_LOAD_PATH, $(repr(Base.DL_LOAD_PATH)))
                cd($(repr(dirname(build_file))))
                include($(repr(build_file)))
                """
            cmd = ```
                $(Base.julia_cmd()) -O0 --color=no --history-file=no
                --startup-file=$(Base.JLOptions().startupfile != 2 ? "yes" : "no")
                --compilecache=$((VERSION < v"0.7.0-DEV.1735" ? Base.JLOptions().use_compilecache :
                                  Base.JLOptions().use_compiled_modules) != 0 ? "yes" : "no")
                --eval $code
                ```
            open(log_file, "w") do log
                success(pipeline(cmd, stdout=log, stderr=log))
            end ? Base.rm(log_file, force=true) :
            warn("Error building `$name`; see log file for further info")
        end
    end
    return
end

function rm(env::EnvCache, pkgs::Vector{PackageSpec})
    drop = UUID[]
    # find manifest-mode drops
    for pkg in pkgs
        pkg.mode == :manifest || continue
        info = manifest_info(env, pkg.uuid)
        if info != nothing
            pkg.uuid in drop || push!(drop, pkg.uuid)
        else
            str = has_name(pkg) ? pkg.name : string(pkg.uuid)
            warn("`$str` not in manifest, ignoring")
        end
    end
    # drop reverse dependencies
    while !isempty(drop)
        clean = true
        for (name, infos) in env.manifest, info in infos
            haskey(info, "uuid") && haskey(info, "deps") || continue
            deps = map(UUID, values(info["deps"]))
            isempty(drop ∩ deps) && continue
            uuid = UUID(info["uuid"])
            uuid ∉ drop || continue
            push!(drop, uuid)
            clean = false
        end
        clean && break
    end
    # find project-mode drops
    for pkg in pkgs
        pkg.mode == :project || continue
        found = false
        for (name::String, uuid::UUID) in env.project["deps"]
            has_name(pkg) && pkg.name == name ||
            has_uuid(pkg) && pkg.uuid == uuid || continue
            !has_name(pkg) || pkg.name == name ||
                error("project file name mismatch for `$uuid`: $(pkg.name) ≠ $name")
            !has_uuid(pkg) || pkg.uuid == uuid ||
                error("project file UUID mismatch for `$name`: $(pkg.uuid) ≠ $uuid")
            uuid in drop || push!(drop, uuid)
            found = true
            break
        end
        found && continue
        str = has_name(pkg) ? pkg.name : string(pkg.uuid)
        warn("`$str` not in project, ignoring")
    end
    # delete drops from project
    n = length(env.project["deps"])
    if VERSION < v"0.7.0-DEV.1393"
        filter!(env.project["deps"]) do _, uuid
            UUID(uuid) ∉ drop
        end
    else
        filter!(env.project["deps"]) do _uuid # (_, uuid) doesn't parse on 0.6
            _, uuid = _uuid
            UUID(uuid) ∉ drop
        end
    end
    if length(env.project["deps"]) == n
        info("No changes")
        return
    end
    # only keep reachable manifest entires
    prune_manifest(env)
    # update project & manifest
    write_env(env)
end

function add(env::EnvCache, pkgs::Vector{PackageSpec})
    # copy added name/UUIDs into project
    for pkg in pkgs
        env.project["deps"][pkg.name] = string(pkg.uuid)
    end
    # if a package is in the project file and
    # the manifest version in the specified version set
    # then leave the package as is at the installed version
    for (name::String, uuid::UUID) in env.project["deps"]
        info = manifest_info(env, uuid)
        info != nothing && haskey(info, "version") || continue
        version = VersionNumber(info["version"])
        for pkg in pkgs
            pkg.uuid == uuid && version ∈ pkg.version || continue
            pkg.version = version
        end
    end
    # resolve & apply package versions
    resolve_versions!(env, pkgs)
    new = apply_versions(env, pkgs)
    write_env(env) # write env before building
    build_versions(env, new)
end

function up(env::EnvCache, pkgs::Vector{PackageSpec})
    # resolve upgrade levels to version specs
    for pkg in pkgs
        pkg.version isa UpgradeLevel || continue
        level = pkg.version
        info = manifest_info(env, pkg.uuid)
        ver = VersionNumber(info["version"])
        if level == UpgradeLevel(:fixed)
            pkg.version = VersionNumber(info["version"])
        else
            r = level == UpgradeLevel(:patch) ? VersionRange(ver.major, ver.minor) :
                level == UpgradeLevel(:minor) ? VersionRange(ver.major) :
                level == UpgradeLevel(:major) ? VersionRange() :
                    error("unexpected upgrade level: $level")
            pkg.version = VersionSpec(r)
        end
    end
    # resolve & apply package versions
    resolve_versions!(env, pkgs)
    new = apply_versions(env, pkgs)
    write_env(env) # write env before building
    build_versions(env, new)
end

function test(env::EnvCache, pkgs::Vector{PackageSpec}; coverage=false)
    # See if we can find the test files for all packages
    missing_runtests = String[]
    testfiles        = String[]
    version_paths    = String[]
    for pkg in pkgs
        info = manifest_info(env, pkg.uuid)
        haskey(info, "hash-sha1") || cmderror("Could not find hash-sha for package $(pkg.name)")
        version_path = find_installed(pkg.uuid, SHA1(info["hash-sha1"]))
        testfile = joinpath(version_path, "test", "runtests.jl")
        if !isfile(testfile)
            push!(missing_runtests, pkg.name)
        end
        push!(version_paths, version_path)
        push!(testfiles, testfile)
    end
    if !isempty(missing_runtests)
        cmderror(length(missing_runtests) == 1 ? "Package " : "Packages ",
                join(missing_runtests, ", "),
                " did not provide a `test/runtests.jl` file")
    end

    pkgs_errored = []
    for (pkg, testfile, version_path) in zip(pkgs, testfiles, version_paths)
        info("Testing $(pkg.name) located at $version_path")
        if env.preview[]
            info("In preview mode, skipping tests for $(pkg.name)")
            continue
        end
        # TODO, cd to test folder (need to be careful with getting the same EnvCache
        # as for this session in that case
        compilemod_opt, compilemod_val = VERSION < v"0.7.0-DEV.1735" ?
                ("compilecache" ,     Base.JLOptions().use_compilecache) :
                ("compiled-modules",  Base.JLOptions().use_compiled_modules)

        testcmd = `"import Pkg3; include(\"$testfile\")"`
        cmd = ```
            $(Base.julia_cmd())
            --code-coverage=$(coverage ? "user" : "none")
            --color=$(Base.have_color ? "yes" : "no")
            --$compilemod_opt=$(Bool(compilemod_val) ? "yes" : "no")
            --check-bounds=yes
            --startup-file=$(Base.JLOptions().startupfile != 2 ? "yes" : "no")
            $testfile
        ```
        try
            run(cmd)
            info("$(pkg.name) tests passed")
        catch err
            push!(pkgs_errored, pkg.name)
        end

    end

    if !isempty(pkgs_errored)
        cmderror(length(pkgs_errored) == 1 ? "Package " : "Packages ",
                 join(pkgs_errored, ", "),
                 " errored during testing")
    end
end
end # module


text:
using FEMUCH
using Test

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


text:
## log density function of N(mu, sd^2)
function logdnorm(y::Array, mu::Float64, sd::Float64)
    return( -0.5 * log(2*pi) .- 1.0 * log(sd) .- 1/(2*sd^2) * (y .- mu).^2 )
end

function logdnorm(y::Float64, mu::Float64, sd::Float64)
    return( -0.5 * log(2*pi) .- 1.0 * log(sd) .- 1/(2*sd^2) * (y .- mu).^2 )
end
function toy_mh(y::Array, S::Int64)
    s2 = 1.0; t2 = 10.0; mu = 5.0;

    theta = 0.0; delta2 = 2; S = 10000; THETA::Array{Float64} = []

    for s = 1:S
        theta_star = randn() * sqrt(delta2) + theta
        log_r = ( sum(logdnorm(y, theta_star, sqrt(s2))) + logdnorm(theta_star, mu, sqrt(t2)) ) -
                     ( sum(logdnorm(y, theta, sqrt(s2))) + logdnorm(theta, mu, sqrt(t2)) )

        if log(rand()) < log_r
            theta = theta_star
        end
        THETA = append!(THETA, theta)
    end
    return(THETA)
end

#y = [9.37, 10.18, 9.16, 11.60, 10.33]
y = randn(100) * sqrt(10)
res = toy_mh(y, Int64(1e10))
# using PyPlot
# plt[:hist](toy_mh(5000))
# show()

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
module Breadth
import Base.~, Base.∈, Base.∪, Base.∩, Base.-#, Base.:(≡)
using Memoize

#include("tree_type3/__init__.jl")
include("tree_type/TreeType.jl")

# macro ex(func)
#     return quote
#         $(esc(func))(x::Vararg) = sum(x)
#         export $(esc(func))
#     end
# end

# @ex ∨

end


text:
using CSV, DataFrames, Distributed, Dates, LinearAlgebra, Jevo, Distributions, DelimitedFiles, LambertW
using SharedArrays, TimerOutputs
# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

# Get number of workers
if length(ARGS) == 1
    addprocs(parse(Int64, ARGS[1]))
elseif length(ARGS) > 1
    throw(ArgumentError("Only one command line argument (cores)."))
end

# Import packages needed for all workers
@everywhere  begin
    using Jevo
    using Distributions
    using DelimitedFiles
    using SharedArrays
    using LambertW
    using LinearAlgebra
end

# Parameters
@everywhere begin
    reps = 120
    steps = 2 * 10^7
    rho = [0, 0.1, 0.5, 1, 2, 4]
    l_0 = 10
    N = 100
    nu = 1/N

    f0 = 50/2N
    ϵ = 2
    n = 20
    emat = ϵ * (ones(n, n) - Matrix{Float64}(I, n, n))

    l0_kappa(kappa, l) = 1/2 * lambertw(2 * ϵ^2 * N * l * f0 * (n-1)/n^2 * exp(10)/(1+kappa))
    #fl(l_opt) = l0_kappa(0, 10)/l_opt^2 * n^2 / (n-1) * 1/ϵ
    fl(l_opt) = l0_kappa(0, 10)^2 / l_opt^2 * n^2 / (n-1) * 1/2

end

E_results = SharedArray{Float64, 2}(length(rho), reps)
l_results = SharedArray{Float64, 2}(length(rho), reps)
rho_list = SharedArray{Float64, 2}(length(rho), reps)

# Function to run one simulation
@everywhere function run(N, rho, nu, l_0, n, emat, steps)
    # Initiate population
    pop = Jevo.mono_pop(N=100, l=l_0, n=n, m=n)
    Jevo.initiate!(pop, opt=true)
    F = Jevo.fermi_fitness(f0=f0, fl=fl(10)/2N)
    for i in 1:steps
        Jevo.bp_substitution!(pop, emat, F)
        if rand() < rho/N
            Jevo.driver_mutation!(pop)
        end
        if rand() < nu
            Jevo.l_substitution!(pop, emat, F)
        end
        # Recover lost sites
        if length(pop.seqs) < 6
            Jevo.initiate!(pop, opt=true)
        end
    end
    return Jevo.get_energy(pop, emat), length(pop.seqs)
end

# Write Metadata
name = "supp_n_20"
open(date * name * "_METADATA.txt", "a") do io
    write(io, "N=$N\n")
    write(io, "f0=$f0\n")
    write(io, "fl=$fl\n")
    write(io, "repetitions=$reps\n")
    write(io, "steps=$steps\n")
    write(io, "rho=$rho\n")
    write(io, "nu=$nu\n")
    write(io, "l_0=$l_0\n")
end

# Run simulations and enjoy speed
@sync @distributed for j in 1:reps
    for r in 1:length(rho)
        E, L = run(N, rho[r], nu, l_0, n, emat, steps)
        E_results[r, j] = E
        l_results[r, j] = L
        rho_list[r, j] = rho[r]
    end
    println("Run $j done.")
end

df = DataFrame(gamma=[(E_results...)...], l=[(l_results...)...], rho=[(rho_list...)...])
CSV.write(name * ".csv", df)


text:
module MRIReco

using Reexport
using FileIO
using ProgressMeter
@reexport using ImageUtils
using NFFT
using Distributions
using LinearOperators
@reexport using RegularizedLeastSquares
using StatsBase
using LinearAlgebra
using Random
@reexport using FFTW
using Distributed
using Graphics: @mustimplement
using ColorTypes
using ColorVectorSpace
using Wavelets
using HDF5
using LightXML
using NIfTI
@reexport using Unitful

include("Trajectories/Trajectories.jl")
include("Sequences/Sequence.jl")
include("Datatypes/Datatypes.jl")
include("Tools/Tools.jl")
include("Operators/Operators.jl")
include("Simulation/Simulation.jl")
include("Reconstruction/Reconstruction.jl")
include("IO/IO.jl")
include("Sampling/Sampling.jl")

end # module


text:
using Test, PartialSvdStoch, Serialization, Distributions

using Logging
using Base.CoreLogging

logger = ConsoleLogger(stdout, CoreLogging.Debug)
global_logger(logger)

include("testvrpca.jl")
include("testifsvd.jl")

# a data file for a matrix(11520, 22499) can be found in ./data/plimat.serialized
# mass spectrometry image

plimat = Some(Matrix{Float64})
dataf = "data/plimat.serialized"
if isfile(dataf)
    plimat = Some(Serialization.deserialize(open(dataf, "r")))
end

# This test compares singular values from exact svd with those extracted from
# incremental svd from LowRankApproxMc

@testset "ifsvd_versus_svd" begin
    @info "========================================================="
    @info "\n\n in test ifsvd_versus_svd"
    d = 500
    n = 5000
    k = 50
    mat = rand(d,n)
    for j in 1:k
    mat[:,j] = mat[:,j] * 20.
    end
    for j in 1:5
        mat[:,j+k] = mat[:,j+k] * 100.
    end
    @test testifastsvd_S(mat, 200)
end

# This test shows how in presence high variability of data
# it is possible to get a competitive approximation with incremental svd.

@testset "ifsvd_lowRankApproxEpsil" begin
    @info "========================================================="
    @info "\n\n in test ifsvd_lowRankApproxEpsil"
    d = 10000
    n = 20000
    mat = zeros(d,n)
    #
    Xlaw = Normal(0, 1000.)
    # multiply some columns
    blocsize = 100
    for bloc in 1:round(Int,n/blocsize)
        val = abs(rand(Xlaw))
        for i in 1:d
            if rand() < 0.01
                for j in 1:blocsize
                    mat[i,(bloc-1)*blocsize+j] = val * (1. - rand()/10.)
                end
            end
        end
    end
    @test testifastsvd_epsil(mat, 100, 0.05)
end

#
# test vrpca_epsil

# This tests asks for an approx at 0.01
# It extracts a rank 100 matrix at 0.005 precision within 3s (after first compilation).
# With LowRankApprox.prange alone we need to go up to rank=500 (initialized with
# opts = LowRankApprox.LRAOptions(rank=500, rtol = 0.0001) to achieve a relative
# Froebonius error of  0.005. It then runs in 1.15s.

# So we see that LowRankApprox is really fast but needs many more singular vectors
# and that VrPCA gives a clear reduction of relative error.
# Of course VrPCA initialized with LowRankApprox would give an analogous reduction
# in relative error. Cf testset "vrpca_withLoxwRankApprox"
#

@testset "vrpca_epsil" begin
    @info "========================================================="
    @info "\n\n in test vrpca_epsil"
    d = 5000
    n = 10000
    k = 100
    maxval = 1000
    mat = rand(d,n)
    # multiply some columns
    for i in 1:k
        mat[:,i] = mat[:,i] * maxval
    end
    @test testvrpca_withepsil(mat, 0.01)
end

@testset "vrpca_withLowRankApprox" begin
    @info "========================================================="
    @info "\n\n in test vrpca_withLowRankApprox"
    d = 5000
    n = 10000
    k = 100
    maxval = 1000
    mat = rand(d,n)
    # multiply some columns. but do not sample them( some will be sampled twice and will affect sketching)
    for i in 1:k
        mat[:,i] = mat[:,i] * maxval
    end
    @test testvrpca_withLowRankApprox(mat, 100)
end


text:
# # Adding more states

using Revise #src
using Test #jl
using Aiyagari
using QuantEcon, Parameters
using Plots #md

PKG_HOME = joinpath(dirname(pathof(Aiyagari)), "..")

# ## Enlarging the state space

# First, let's add an additional exogenous state variable

z_grid = [0.5; 1.0; 1.5]
z_prob = [0.7 0.15 0.15;
          0.2 0.6 0.2;
          0.15 0.15 0.7]
z_MC = MarkovChain(z_prob, z_grid, :z)

x_grid = [1, 2, 3]
x_prob = [0.7 0.3 0.0;
          0.0 0.99 0.01;
          1.0 0.0 0.0]

x_MC = MarkovChain(x_prob, x_grid, :x)

exo1 = ExogenousStateSpace([z_MC])
exo2 = ExogenousStateSpace([z_MC, x_MC])

# Second, let's also add an additional endogenous state variable

w_grid = LinRange(0.0, 3.0, 20)
h_grid = LinRange(0.0, 1.0, 15)

endo1 = EndogenousStateSpace((w=w_grid, ))
endo2 = EndogenousStateSpace((w=w_grid, h=h_grid))

# ## Define some other things

include(joinpath(PKG_HOME, "test", "housing-simple-nlopt.jl"))

u(c; γ) = c^(1-γ) / (1-γ)

function u(c,h; ξ=0.8159, ρ=map(s -> (s-1)/s, 0.13), γ=2.0)
  C = (ξ * h^ρ + (1-ξ) * c^ρ)^(1/ρ)
  u(C, γ=γ)
end

mutable struct HousingAS{T1,T2,T3,T4} <: AggregateState
  r::T1
  p::T2 # house price
  ρ::T3 # rent
  dist::T4 # the distribution over idiosynchratic states
end

function HousingAS(r, p, endo, exo, param; ρ=p * (param.δ + r))
  dist_proto = zeros((length(endo), length(exo)))
  HousingAS(r, p, ρ, dist_proto)
end

# ## Solve the model with and without redundant states

r = 0.15
p = 0.9
param  = (β = 0.7, θ = 0.9, δ = 0.1, h_thres = eps())
agg_state11 = HousingAS(r, p, endo1, exo1, param)
agg_state21 = HousingAS(r, p, endo2, exo1, param)
agg_state12 = HousingAS(r, p, endo1, exo2, param)
agg_state22 = HousingAS(r, p, endo2, exo2, param)

@time out11 = solve_bellman(endo1, exo1, agg_state11, param, Owner(Aiyagari.Unconditional(), Aiyagari.NoState()), rtol=√eps())
#0.43 s (49 it)
@time out21 = solve_bellman(endo2, exo1, agg_state21, param, Owner(Aiyagari.Unconditional(), Aiyagari.IsState()), rtol=√eps())
#4.85 s (49 it)
@time out12 = solve_bellman(endo1, exo2, agg_state12, param, Owner(Aiyagari.Unconditional(), Aiyagari.NoState()), rtol=√eps())
#1.03 s (49 it)
@time out22 = solve_bellman(endo2, exo2, agg_state22, param, Owner(Aiyagari.Unconditional(), Aiyagari.IsState()), rtol=√eps())
#14.41 s (49 it)

# Check if results are the same
val11 = out11.val
plt11 = plot(w_grid, val11, alpha=0.35, title="no redundant states") #md

plt21 = plot!(plt11, legend=false, title="redundant endogenous state") #md
@testset "irrelevant endogenous state" begin #jl
  val21 = reshape(out21.val, (size(endo2)..., length(exo1)))
  let #jl
    max_diff = Inf #jl
    for i_h in 1:length(h_grid)
      val_i = val21[:,i_h,:]
      plot!(plt21, w_grid, val_i, alpha = 0.35 / length(h_grid)) #md
      @test all(val11 .≈ val_i) #jl
      max_i = maximum(abs, val11 .- val_i) #jl
      max_diff = max_diff < max_i ? max_diff : max_i  #jl
     end
    println(max_diff) #jl
  end #jl
end #jl

@testset "irrelevant exogenous state" begin #jl
  val12 = reshape(out12.val, (length(endo1), size(exo2)...)) #jl
  let #jl
    max_diff = Inf #jl
    for i_x in 1:length(x_grid) #jl
      val_i = val12[:,:,i_x] #jl
      @test all(val11 .≈ val_i) #jl
      max_i = maximum(abs, val11 .- val_i) #jl
      max_diff = max_diff < max_i ? max_diff : max_i  #jl
    end #jl
    println(max_diff) #jl
  end #jl
end #jl
plt12 = plot!(plt11, w_grid, out12.val, alpha=0.35, legend=false,  title="redundant exogenous state") #md

plt22 = plot!(plt11, legend=false, title="redundant endo + exo state")  #md
@testset "irrelevant endogenous and exogenous states" begin #jl
   val22 = reshape(out22.val, (size(endo2)..., size(exo2)...))
   let #jl
     max_diff = Inf #jl
     for i_x in 1:length(x_grid)
       for i_h in 1:length(h_grid)
        val_i = val22[:,i_h,:,i_x]
        plot!(plt22, w_grid, val_i, alpha = 0.35 / length(x_grid) / length(h_grid)) #md
        @test all(val11 .≈ val_i) #jl
        max_i = maximum(abs, val11 .- val_i) #jl
        max_diff = max_diff < max_i ? max_diff : max_i  #jl
      end
    end
    println(max_diff) #jl
  end #jl
end #jl

display(plt11) #md
#-
display(plt12) #md
#-
display(plt21) #md
#-
display(plt22) #md
#plot(plt11, plt12, plt21, plt22) #md

# TODO: simple (linear or binary adjustment costs)
# TODO: house quality shock


text:
alpha = 0.5
accuracy = 0.0000001
edgelist = "test/fixtures/simple.adj"

rank = Rank{Float64,Int64}(alpha, accuracy, edgelist)
r = stationary_distribution(rank)
count = length(r)

expected = [float(5 / 18), float(4 / 9), float(5 / 18)]
precision = 5

@assert isequal(floor(expected, precision), floor(r, precision))


text:
export createDirectory

# has to be split from the declaration of Directory due to dependancy on files that require Directory

function getDirectoryEntries(directory::Directory{GID, PID, LID}, map::BlockMap{GID, PID, LID},
        globalEntries::AbstractArray{Number}, high_rank_sharing_procs::Bool=false)::Tuple{AbstractArray{PID}, AbstractArray{LID}} where GID <: Integer where PID <: Integer where LID <: Integer
    getDirectoryEntries(directory, map, Array{GID, 1}(undef, globalEntries), high_rank_sharing_procs)
end

function getDirectoryEntries(directory::Directory{GID, PID, LID}, map::BlockMap{GID, PID, LID},
        globalEntries::AbstractArray{GID})::Tuple{AbstractArray{PID}, AbstractArray{LID}} where GID <: Integer where PID <: Integer where LID <: Integer
    getDirectoryEntries(directory, map, globalEntries, false)
end

"""
    createDirectory(comm::Comm, map::BlockMap)
Create a directory object for the given Map
"""
function createDirectory(comm::Comm{GID, PID, LID}, map::BlockMap{GID, PID, LID})::BasicDirectory{GID, PID, LID} where GID <: Integer where PID <: Integer where LID <: Integer
    BasicDirectory{GID, PID, LID}(map)
end


text:
# Translation utilities

export PhaseShiftLinOp, JacobianPhaseShiftEvaluated, PhaseShiftParametericLinOp, PhaseShiftParametericDelayedEval
export phase_shift_linop, phase_shift, Jacobian, ∂

## Linear operator

struct PhaseShiftLinOp{T}<:AbstractLinearOperator{Complex{T},2,2}
    phase_shift::AbstractArray{Complex{T},2}
end

function phase_shift_linop(K::AbstractKSpaceSampling{T}, τ::AbstractArray{T,2}) where {T<:Real}
    phase_shift  = exp.(-im*sum(coord(K).*reshape(τ,:,1,3); dims=3)[:,:,1])
    return PhaseShiftLinOp{T}(phase_shift)
end

AbstractLinearOperators.domain_size(P::PhaseShiftLinOp) = size(P.phase_shift)
AbstractLinearOperators.range_size(P::PhaseShiftLinOp) = domain_size(P)

AbstractLinearOperators.matvecprod(P::PhaseShiftLinOp{T}, d::AbstractArray{Complex{T},2}) where {T<:Real} = P.phase_shift.*d
Base.:*(P::PhaseShiftLinOp{T}, d::AbstractArray{Complex{T},3}) where {T<:Real} = P.phase_shift.*d
AbstractLinearOperators.matvecprod_adj(P::PhaseShiftLinOp{T}, d::AbstractArray{Complex{T},2}) where {T<:Real} = conj(P.phase_shift).*d

## Parameteric linear operator

struct PhaseShiftParametericLinOp{T}
    K::AbstractKSpaceSampling{T}
end

phase_shift(K::AbstractKSpaceSampling{T}) where {T<:Real} = PhaseShiftParametericLinOp{T}(K)

(P::PhaseShiftParametericLinOp{T})(τ::AbstractArray{T,2}) where {T<:Real} = phase_shift_linop(P.K, τ)

(P::PhaseShiftParametericLinOp)() = P

## Parameteric delayed evaluation

struct PhaseShiftParametericDelayedEval{T}
    P::PhaseShiftParametericLinOp{T}
    d::AbstractArray{Complex{T},2}
end

Base.:*(P::PhaseShiftParametericLinOp{T}, d::AbstractArray{Complex{T},2}) where {T<:Real} =  PhaseShiftParametericDelayedEval{T}(P, d)

(Pd::PhaseShiftParametericDelayedEval{T})(τ::AbstractArray{T,2}) where {T<:Real} = phase_shift_linop(τ, Pd.P.K)*Pd.d

## Jacobian of phase-shift evaluated

struct JacobianPhaseShiftEvaluated{T}<:AbstractLinearOperator{Complex{T},2,2}
    J::AbstractArray{Complex{T},3}
end

function Jacobian(Pd::PhaseShiftParametericDelayedEval{T}, τ::AbstractArray{T,2}) where {T<:Real}
    Pτ = phase_shift_linop(Pd.P.K, τ)
    Pτd = Pτ*Pd.d
    J = -im*coord(Pd.P.K).*reshape(Pτd, size(Pτd)..., 1)
    return Pτd, Pτ, JacobianPhaseShiftEvaluated{T}(J)
end

∂(Pd::PhaseShiftParametericDelayedEval{T}, τ::AbstractArray{T,2}) where {T<:Real} = Jacobian(Pd, τ)

AbstractLinearOperators.domain_size(∂Pd::JacobianPhaseShiftEvaluated) = (size(∂Pd.J,1),3)
AbstractLinearOperators.range_size(∂Pd::JacobianPhaseShiftEvaluated) = size(∂Pd.J)[1:2]

AbstractLinearOperators.matvecprod(∂Pd::JacobianPhaseShiftEvaluated{T}, Δτ::AbstractArray{CT,2}) where {T<:Real,CT<:RealOrComplex{T}} = sum(∂Pd.J.*reshape(Δτ, :, 1, 3); dims=3)[:,:,1]
Base.:*(∂Pd::JacobianPhaseShiftEvaluated{T}, Δτ::AbstractArray{T,2}) where {T<:Real} = AbstractLinearOperators.matvecprod(∂Pd, Δτ)
AbstractLinearOperators.matvecprod_adj(∂Pd::JacobianPhaseShiftEvaluated{T}, Δd::AbstractArray{Complex{T},2}) where {T<:Real} = sum(conj(∂Pd.J).*reshape(Δd, size(Δd)..., 1); dims=2)[:,1,:]

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
"""
IntegRuleModule
Module for numerical integration rules.
"""
module IntegRuleModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

abstract type IntegRule end

"""
TriRule

Triangular quadrature rule.

Used for integration on the standard triangle.
"""
struct TriRule <: IntegRule
    npts::FInt
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 2}
end

"""
    TriRule(npts=1)

Type for triangular quadrature rule.  Used for integration of the standard
triangle, which is between 0 and 1 in both parametric coordinates.
`npts` = number of points (1-- one-point rule, 3 -- three-point rule,
 6 -- six point rule, 10 -- Strang 10 point, order 13, degree of precision 7, rule).
"""
function TriRule(npts=1)
    if npts == 1 # integrates exactly linear polynomials
        param_coords = [1.0/3. 1.0/3.];
        weights = reshape([1.0]/2.0,1,1);
    elseif npts == 3 # integrates exactly quadratic polynomials
        param_coords = [ 2.0/3 1.0/6; 1.0/6 2.0/3; 1.0/6 1.0/6 ];
        weights = [1.0/3 1.0/3 1.0/3]/2;
    elseif npts == 6 # integrates exactly quartic polynomials
        param_coords = [ 0.816847572980459 0.091576213509771;
                        0.091576213509771 0.816847572980459;
                        0.091576213509771 0.091576213509771;
                        0.108103018168070 0.445948490915965;
                        0.445948490915965 0.108103018168070;
                        0.445948490915965 0.445948490915965];
        weights = [0.109951743655322*[1, 1, 1] 0.223381589678011*[1, 1, 1] ]/2;
    elseif npts == 13 # integrates exactly quartic polynomials
        param_coords = [0.333333333333333  0.333333333333333
                        0.479308067841923  0.260345966079038
                        0.260345966079038  0.479308067841923
                        0.260345966079038  0.260345966079038
                        0.869739794195568  0.065130102902216
                        0.065130102902216  0.869739794195568
                        0.065130102902216  0.065130102902216
                        0.638444188569809  0.312865496004875
                        0.638444188569809  0.048690315425316
                        0.312865496004875  0.638444188569809
                        0.312865496004875  0.048690315425316
                        0.048690315425316  0.638444188569809
                        0.048690315425316  0.312865496004875
        ];
        weights = [ -0.149570044467670
                    0.175615257433204
                    0.175615257433204
                    0.175615257433204
                    0.053347235608839
                    0.053347235608839
                    0.053347235608839
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
        ]'/2;
    else
        #nothing doing: this input is wrong
        error( "Unknown number of integration points" )
    end
    rule  = TriRule(npts, param_coords, weights)
    return  rule
end

"""
    GaussRule

The Gauss rule.

The rule is applicable for a tensor product of  intervals -1 <=x<= +1.
"""
struct GaussRule <: IntegRule
    dim::FInt
    order::FInt
    npts::FInt
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 2}
end

"""
    GaussRule(dim=1, order=1)

Gauss rule.
"""
function GaussRule(dim=1, order=1)
    @assert 1 <= dim <= 3 "Gauss rule of dimension $(dim) not available"
    if     (order==1)
        param_coords = [ 0.0 ];
        weights  = [ 2.0 ];
    elseif (order==2)
        param_coords = [ -0.577350269189626 0.577350269189626 ];
        weights  = [ 1.0 1.0 ];
    elseif (order==3)
        param_coords = [ -0.774596669241483  0.0  0.774596669241483 ];
        weights  = [ 0.5555555555555556 0.8888888888888889 0.5555555555555556 ];
    elseif (order==4)
        param_coords = [ -0.86113631159405  -0.33998104358486   0.33998104358486   0.86113631159405];
        weights  = [ 0.34785484513745   0.65214515486255   0.65214515486255   0.34785484513745];
    else
        error("Gauss rule of order $(order) not available")
    end
    if (dim==1)
        param_coords = transpose(param_coords);
        weights = weights;
        npts=order
        param_coords= reshape(param_coords',npts,1)
        weights= reshape(weights,npts,1)
    elseif (dim==2)
        pc=param_coords;
        w=weights;
        npts=order^2
        param_coords=zeros(FFlt,npts,2);
        weights=zeros(FFlt,npts,1);
        r=1
        for i=1:order
            for j=1:order
                param_coords[r,:]=[pc[i] pc[j]];
                weights[r]=w[i]*w[j];
                r=r+1
            end
        end
    else # (dim==3)
        pc=param_coords;
        w=weights;
        npts=order^3
        param_coords=zeros(FFlt,npts,3);
        weights=zeros(FFlt,npts,1);
        r=1
        for i=1:order
            for j=1:order
                for k=1:order
                    param_coords[r,:]=[pc[i] pc[j] pc[k]];
                    weights[r]=w[i]*w[j]*w[k];
                    r=r+1
                end
            end
        end
    end
    return GaussRule(order,dim,npts,param_coords,weights)
end

"""
    TetRule

Tetrahedral quadrature rule.

Used for integration on the standard tetrahedron.
"""
struct TetRule <: IntegRule
    npts::FInt
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 2}
end

"""
    TetRule(npts=1)

Tetrahedral integration rule.
npts=number of points (1-- one-point rule, 4 -- four-point rule, 5 -- five point
rule).
"""
function TetRule(npts=1)
    if npts == 1 # integrates exactly linear polynomials
        param_coords = reshape([0.25,0.25,0.25],1,3);
        weights = reshape([1.0]/6.0,1,1);
    elseif npts == 4 # integrates exactly quadratic polynomials
        param_coords = [[0.13819660 0.13819660 0.13819660];
        [0.58541020 0.13819660 0.13819660];
        [0.13819660 0.58541020 0.13819660];
        [0.13819660 0.13819660 0.58541020]];;
        weights = [ 0.041666666666666666667   0.041666666666666666667   0.041666666666666666667   0.041666666666666666667];
    elseif npts == 5 #  Zienkiewicz #3.
        a =   1.0 / 6.0;
        b =   0.25;
        c =   0.5;
        d = - 0.8;
        e =   0.45;
        param_coords = [[b b b];
        [c a a];
        [a c a];
        [a a c];
        [a a a]];
        weights = [d  e  e  e  e]/6;
    else
        #nothing doing: this input is wrong
        error( "Unknown number of integration points" )
    end

    return  TetRule(npts, param_coords, reshape(weights,length(weights),1))
end

"""
    PointRule

Class of point quadrature rule.

Used for integration on the standard "point" shape.
"""
struct PointRule <: IntegRule
    npts::FInt
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 2}
end

"""
    PointRule()

POINT integration rule.
"""
function PointRule()
    return PointRule(1, reshape([1.0], 1, 1), reshape([1.0], 1, 1))
end

"""
    SimplexRule

Simplex quadrature rule.

Used for integration on the standard triangle or the standard tetrahedron.
"""
struct SimplexRule <: IntegRule
    npts::FInt
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 2}
end

function SimplexRule(dim=1, npts=1)
    @assert 0 <= dim <= 3 "Simplex rule of dimension $(dim) not available"

    if dim == 0
        r = PointRule()
        return SimplexRule(r.npts, r.param_coords, r.weights)
    elseif dim == 1
        r = GaussRule(dim, npts)
        return SimplexRule(r.npts, r.param_coords, r.weights)
    elseif dim == 2
        r = TriRule(npts)
        return SimplexRule(r.npts, r.param_coords, r.weights)
    else
        r = TetRule(npts)
        return SimplexRule(r.npts, r.param_coords, r.weights)
    end
end

"""
    TrapezoidalRule

The trapezoidal rule.

The rule is applicable for a tensor product of  intervals -1 <=x<= +1.
"""
struct TrapezoidalRule <: IntegRule
    dim::FInt
    npts::FInt
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 2}
end

"""
    TrapezoidalRule(dim=1)

Trapezoidal rule.
"""
function TrapezoidalRule(dim=1)
    @assert 1 <= dim <= 3 "Trapezoidal rule of dimension $(dim) not available"
    order = 2 # order of the trapezoidal rule: fixed
    param_coords = [ -1.0 1.0 ];
    weights  = [ 1.0 1.0 ];
    if (dim==1)
        param_coords = transpose(param_coords);
        weights = weights;
        npts=order^1
        param_coords= reshape(param_coords',npts,1)
        weights= reshape(weights,npts,1)
    elseif (dim==2)
        pc=param_coords;
        w=weights;
        npts=order^2
        param_coords=zeros(FFlt,npts,2);
        weights=zeros(FFlt,npts,1);
        r=1
        for i=1:order
            for j=1:order
                param_coords[r,:]=[pc[i] pc[j]];
                weights[r]=w[i]*w[j];
                r=r+1
            end
        end
    else # (dim==3)
        pc=param_coords;
        w=weights;
        npts=order^3
        param_coords=zeros(FFlt,npts,3);
        weights=zeros(FFlt,npts,1);
        r=1
        for i=1:order
            for j=1:order
                for k=1:order
                    param_coords[r,:]=[pc[i] pc[j] pc[k]];
                    weights[r]=w[i]*w[j]*w[k];
                    r=r+1
                end
            end
        end
    end
    return TrapezoidalRule(dim,npts,param_coords,weights)
end

end


text:
letters = ['a':'z';
		'A':'Z';
		'0':'9';
		' ';
		'	';
		'.';
		'"';
		'!';
		'?';
		'/';
		'(';
		')';
		'=';
		'-';
		'[';
		']';
		'{';
		'}';
		'\\']

function setupFontCharacters()
	glPixelStorei(GL_UNPACK_ALIGNMENT, 1)

	tex = GLuint[0]
	glGenTextures(1, tex)
	glActiveTexture(GL_TEXTURE1)
	glBindTexture(GL_TEXTURE_2D, tex[1])

	fontState.textAtlas = tex[1]

	w = 0
	h = 0
	@inbounds for c in letters
		img, metric = renderface(fontState.face, c)
		w += size(img, 1)
		h = max(h, size(img, 2))
	end
	fontState.atlasWidth = w
	fontState.atlasHeight = h

	# glTexImage2D(GL_TEXTURE_2D, 0, GL_RED, w, h, 0, GL_RED, GL_FLOAT, C_NULL)
	glTexImage2D(GL_TEXTURE_2D, 0, GL_RED, w, h, 0, GL_RED, GL_UNSIGNED_BYTE, C_NULL)
	glGenerateMipmap(GL_TEXTURE_2D)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_NEAREST)
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
	max_aniso = [Float32(0.0)]
	glGetFloatv(GL_TEXTURE_MAX_ANISOTROPY_EXT, max_aniso)
	glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY_EXT, 16)

	x = 0
	@inbounds for c in letters
		img, metric = renderface(fontState.face, c)

		w = size(img, 1)
		h = size(img, 2)
		# w1 = w
		# h1 = h

		# the following uses a modification of the signed distance transform
		# that comes from GLVisualize.jl. it follows some of the ideas from
		# https://www.mapbox.com/blog/text-signed-distance-fields/.
		# what we need to do is turn the glyph image into an array of Boolean
		# values. 'true' where the character is and 'false' where it is not.
		# then, we can run the sdf routine on that and save that as our
		# texture. the font shader then takes care of turning this into a
		# textual element.

		# this pads the image so that we can later downsample it properly
		# restrict_steps = 2
		# halfpad = 2*(2^restrict_steps)
		# w = w + 2*halfpad
		# h = h + 2*halfpad
		# in_or_out = zeros(Bool, w, h)

		# for x = 1:w
		#	for y = 1:h
		#		i = x-halfpad
		#		j = y-halfpad
		#		if checkbounds(Bool, size(img), i, j)
		#			in_or_out[i, j] = (img[i, j] >  0.5)
		#		else
		#			in_or_out[i, j] = false
		#		end
		#	end
		# end

		# sd = sdf(in_or_out)

		# for x = 1:restrict_steps
		#	w1, h1 = Images.restrict_size(w1), Images.restrict_size(h1)
		#	sd = Images.restrict(sd)
		# end

		# maxlen = norm([size(img, 1), size(img, 2)])

		# distfield = zeros(Float32, size(sd, 1), size(sd, 2))
		# for x = 1:size(sd, 1)
		#	for y = 1:size(sd, 2)
		#		distfield[x, y] = clamp(sd[x, y]/maxlen, -1, 1)
		#	end
		# end

		# glTexSubImage2D(GL_TEXTURE_2D, 0, x, 0, w, h, GL_RED, GL_FLOAT, distfield)
		glTexSubImage2D(GL_TEXTURE_2D, 0, x, 0, w, h, GL_RED, GL_UNSIGNED_BYTE, img)

        character = textCharacter([w, h],
					[metric.horizontal_bearing[1], metric.horizontal_bearing[2]],
					[metric.advance[1], metric.advance[2]],
					x / fontState.atlasWidth)
        fontState.characters[c] = character
        x += w
	end

    glPixelStorei(GL_UNPACK_ALIGNMENT, 0)
    FT_Done_Face(fontState.face[1])
end


text:
#=
BALL:
- Julia version:
- Author: Dan
- Date: 2021-08-19
=#
module BALL
using Reexport
        #require bijections, DataFrames
    module COMMON
        include("./COMMON/global.jl")
        include("./COMMON/common.jl")
    end

    module CONCEPT
        using ..COMMON
        include("./CONCEPT/selectable.jl")
        include("./CONCEPT/composite_interface.jl")
        include("./CONCEPT/timeStamp.jl")
        include("./CONCEPT/composite.jl")
        include("./CONCEPT/kernel_interfaces.jl")
    end

    module KERNEL
        using ..COMMON
        using ..CONCEPT
        using Bijections
        using StaticArrays
        using DataFrames
        include("./KERNEL/bond.jl")
        include("./KERNEL/PTE.jl")
        include("./KERNEL/atom.jl")
        include("./KERNEL/atom_bijection.jl")
        include("./KERNEL/chain.jl")
        include("./KERNEL/residue.jl")
        include("./KERNEL/system.jl")
        include("./KERNEL/composite_iterator.jl")
        include("./KERNEL/dataframesystem.jl")
        include("./KERNEL/getter_setter.jl")
        include("./KERNEL/aux_functions.jl")
    end

    module SPATIAL
        using ..COMMON
        using ..CONCEPT
        using ..KERNEL
        using StaticArrays
        using DataFrames
        using LinearAlgebra
        include("./SPATIAL/spatial.jl")
    end

    module MOLMEC
        using ..COMMON
        using ..CONCEPT
        using ..KERNEL
        include("./MOLMEC/MMFF94Parameters.jl")
    end

    module FILEFORMATS
        using ..COMMON
        using ..CONCEPT
        using ..KERNEL
        #include("./BioStructures_interface")#slipt dataformat up
        include("./FILEFORMATS/PDBParser.jl")#split dataformat up
    end

    module STRUCTURE
        using DataFrames
        using ..COMMON
        using ..CONCEPT
        using ..KERNEL
        include("./STRUCTURE/simple_molecular_graph.jl")
        include("./STRUCTURE/minimum_cycle_basis.jl")

    end

    module QSAR
        using ..COMMON
        using ..CONCEPT
        using ..KERNEL
        using ..STRUCTURE
        using DataFrames
        include("./QSAR/add_hydrogen_processor.jl")
        include("./QSAR/ring_perception.jl")
    end

    @reexport using .COMMON
    @reexport using .CONCEPT
    @reexport using .KERNEL
    @reexport using .MOLMEC
    @reexport using .FILEFORMATS
    @reexport using .STRUCTURE
    @reexport using .QSAR
    @reexport using .SPATIAL
    export COMMON, CONCEPT, KERNEL, MOLMEC, FILEFORMATS, STRUCTURE, QSAR, SPATIAL

end