text:
"""
    QuickHeaps.AbstractNode{K,V}

is the super-type of nodes with a key of type `K` and a value of type `V`.
Nodes can be used in binary heaps and priority queues where they are used to
store key-value pairs and, by specializing the `Base.lt` method, to implement
how the nodes are ordered.

A node `nd` can be iterated and converted into a pair or a 2-tuple and
conversely:

    k, v = nd               # k is the key and v is the value of the node nd
    Pair(nd)                # yields k=>v
    Tuple(nd)               # yields (k,v)
    QuickHeaps.Node((k, v)) # is the same as QuickHeaps.Node(k, v)
    QuickHeaps.Node(k => v) # is the same as QuickHeaps.Node(k, v)

The `getkey` and (unexported) [`QuickHeaps.getval`](@ref) methods can be used
to respectively retrieve the key and value of a node.  These two methods may be
specialized for the sub-types of [`QuickHeaps.AbstractNode`](@ref).  For
example, a *key-only* node type can be fully implemented by:

    struct KeyOnlyNode{K} <: QuickHeaps.AbstractNode{K,Nothing}
        key::K
    end
    QuickHeaps.getkey(x::KeyOnlyNode) = getfield(x, :key)
    QuickHeaps.getval(x::KeyOnlyNode) = nothing
    KeyOnlyNode(x::Tuple{K,Nothing}) where {K} = KeyOnlyNode{K}(x[1])
    KeyOnlyNode(x::Pair{K,Nothing}) where {K} = KeyOnlyNode{K}(x.first)

"""
abstract type AbstractNode{K,V} end

"""
    QuickHeaps.AbstractPriorityQueue{K,V,T,O}

is the super type of priority queues with nodes consisting in pairs of keys of
type `K` and priority values of type `V`.  Priority queues implement an API
similar to dictionaries with the additional feature of maintaining an ordered
structure so that getting the node of highest priority costs `O(1)` while
pushing a node costs `O(log(n))` with `n` the size of the queue.

Type parameter `T<:AbstractNode{K,V}` is the type of the nodes stored in the
queue and type parameter `O` is the type of the ordering of the queue.  How are
ordered the nodes is completely customizable by specializing the `Base.lt`
method with the following signature:

    lt(o::CustomOrderingType, a::CustomNodeType, b::CustomNodeType)

which shall yield whether node `a` has (strictly) higher priority than node `b`
in the queue and where `CustomOrderingType<:Base.Ordering` and
`CustomNodeType<:QuickHeaps.AbstractNode` are the respective types of the
ordering and of the node of the priority queue.

For the default node type, `QuickHeaps.Node{K,V}`, the implementation is:

    lt(o::Ordering, a::T, b::T) where {T<:Node} = lt(o, a.val, b.val)

In other words, nodes are sorted by their value according to ordering `o`.

"""
abstract type AbstractPriorityQueue{
    K,V,T<:AbstractNode{K,V},O<:Ordering} <: AbstractDict{K,V} end

default_ordering(::Type{<:AbstractPriorityQueue}) = Forward

typename(::Type{<:AbstractPriorityQueue}) = "priority queue"

"""
    QuickHeaps.Node{K,V}(k,v)

yields a node storing key `k` and value `v`.  Type parameters `K` and `V` are
the respective types of the key and of the value.  If omitted the defaults are
`K=typeof(k)` and `V=typeof(v)`.

""" Node

# The following structure was introduced to be more specific than Pair{K,V}.
# So Base.lt can be specialized without type-piracy.
struct Node{K,V} <: AbstractNode{K,V}
    key::K
    val::V
end

# These methods can be specialized.
import Base: getkey

"""
    getkey(x::QuickHeaps.AbstractNode) -> k

yields the key `k` of node `x`.  This method may be specialized for any
sub-types of [`QuickHeaps.AbstractNode`](@ref).

Also see [`QuickHeaps.getval`](@ref).

"""
getkey(x::Node) = getfield(x, :key)

"""
    QuickHeaps.getval(x::QuickHeaps.AbstractNode) -> v

yields the value `v` of node `x`.  This method may be specialized for any
sub-types of [`QuickHeaps.AbstractNode`](@ref).

Also see [`getkey(::QuickHeaps.AbstractNode)`](@ref).

"""
getval(x::Node) = getfield(x, :val)

for type in (:AbstractNode, :Node)
    @eval begin
        $type(x::$type) = x
        $type{K}(x::$type{K}) where {K} = x
        $type{K,V}(x::$type{K,V}) where {K,V} = x
    end
end
Node(x::AbstractNode) = Node(getkey(x), keyval(x))
Node{K}(x::AbstractNode) where {K} = Node{K}(getkey(x), keyval(x))
Node{K,V}(x::AbstractNode{K,V}) where {K,V} = Node{K,V}(getkey(x), keyval(x))

Node(x::Tuple{Any,Any}) = Node(x[1], x[2])
Tuple(x::AbstractNode) = (getkey(x), getval(x))

Node(x::Pair) = Node(x.first, x.second)
Pair(x::AbstractNode) = getkey(x) => getval(x)

iterate(x::AbstractNode) = (getkey(x), first)
iterate(x::AbstractNode, ::typeof(first)) = (getval(x), last)
iterate(x::AbstractNode, ::typeof(last)) = nothing

# Nodes are sorted according to their values.
for O in (:Ordering, :ForwardOrdering, :ReverseOrdering, :FastForwardOrdering)
    @eval begin
        lt(o::$O, a::T, b::T) where {T<:AbstractNode} =
            lt(o, getval(a), getval(b))
    end
end

"""
    PriorityQueue{K,V}([o=Forward,] T=Node{K,V})

yields a priority queue with ordering specified by `o::Ordering` and for nodes
of type `T<:AbstractNode{K,V}` with `K` the type of the keys and `V` the type
of the values encoding the priority.

Type parameters `K` and `V` may be omitted if the node type `T` is specified.

Parameters:

- `K` is key type;
- `V` is value type;
- `T<:AbstractNode{K,V}` is node type;
- `O<:Ordering` is ordering type;
- `I<:Index{K}` is the type of the index;

Keys can be:

- A linear index or a Cartesian index in another array.  The index of the
  priority queue is a, possibly, multi-dimensional array of integers with the
  same dimensions as the other array.  The stored values are the linear indices
  in the vector of nodes of the priority queue or 0 if there no such node.
  This kind of index must keep a track of the number of nodes.

- Anything else: a dictionary will be used to associate keys with indices in
  the vector of nodes of the priority queue.

To enqueue key `k` with value `v` in priority queue `pq`, all the following
are equivalent:

    pq[k] = v
    push!(pq, k => v)
    enqueue!(pq, k, v)
    enqueue!(pq, k => v)

Note that key `k` may already exists in `pq`; in that case, the value
associated with the key is updated and the queue reorderd as needed.  This is
faster than (FIXME: check this) first deleting the key and then enqueuing the
key with the new value.

To extract the key `k` and priority `v` of the node of highest priority from
the queue `pq`, call one of:

    k = dequeue!(pq)
    k, v = pop!(pq)

Method `dequeue!` may be called with a type argument `T` as `dequeue!(T,pq)` to
extract the node of highest priority and apply constructor `T` to convert it to
something else.  For example, `pop!(pq)` is implemented as `dequeue!(Pair,pq)`
and you may call `dequeue!(QuickHeaps.AbstractNode,pq)` to extract the node.
If you want to apply a function `f` taht is not a type constructor to the node
of highest priority, call:

    f(dequeue!(QuickHeaps.AbstractNode,pq))

To just examine the node of highest priority, call one of:

    k, v = peek(pq)
    k, v = first(pq)

Like the `dequeue!` method, the `peek` method may also be called with a type
argument.

A priority queue `pq` behaves like a dictionary:

    length(pq)                # yields number of nodes
    isempty(pq)               # yields whether priority queue is empty
    empty!(pq)                # empties priority queue
    keytype(pq)               # yields key type `K`
    valtype(pq)               # yields value type `V`
    pq[k] = v                 # set value `v` of key `k`
    push!(pq, k => v)         # idem.
    haskey(pq, k)             # yields whether key exists
    get(pq, k, def)           # query value at key, with default
    delete!(pq, k)            # delete node at key `k`
    keys(pq)                  # yields iterator over keys
    vals(pq)                  # yields iterator over values
    for (k,v) in pq; ...; end # loop over key-value pairs

Note that the 3 above iterators yield nodes in their storage order which is not
that of their priority.  The order is however the same for the 3 cases.

"""
struct PriorityQueue{K,V,T,O} <: AbstractPriorityQueue{K,V,T,O}
    order::O
    nodes::Vector{T}
    index::Dict{K,Int}
end

# Copy constructor.  The copy is independent from the original.
copy(pq::PriorityQueue{K,V,T,O}) where {K,V,T,O} =
    PriorityQueue{K,V,T,O}(ordering(pq), copy(nodes(pq)), copy(index(pq)))

"""
    FastPriorityQueue{V}([o=Forward,] T=Node{Int,V}, dims...)

yields a priority queue with ordering specified by `o::Ordering` and for nodes
of type `T<:AbstractNode{Int,V}` with `V` the type of the values encoding the
priority. Type parameter `V` may be omitted if the node type `T` is specified.

The keys in this specialized priority queue are the linear or Cartesian indices
in an array of size `dims...`.  Internally, the keys are stored as linear
indices of type `Int` but the priority queue may be indexed with linear or
Cartesian indices.  For example, if `pq` is a priority queue of this kind built
with `dims = (3,4,5)`, then all the following expressions refer to the same
key:

    pq[44]
    pq[2,3,4]
    pq[CartesianIndex(2,3,4)]

With `n` nodes, the storage of this kind of priority queue is
`prod(dims)*sizeof(Int) + n*sizeof(T)` bytes.

The method

    haskey(pq, k)

yields whether key `k` exists in priority queue `pq`.  Note that `false` is
returned, if key `k` is an invalid key (out of bounds).  This is unlike the
syntax `pq[k]` which would throw an exception in that case.

"""
struct FastPriorityQueue{V,T<:AbstractNode{Int,V},
                         O,N} <: AbstractPriorityQueue{Int,V,T,O}
    order::O
    nodes::Vector{T}
    index::Array{Int,N}
end

# Copy constructor.  The copy is independent from the original.
copy(pq::FastPriorityQueue{V,T,O,N}) where {V,T,O,N} =
    FastPriorityQueue{V,T,O,N}(ordering(pq), copy(nodes(pq)), copy(index(pq)))

# Constructors for PriorityQueue instances.

function PriorityQueue{K,V}(o::O = default_ordering(PriorityQueue),
                            ::Type{T} = Node{K,V}) where {K,V,
                                                          T<:AbstractNode{K,V},
                                                          O<:Ordering}
    return PriorityQueue{K,V,T,O}(o, T[], Dict{K,V}())
end

PriorityQueue{K,V}(::Type{T}) where {K,V,T<:AbstractNode{K,V}} =
    PriorityQueue{K,V}(default_ordering(PriorityQueue), T)

PriorityQueue{K}(::Type{T}) where {K,V,T<:AbstractNode{K,V}} =
    PriorityQueue{K,V}(T)

PriorityQueue(::Type{T}) where {K,V,T<:AbstractNode{K,V}} =
    PriorityQueue{K,V}(T)

PriorityQueue{K}(o::Ordering, ::Type{T}) where {K,V,T<:AbstractNode{K,V}} =
    PriorityQueue{K,V}(o, T)

PriorityQueue(o::Ordering, ::Type{T}) where {K,V,T<:AbstractNode{K,V}} =
    PriorityQueue{K,V}(o, T)

# Constructors for FastPriorityQueue instances.

FastPriorityQueue{V}(dims::Integer...) where {V} =
    FastPriorityQueue{V}(dims)

FastPriorityQueue{V}(o::Ordering, dims::Integer...) where {V} =
    FastPriorityQueue{V}(o, dims)

FastPriorityQueue{V}(T::Type{<:AbstractNode{Int,V}}, dims::Integer...) where {V} =
    FastPriorityQueue{V}(T, dims)

FastPriorityQueue(T::Type{<:AbstractNode{Int,<:Any}}, dims::Integer...) =
    FastPriorityQueue(T, dims)

FastPriorityQueue{V}(o::Ordering, T::Type{<:AbstractNode{Int,V}}, dims::Integer...) where {V} =
    FastPriorityQueue{V}(o, T, dims)

FastPriorityQueue(o::Ordering, T::Type{<:AbstractNode{Int,<:Any}}, dims::Integer...) =
    FastPriorityQueue(o, T, dims)

FastPriorityQueue{V}(dims::Tuple{Vararg{Integer}}) where {V} =
    FastPriorityQueue(Node{Int,V}, dims)

FastPriorityQueue{V}(o::Ordering, dims::Tuple{Vararg{Integer}}) where {V} =
    FastPriorityQueue(o, Node{Int,V}, dims)

FastPriorityQueue{V}(o::Ordering, T::Type{<:AbstractNode{Int,V}}, dims::Tuple{Vararg{Integer}}) where {V} =
    FastPriorityQueue(o, T, dims)

FastPriorityQueue{V}(T::Type{<:AbstractNode{Int,V}}, dims::Tuple{Vararg{Integer}}) where {V} =
    FastPriorityQueue(T, dims)

FastPriorityQueue(T::Type{<:AbstractNode{Int,V}}, dims::Tuple{Vararg{Integer}}) where {V} =
    FastPriorityQueue(default_ordering(FastPriorityQueue), T, dims)

FastPriorityQueue(o::O, T::Type{<:AbstractNode{Int,V}}, dims::NTuple{N,Integer}) where {O<:Ordering,V,N} =
    FastPriorityQueue{V,T,O,N}(o, T[], zeros(Int, dims))

#show(io::IO, ::MIME"text/plain", pq::AbstractPriorityQueue) =
#    print(io, "priority queue of type ", nameof(typeof(pq)),
#          " with ", length(pq), " node(s)")

show(io::IO, ::MIME"text/plain", pq::PriorityQueue{K,V}) where {K,V} =
    print(io, typename(pq), " of type ", nameof(typeof(pq)), "{", nameof(K),
          ",", nameof(V), "} with ", length(pq), " node(s)")

show(io::IO, ::MIME"text/plain", pq::FastPriorityQueue{V}) where {V} =
    print(io, typename(pq), " of type ", nameof(typeof(pq)), "{", nameof(V),
          "} with ", length(pq), " node(s)")

ordering(pq::AbstractPriorityQueue)  = getfield(pq, :order)
nodes(pq::AbstractPriorityQueue) = getfield(pq, :nodes)
index(pq::AbstractPriorityQueue) = getfield(pq, :index)

length(pq::AbstractPriorityQueue) = length(nodes(pq))

isempty(pq::AbstractPriorityQueue) = (length(pq) ≤ 0)

keytype(pq::AbstractPriorityQueue) = keytype(typeof(pq))
keytype(::Type{<:AbstractPriorityQueue{K,V}}) where {K,V} = K

valtype(pq::AbstractPriorityQueue) = valtype(typeof(pq))
valtype(::Type{<:AbstractPriorityQueue{K,V}}) where {K,V} = V

haskey(pq::AbstractPriorityQueue, key) = (heap_index(pq, key) != 0)

function get(pq::AbstractPriorityQueue, key, def)
    n = length(pq)
    if n > 0
        i = heap_index(pq, key)
        if in_range(i, n) # FIXME: Testing that i > 0 should be sufficient.
            @inbounds x = getindex(nodes(pq), i)
            return Tuple(x)
        end
    end
    return def
end

function delete!(pq::AbstractPriorityQueue, key)
    n = length(pq)
    if n > 0
        i = heap_index(pq, key)
        if in_range(i, n) # FIXME: Testing that i > 0 should be sufficient.
            A = nodes(pq)
            @inbounds k = getkey(A[i]) # key to be deleted
            if i < n
                # Replace the deleted node by the last node in the heap and
                # up-/down-heapify to restore the binary heap structure.  We
                # cannot assume that the deleted node data be accessible nor
                # valid, so we explicitely replace it before deciding in which
                # direction to go and reheapify.  Also see `unsafe_enqueue!`.
                @inbounds x = A[n] # last node
                @inbounds A[i] = x # do replace deleted node
                o = ordering(pq)
                if i ≤ 1 || lt(o, (@inbounds A[heap_parent(i)]), x)
                    unsafe_heapify_down!(pq, i, x, n - 1)
                else
                    unsafe_heapify_up!(pq, i, x)
                end
            end
            unsafe_shrink!(pq, n - 1)
            unsafe_delete_key!(pq, k)
        end
    end
    return pq
end

first(pq::AbstractPriorityQueue) = peek(pq)
peek(pq::AbstractPriorityQueue) = peek(Pair, pq)

# FIXME: Same code as for binary heaps.
function peek(T::Type, pq::AbstractPriorityQueue)
    isempty(pq) && throw_argument_error(typename(pq), " is empty")
    @inbounds x = getindex(nodes(pq), 1)
    return T(x)
end

function empty!(pq::PriorityQueue)
    empty!(nodes(pq))
    empty!(index(pq))
    return pq
end

function empty!(pq::FastPriorityQueue)
    empty!(nodes(pq))
    fill!(index(pq), 0)
    return pq
end

# Private structure used by iterators on priority queues.
struct _PriorityQueueIterator{F,T<:AbstractPriorityQueue}
    f::F
    pq::T
end

IteratorEltype(itr::_PriorityQueueIterator) = IteratorEltype(typeof(itr))
IteratorEltype(::Type{<:AbstractPriorityQueue}) = HasEltype()
eltype(itr::_PriorityQueueIterator) = eltype(typeof(itr))
function eltype(::Type{<:_PriorityQueueIterator{F,T}}
                ) where {K,V,F<:typeof(getkey),T<:AbstractPriorityQueue{K,V}}
    return K
end
function eltype(::Type{<:_PriorityQueueIterator{F,T}}
                ) where {K,V,F<:typeof(getval),T<:AbstractPriorityQueue{K,V}}
    return V
end
eltype(::Type{<:_PriorityQueueIterator}) = Any

IteratorSize(itr::_PriorityQueueIterator) = IteratorSize(typeof(itr))
IteratorSize(::Type{<:_PriorityQueueIterator}) = HasLength()
length(itr::_PriorityQueueIterator) = length(itr.pq)

# Unordered iterators.  NOTE: Both `keys` and `values` shall however return the
# elements in the same order.
function iterate(pq::AbstractPriorityQueue, i::Int = 1)
    i ≤ length(pq) || return nothing
    @inbounds x = getindex(nodes(pq), i)
    return Pair(x), i + 1
end
keys(pq::AbstractPriorityQueue) = _PriorityQueueIterator(getkey, pq)
values(pq::AbstractPriorityQueue) = _PriorityQueueIterator(getval, pq)
function Base.iterate(itr::_PriorityQueueIterator, i::Int = 1)
    i ≤ length(itr.pq) || return nothing
    @inbounds x = getindex(nodes(itr.pq), i)
    return itr.f(x), i + 1
end

pop!(pq::AbstractPriorityQueue) = dequeue!(Pair, pq)

dequeue_pair!(pq::AbstractPriorityQueue) = dequeue!(Pair, pq)

dequeue!(pq::AbstractPriorityQueue) = getkey(dequeue!(AbstractNode, pq))

# This is almost the same code as pop! for a binary heap.
function dequeue!(T::Type, pq::AbstractPriorityQueue)
    n = length(pq)
    n ≥ 1 || throw_argument_error(typename(pq), " is empty")
    A = nodes(pq)
    @inbounds x = A[1]
    if n > 1
        # Peek the last node and down-heapify starting at the root of the
        # binary heap to insert it.
        @inbounds y = A[n]
        unsafe_heapify_down!(pq, 1, y, n - 1)
    end
    unsafe_delete_key!(pq, getkey(x))
    unsafe_shrink!(pq, n - 1)
    return T(x)
end

# FIXME: Same code as for a binary heap.
push!(pq::AbstractPriorityQueue, ::Tuple{}) = pq

# FIXME: Same code as for a binary heap.
function push!(pq::AbstractPriorityQueue, args...)
    for arg in args
        push!(pq, arg)
    end
    return pq
end

push!(pq::AbstractPriorityQueue, pair::Pair) = enqueue!(pq, pair)

push!(pq::AbstractPriorityQueue{<:Any,<:Any,T}, node::T) where {T} =
    enqueue!(pq, node)

getindex(pq::AbstractPriorityQueue, ::Tuple{}) = throw_missing_key()

setindex!(pq::AbstractPriorityQueue, val, ::Tuple{}) = throw_missing_key()

throw_missing_key() = throw_argument_error("missing key")

function getindex(pq::PriorityQueue, key)
    i = heap_index(pq, key)
    # FIXME: Testing that i > 0 should be sufficient.
    in_range(i, length(pq)) || throw_argument_error(
        typename(pq), " has no node with key ", key)
    @inbounds r = getindex(nodes(pq), i)
    return getval(r)
end

setindex!(pq::PriorityQueue, val, key) = enqueue!(pq, key, val)

# Union of types that can be used to index fast priority queues.
const FastIndex = Union{Integer,CartesianIndex}

# For indexing fast priority queues, we first convert the key into a linear
# index (using the current bounds checking state).
for keytype in (:Integer, :(FastIndex...))
    @eval begin
        @inline @propagate_inbounds function getindex(pq::FastPriorityQueue,
                                                      key::$keytype)
            k = linear_index(pq, key)
            @inbounds i = getindex(index(pq), k)
            A = nodes(pq)
            if in_range(i, A)
                @inbounds x = A[i]
                return getval(x)
            end
            throw_argument_error(typename(pq), " has no node with key ",
                                 normalize_key(pq, key))
        end
        @inline @propagate_inbounds function setindex!(pq::FastPriorityQueue,
                                                       val,
                                                       key::$keytype)
            return enqueue!(pq, key, val)
        end
    end
end

normalize_key(pq::FastPriorityQueue, key::Integer) = to_int(key)
normalize_key(pq::FastPriorityQueue, key::Tuple{Vararg{FastIndex}}) =
    to_indices(index(pq), key)

"""
    to_key(pq, k)

converts the key `k` to the type suitable for priority queue `pq`.

"""
to_key(pq::AbstractPriorityQueue{K,V}, key::K) where {K,V} = key
to_key(pq::AbstractPriorityQueue{K,V}, key) where {K,V} = to_type(K, key)

"""
    to_val(pq, v)

converts the value `v` to the type suitable for priority queue `pq`.

"""
to_val(pq::AbstractPriorityQueue{K,V}, val::V) where {K,V} = val
to_val(pq::AbstractPriorityQueue{K,V}, val) where {K,V} = to_type(V, val)

"""
    to_node(pq, k, v)

converts the the key `k` and the value `v` into a node type suitable for
priority queue `pq`.

"""
to_node(pq::AbstractPriorityQueue{K,V,T}, key, val) where {K,V,T} =
    T(to_key(pq, key), to_val(pq, val))

"""
    heap_index(pq, k) -> i

yields the index in the binary heap backing the storage of the nodes of the
priority queue `pq` of the key `k`.  If the key is not in priority queue, `i =
0` is returned, otherwise `i ∈ 1:n` with `n = length(pq)` is returned.

The `heap_index` method is used to implement `haskey`, `get`, and `delete!`
methods for priority queues.  The `heap_index` method shall be specialized for
any concrete sub-types of `QuickHeaps.AbstractPriorityQueue`.

""" heap_index # NOTE: `heap_index` ~ `ht_keyindex` in `base/dict.jl`

# By default, pretend that the key is missing.
heap_index(pq::AbstractPriorityQueue, key) = 0

heap_index(pq::PriorityQueue, key) = get(index(pq), key, 0)

function heap_index(pq::FastPriorityQueue, key::Integer)
    k = to_int(key)
    I = index(pq)
    in_range(k, I) || return 0
    @inbounds i = I[k]
    return i
end

function heap_index(pq::FastPriorityQueue,
                    key::CartesianIndex)
    I = index(pq)
    if checkbounds(Bool, I, key)
        @inbounds i = I[key]
        return i
    end
    return 0
end

function heap_index(pq::FastPriorityQueue,
                    key::Tuple{Vararg{FastIndex}})
    I = index(pq)
    if checkbounds(Bool, I, key...)
        @inbounds i = I[key...]
        return i
    end
    return 0
end

"""
    linear_index(pq, k)

converts key `k` into a linear index suitable for the fast priority queue `pq`.
The key can be a linear index or a multi-dimensional index (anything accepted
by `to_indices`).  The current settings for bounds checking are used.

"""
@inline @propagate_inbounds linear_index(pq::FastPriorityQueue, key::Integer) =
    # Convert to Int, then re-call linear_index for bound checking.  Note that
    # the type assertion performed by `to_type` avoids infinite recursion.
    linear_index(pq, to_type(Int, key))

@inline function linear_index(pq::FastPriorityQueue, key::Int)
    @boundscheck checkbounds(index(pq), key)
    return key
end

@inline @propagate_inbounds function linear_index(pq::FastPriorityQueue,
                                                  key::Tuple{Vararg{FastIndex}})
    # FIXME: Shall we store the linear_indices (a small object) in the priority
    #        queue directly?
    return LinearIndices(index(pq))[key...] # also does the bound checking
end

linear_index(pq::FastPriorityQueue, key) = throw_invalid_key(pq, key)

@noinline throw_invalid_key(pq::AbstractPriorityQueue, key) = throw_argument_error(
    "invalid key of type ", typeof(key), " for ", nameof(typeof(pq)))

@noinline throw_invalid_key(pq::FastPriorityQueue, key) = throw_argument_error(
    "invalid key of type ", typeof(key), " for ", nameof(typeof(pq)),
    " expecting a linear index, an ", ndims(index(pq)),
    "-dimensional Cartesian index")

# The following is to allow the syntax enqueue!(pq, key=>val)
enqueue!(pq::AbstractPriorityQueue, pair::Pair) =
    enqueue!(pq, pair.first, pair.second)

# For a general purpose priority queue, build the node then enqueue.
enqueue!(pq::PriorityQueue, key, val) = enqueue!(pq, to_node(pq, key, val))
enqueue!(pq::PriorityQueue{K,V,T}, x::T) where {K,V,T} =
    unsafe_enqueue!(pq, x, get(index(pq), getkey(x), 0))

# For a fast priority queue, converts the key into a linear index, then enqueue.
@inline @propagate_inbounds function enqueue!(pq::FastPriorityQueue, key, val)
    k = linear_index(pq, key) # not to_key
    v = to_val(pq, val)
    x = to_node(pq, k, v)
    @inbounds i = getindex(index(pq), k)
    return unsafe_enqueue!(pq, x, i)
end

enqueue!(pq::FastPriorityQueue{V,T}, x::T) where {V,T} =
    enqueue!(pq, getkey(x), getval(x))

"""
    unsafe_enqueue!(pq, x, i) -> pq

stores node `x` in priority queue `pq` at index `i` and returns the priority
queue.  The argument `i` is an index in the binary heap backing the storage of
the nodes of the priority queue.  Index `i` is determined by the key `k` of the
node `x` and by the current state of the priority queue.  If `i` is not a valid
index in the binary heap, a new node is added; otherwise, the node at index `i`
in the binary heap is replaced by `x`.  In any cases, the binary heap is
reordered as needed.

This function is *unsafe* because it assumes that the key `k` of the node `x`
is valid (e.g. it is not out of bounds for fast priority queues) in the sense
that `I[k]` is valid for the index `I` of the priority queue.

"""
function unsafe_enqueue!(pq::AbstractPriorityQueue{K,V,T},
                         x::T, i::Int) where {K,V,T}
    A = nodes(pq)
    if in_range(i, A)
        # The key alreay exists.  Replace the node in the heap by the new node
        # and up-/down-heapify to restore the binary heap structure.  We cannot
        # assume that the replaced node data be accessible nor valid, so we
        # explicitely replace it before deciding in which direction to go and
        # reheapify.  Also see `delete!`.
        @inbounds A[i] = x # do replace deleted node
        o = ordering(pq)
        if i ≤ 1 || lt(o, (@inbounds A[heap_parent(i)]), x)
            unsafe_heapify_down!(pq, i, x)
        else
            unsafe_heapify_up!(pq, i, x)
        end
    else
        # No such key already exists.  Create a new slot at the end of the node
        # list and up-heapify to fix the structure and insert the new node.
        n = length(pq) + 1
        unsafe_heapify_up!(unsafe_grow!(pq, n), n, x)
    end
    return pq
end

"""
    unsafe_enqueue!(dir, pq, k, v) -> pq

requeues key `k` at priority `v` in priority queue `pq` forcing the
heapification of the binary heap backing the storage of the nodes of `pq` in
the direction `dir`.  A node with the same key `k` must already exists in the
queue.  If `dir = Val(:down)`, it is assumed that the new priority `v` of the
key `k` is less than the former priority; if `dir = Val(:up)`, it is assumed
that the new priority is greater than the former one.

This specialization of the `unsafe_enqueue!` method is *unsafe* because the
binary heap backing the storage of the nodes may be left with an invalid
structure if `dir` is wrong.

"""
@inline @propagate_inbounds function unsafe_enqueue!(dir::Union{Val{:down},
                                                                Val{:up}},
                                                     pq::AbstractPriorityQueue,
                                                     k, v)
    i = heap_index(pq, k)
    in_range(i, length(pq)) || throw_argument_error(
        "key ", key, " does not exists in ", typename(pq))
    @inbounds x = getindex(nodes(pq), i)
    unsafe_heapify!(dir, pq, x, i)
    return pq
end

function unsafe_heapify!(::Val{:down},
                         pq::AbstractPriorityQueue{K,V,T},
                         x::T, i::Int) where {K,V,T}
    unsafe_heapify_down!(pq, i, x)
end

function unsafe_heapify!(::Val{:up},
                         pq::AbstractPriorityQueue{K,V,T},
                         x::T, i::Int, ::Val{:up}) where {K,V,T}
    unsafe_heapify_up!(pq, i, x)
end

"""
    unsafe_grow!(pq, n) -> pq

grows the size of the binary heap backing the storage of the nodes of the
priority queue `pq` to be `n` and returns the priority queue object.

"""
unsafe_grow!(pq::Union{PriorityQueue,FastPriorityQueue}, n::Int) = begin
    resize!(nodes(pq), n)
    return pq
end

"""
    unsafe_shrink!(pq, n)

shrinks the size of the binary heap backing the storage of the nodes of the
priority queue `pq` to be `n`.

"""
unsafe_shrink!(pq::Union{PriorityQueue,FastPriorityQueue}, n::Int) =
    resize!(nodes(pq), n)

"""
    unsafe_delete_key!(pq, k)

deletes key `k` from the index of the priority queue `pq` assuming `k`
is valid.

"""
unsafe_delete_key!(pq::AbstractPriorityQueue{K}, key::K) where {K} =
    unsafe_delete_key!(index(pq), key)

# Specialized version for the type of the index.
unsafe_delete_key!(I::Array{Int}, key::Int) = @inbounds I[key] = 0
unsafe_delete_key!(I::AbstractDict, key) = delete!(I, key)

@inline function unsafe_heapify_down!(pq::AbstractPriorityQueue{K,V,T},
                                      i::Int, x::T,
                                      n::Int = length(pq)) where {K,V,T}
    o = ordering(pq)
    A = nodes(pq)
    I = index(pq)
    # FIXME: @inbounds
    begin
        while (l = heap_left(i)) ≤ n
            j = (r = heap_right(i)) > n || lt(o, A[l], A[r]) ? l : r
            lt(o, A[j], x) || break
            I[getkey(A[j])] = i
            A[i] = A[j]
            i = j
        end
        I[getkey(x)] = i
        A[i] = x
    end
end

@inline function unsafe_heapify_up!(pq::AbstractPriorityQueue{K,V,T},
                                    i::Int, x::T) where {K,V,T}
    o = ordering(pq)
    A = nodes(pq)
    I = index(pq)
    # FIXME: @inbounds
    begin
        while (j = heap_parent(i)) ≥ 1 && lt(o, x, A[j])
            I[getkey(A[j])] = i
            A[i] = A[j]
            i = j
        end
        I[getkey(x)] = i
        A[i] = x
    end
end

ext: jl


text:
using LinearAlgebra
using Fermi.DIIS

# Functions specific to DF-CTF
#
############ PREPROCESSING ################
"""
    Fermi.CoupledCluster.RCCSD{T}(Alg::CTF)

Compute a RCCSD wave function using the Compiled time factorization algorithm (CTF)
"""
function BCCD{T}(guess::BCCD{Tb},Alg::DFCTF) where { T <: AbstractFloat,
                                                    Tb <: AbstractFloat }
    refwfn = Fermi.HartreeFock.RHF()

    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    ints = refwfn.ints
    Fermi.Integrals.aux_ri!(ints)
    BCCD{T}(refwfn, guess, ints, Alg)
end

function BCCD{T}(refwfn::RHF, guess::BCCD{Tb}, ints::IntegralHelper{Tc}, Alg::DFCTF) where { T <: AbstractFloat,
                                                                                Tb <: AbstractFloat,
                                                                                Tc <: AbstractFloat }
    d = [i - a for i = diag(ints["FOO"]), a = diag(ints["FVV"])]
    D = [i + j - a - b for i = diag(ints["FOO"]), j = diag(ints["FOO"]), a = diag(ints["FVV"]), b = diag(ints["FVV"])]
    newT1 = ints["FOV"]./d
    Bov = ints["BOV"]
    @tensor oovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
    newT2 = oovv ./ D
    BCCD{T}(refwfn, ints, newT1, newT2, Alg)
end
############ END PREPROCESSING ################
#
#
#
############# KERNEL FUNCTIONS ###################

function print_bcc_alg(Alg::DFCTF)
    @output "\n    • Computing BCC with the DF-CTF algorithm .\n\n"
end

function compute_integrals(ints,Alg::DFCTF)
    Fermi.Integrals.aux_ri!(ints)
    @output "Aux basis: {}\n" ints.bname["aux"]
    ints["BOV"]
    ints["BOO"]
    ints["BVV"]
end

function delete_integrals(ints,alg::DFCTF)
    delete!(ints.cache,"BOV")
    delete!(ints.cache,"BOO")
    delete!(ints.cache,"BVV")
    delete!(ints.cache,"F")
    delete!(ints.cache,"FOO")
    delete!(ints.cache,"FVV")
    delete!(ints.cache,"FOV")
end

function compute_oovv(ints,alg::DFCTF)
    Bov = ints["BOV"]
    @tensor oovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
    oovv
end

function update_T1(T2::Array{T,4}, newT1::Array{T,2}, foo, fov, fvv, ints::IntegralHelper, alg::DFCTF) where { T <: AbstractFloat }
    newT1 .= 0
    Bov = ints["BOV"]
    Boo = ints["BOO"]
    Bvv = ints["BVV"]
    @tensor begin
        Voooo[i,j,k,l] := Boo[Q,i,k]*Boo[Q,j,l]
        Vooov[i,j,k,a] := Boo[Q,i,k]*Bov[Q,j,a]
        Voovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
        Vovov[i,a,j,b] := Boo[Q,i,j]*Bvv[Q,a,b]
    end
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x,Q=>50x) begin
        newT1[i,a] += fov[i,a]
        newT1[i,a] += 2.0*fov[k,c]*T2[i,k,a,c]
        newT1[i,a] -= fov[k,c]*T2[k,i,a,c]
        newT1[i,a] -= T2[k,i,c,d]*Bov[Q,k,d]*Bvv[Q,a,c]
        newT1[i,a] += 2.0*T2[i,k,c,d]*Bov[Q,k,d]*Bvv[Q,a,c]
        newT1[i,a] += -2.0*T2[k,l,a,c]*Vooov[k,l,i,c]
        newT1[i,a] += T2[l,k,a,c]*Vooov[k,l,i,c]
    end
end

function update_T2(T2::Array{T,4},newT2::Array{T,4},foo,fov,fvv,ints::IntegralHelper,alg::DFCTF) where T <: AbstractFloat
    newT2 .= 0
    Bov = ints["BOV"]
    Boo = ints["BOO"]
    Bvv = ints["BVV"]
    @tensor begin
        Voooo[i,j,k,l] := Boo[Q,i,k]*Boo[Q,j,l]
        Vooov[i,j,k,a] := Boo[Q,i,k]*Bov[Q,j,a]
        Voovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
        Vovov[i,a,j,b] := Boo[Q,i,j]*Bvv[Q,a,b]
    end
    dfsz,nocc,nvir = size(Bov)
    T2_inter = zeros(T,dfsz,nvir,nvir)
    T2_slice = zeros(T,nvir,nvir)
    Tslice = zeros(T,nvir,nvir)
    for i=1:nocc, j=1:nocc
        T2_inter .= 0
        T2_slice .= 0
        @views Tslice .= T2[i,j,:,:]
        Fermi.contract!(T2_inter,Tslice,Bvv,"Qad","cd","Qca")
        Fermi.contract!(T2_slice,T2_inter,Bvv,"ab","Qad","Qdb")
        newT2[i,j,:,:] += T2_slice
    end

    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>50x) begin
        newT2[i,j,a,b] += Voovv[i,j,a,b]
        newT2[i,j,a,b] += T2[k,l,a,b]*Voooo[i,j,k,l]
        newT2[i,j,a,b] += T2[k,l,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[l,k,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[i,k,a,c]*T2[l,j,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[k,i,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,a,c]*T2[l,k,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,j,a,c]*T2[k,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,j,a,c]*T2[i,l,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += 4.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,d,c]*T2[l,k,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] := -1.0*foo[i,k]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += fvv[c,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,i,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += 2.0*T2[i,k,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[i,k,a,c]*Vovov[j,c,k,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,j,a,c]*Vovov[i,c,k,b]
        P_OoVv[i,j,a,b] += T2[j,k,c,d]*T2[i,l,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T2[i,k,d,c]*T2[l,j,a,b]*Voovv[k,l,c,d]

        newT2[i,j,a,b] += P_OoVv[i,j,a,b] + P_OoVv[j,i,b,a]
    end
end
#
#
############# END KERNEL FUNCTIONS ###################


text:
# This file was generated, do not modify it. # hide
X = copy(dfX)
y = copy(dfYbinary)

coerce!(X, autotype(X, :string_to_multiclass))
yc = CategoricalArray(y[:, 1])
yc = coerce(yc, OrderedFactor)

LinearBinaryClassifierPipe = Pipeline(
    Standardizer(),
    OneHotEncoder(drop_last = true),
    LinearBinaryClassifier()
)

LogisticModel = machine(LinearBinaryClassifierPipe, X, yc)
fit!(LogisticModel)
fp = fitted_params(LogisticModel)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using Documenter, DocumenterMarkdown, PuiseuxPolynomials, LaurentPolynomials

makedocs(sitename="Mvps.jl documentation",format=Markdown(),modules=[PuiseuxPolynomials])


text:
module TestUpdatableCholeskyFactorizations
using UpdatableCholeskyFactorizations
using LinearAlgebra
using Test

@testset "updating Cholesky" begin
    element_types = (Float32, Float64, ComplexF32, ComplexF64)
    for elty in element_types
        n = 4
        A_full = randn(elty, 2n, 2n)
        A_full = A_full'A_full + I # adding diagonal for numerical stability with low-precision floats
        A = A_full[1:n, 1:n]

        UC = updatable_cholesky(A)
        @test UC isa UpdatableCholesky{elty}
        C = cholesky(A) # compare against regular cholesky
        @test issuccess(UC)
        @test UC' == UC
        @test UC.U ≈ C.U
        @test UC.L ≈ C.L
        @test UC.info == C.info
        @test size(UC) == size(A)
        for i in 1:3
            @test size(UC, i) == size(A, i)
        end
        # testing solves
        a = A_full[1:n, n]
        @test UC \ a ≈ A \ a
        @test a ≈ @view A_full[1:n, n] # checking that a is not mutated
        @test UC \ A ≈ (one(elty)*I)(n)
        @test A ≈ @view A_full[1:n, 1:n] # checking that A is not mutated
        @test a' / UC ≈ a' / A
        @test a ≈ @view A_full[1:n, n] # checking that a is not mutated
        @test A / UC ≈ (one(elty)*I)(n)
        @test A ≈ @view A_full[1:n, 1:n] # checking that A is not mutated

        # testing copy
        CUC = copy(UC)
        @test CUC.U == UC.U
        @test CUC.info == UC.info
        @test CUC.n == UC.n

        # conversion to Cholesky type
        C2 = Cholesky(UC)
        @test C2 isa Cholesky{elty}
        @test C2.U ≈ C.U
        @test C2.L ≈ C.L
        @test C2.info == C.info
        @test size(C2) == size(C)

        # conversion from cholesky type
        # TODO: test for CholeskyPivoted
        UC2 = updatable_cholesky(C2)
        @test UC2 isa UpdatableCholesky
        @test UC2.U ≈ C2.U

        # testing constructor with arbitrary factorization input
        @test Matrix(updatable_cholesky(qr(A))) ≈ A

        # adding rows and columns
        for i in 1:n
            ai = @view A_full[1:n+i, n+i]
            Ai = @view A_full[1:n+i, 1:n+i]
            add_column!(UC, ai)
            @test Matrix(UC) ≈ Ai
            @test size(UC) == (n+i, n+i)
            @test UC \ ai ≈ Ai \ ai # testing solves
            @test UC \ Ai ≈ (one(elty)*I)(n+i)
            @test ai' / UC ≈ ai' / Ai
            @test Ai / UC ≈ (one(elty)*I)(n+i)
        end
        @test AbstractMatrix(UC) ≈ A_full

        # removing rows and columns until nothing is left
        for i in reverse(1:2n)
            remove_column!(UC, i)
            @test Matrix(UC) ≈ @view A_full[1:i-1, 1:i-1]
            @test size(UC) == (i-1, i-1)
        end

        # add columns in again
        for i in 1:n
            ai = @view A_full[1:i, i]
            add_column!(UC, ai)
            @test Matrix(UC) ≈ @view A_full[1:i, 1:i]
            @test size(UC) == (i, i)
        end

        # and remove columns, this time in different order
        for i in 1:n
            remove_column!(UC, 1) # keep removing from the front
            @test Matrix(UC) ≈ @view A_full[1+i:n, 1+i:n]
            @test size(UC) == (n-i, n-i)
        end

        # testing adaptive memory management
        UC = updatable_cholesky(A, n) # only pre-allocating enough for A
        @test size(UC.U) == (n, n)
        a = @view A_full[1:n+1, n+1]
        UC = add_column!(UC, a)
        @test size(UC.U_full) == (2n, 2n) # memory got expanded
        @test Matrix(UC) ≈ @view A_full[1:n+1, 1:n+1]
    end
end

end


text:
export ChimeraLattice

struct ChimeraLattice <: AbstractLattice
    Nx::Int
    Ny::Int
end

vertices(lt::ChimeraLattice) = 1:8*lt.Nx*lt.Ny
bravais_size(cl::ChimeraLattice) = (cl.Nx, cl.Ny)
function bravais_size(cl::ChimeraLattice, i::Int)
    if i==1
        cl.Nx
    elseif i==2
        cl.Ny
    else
        DimensionMismatch("expect dimension 1 or 2, got $i")
    end
end
Base.size(lt::ChimeraLattice) = (8, lt.Nx, lt.Ny)

gap_x(lt::ChimeraLattice) = 0.5
gap_y(lt::ChimeraLattice) = 0.5
unit(lt::ChimeraLattice) = 1/(max(lt.Nx*(2+gap_x(lt)), lt.Ny*(4+gap_y(lt))))

# for beautiful printing.
function Base.getindex(lt::ChimeraLattice, k::Int, i::Int, j::Int)
    ki = (k-1) ÷ 4
    kj = mod1(k, 4)-1
    x = ((i-1)*(gap_x(lt) + 2)+ki+0.5)*unit(lt)
    y = ((j-1)*(gap_y(lt) + 4)+kj+0.5)*unit(lt)
    return (x+0.01*kj, y+0.02*ki)
end

function Base.getindex(lt::ChimeraLattice, i::Int)
    CI = CartesianIndices(size(lt))
    lt[CI[i].I...]
end

function bonds(lt::ChimeraLattice)
    bonds = Tuple{Int,Int}[]
    red_bond!(lt, 1, bonds)
    redblack_bond!(lt, 1, bonds)
    for i=2:lt.Nx
        # BLACK
        black_bond!(lt, i-1, bonds)
        # Contract with RED
        red_bond!(lt, i, bonds)
        redblack_bond!(lt, i, bonds)
    end
    return bonds
end

function red_bond!(lt, i::Int, bonds)
    CI = LinearIndices(lt |> size)
    for j=1:lt.Ny-1
        for k = 1:4
            push!(bonds, (CI[k, i, j], CI[k, i, j+1]))
        end
    end
    return bonds
end

function redblack_bond!(lt, i::Int, bonds)
    CI = LinearIndices(lt |> size)
    for j=1:lt.Ny
        for k2=5:8
            for k1=1:4
                push!(bonds, (CI[k2, i, j], CI[k1, i, j]))
            end
        end
    end
    return bonds
end

function black_bond!(lt, i::Int, bonds)
    CI = LinearIndices(lt |> size)
    for j=1:lt.Ny
        for k = 5:8
            push!(bonds, (CI[k,i,j], CI[k,i+1,j]))
        end
    end
    return bonds
end


text:
#
# Implements Real Time (RT) remote world showing. It has to track only
# differences in remote organism's world and show them in a canvas
# of a window. This is why it's different from AppRemoteWorld app.
# Second difference is requests direction - from server to client,
# but not otherwise (like in AppRemoteWorld).
#
# @author DeadbraiN
#
include("global/ImportFolders.jl")
import RemoteWorldRT
#
# Analog of main() in C lang
#
function main()
  local cfg::Config.ConfigData = Config.create()

  if Config.isEmpty(cfg) quit() end
  rw = RemoteWorldRT.create(cfg, cfg.conServerIp, cfg.conServerPort, cfg.conFastServerPort)
  RemoteWorldRT.start(rw)
end
#
# Application entry point
#
main()

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# cartesian distance

distance(a::ENU, b::ENU) = _distance(a.east, a.north, a.up,
                                     b.east, b.north, b.up)

distance(a::ECEF, b::ECEF) = _distance(a.x, a.y, a.z,
                                       b.x, b.y, b.z)

function _distance(x1, y1, z1, x2, y2, z2)
    Δx = x2 - x1
    Δy = y2 - y1
    Δz = z2 - z1

    return hypot(hypot(Δx,  Δy), Δz)
end

# haversine (spherical)

function haversine_distance{T}(a::LL{T}, b::LL{T})
    r = 6_371_009 # reduces average rather than maximum error

    # Using deg2rad is faster and a tiny bit less accurate than sind/cosd --
    #   already very rough, might as well keep it fast.
    ϕ₁, λ₁ = deg2rad(a.lat), deg2rad(a.lon)
    ϕ₂, λ₂ = deg2rad(b.lat), deg2rad(b.lon)

    hϕ = haversin(ϕ₂ - ϕ₁)
    hλ = haversin(λ₂ - λ₁)

    h = hϕ + cos(ϕ₁)*cos(ϕ₂)*hλ

    return 2r * asin(sqrt(h))
end

function haversin(θ)
    s = sin(θ/2)
    return s*s
end

# vicenty's (ellipsoidal)

distance{T}(a::LL{T}, b::LL{T}) = vicentys_inverse(a, b)[1]


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

abstract type AbstractCmd end

# libuv process option flags
const UV_PROCESS_WINDOWS_VERBATIM_ARGUMENTS = UInt8(1 << 2)
const UV_PROCESS_DETACHED = UInt8(1 << 3)
const UV_PROCESS_WINDOWS_HIDE = UInt8(1 << 4)

struct Cmd <: AbstractCmd
    exec::Vector{String}
    ignorestatus::Bool
    flags::UInt32 # libuv process flags
    env::Union{Vector{String},Nothing}
    dir::String
    cpus::Union{Nothing,Vector{UInt16}}
    Cmd(exec::Vector{String}) =
        new(exec, false, 0x00, nothing, "", nothing)
    Cmd(cmd::Cmd, ignorestatus, flags, env, dir, cpus = nothing) =
        new(cmd.exec, ignorestatus, flags, env,
            dir === cmd.dir ? dir : cstr(dir), cpus)
    function Cmd(cmd::Cmd; ignorestatus::Bool=cmd.ignorestatus, env=cmd.env, dir::AbstractString=cmd.dir,
                 cpus::Union{Nothing,Vector{UInt16}} = cmd.cpus,
                 detach::Bool = 0 != cmd.flags & UV_PROCESS_DETACHED,
                 windows_verbatim::Bool = 0 != cmd.flags & UV_PROCESS_WINDOWS_VERBATIM_ARGUMENTS,
                 windows_hide::Bool = 0 != cmd.flags & UV_PROCESS_WINDOWS_HIDE)
        flags = detach * UV_PROCESS_DETACHED |
                windows_verbatim * UV_PROCESS_WINDOWS_VERBATIM_ARGUMENTS |
                windows_hide * UV_PROCESS_WINDOWS_HIDE
        new(cmd.exec, ignorestatus, flags, byteenv(env),
            dir === cmd.dir ? dir : cstr(dir), cpus)
    end
end

has_nondefault_cmd_flags(c::Cmd) =
    c.ignorestatus ||
    c.flags != 0x00 ||
    c.env !== nothing ||
    c.dir !== "" ||
    c.cpus !== nothing

"""
    Cmd(cmd::Cmd; ignorestatus, detach, windows_verbatim, windows_hide, env, dir)

Construct a new `Cmd` object, representing an external program and arguments, from `cmd`,
while changing the settings of the optional keyword arguments:

* `ignorestatus::Bool`: If `true` (defaults to `false`), then the `Cmd` will not throw an
  error if the return code is nonzero.
* `detach::Bool`: If `true` (defaults to `false`), then the `Cmd` will be run in a new
  process group, allowing it to outlive the `julia` process and not have Ctrl-C passed to
  it.
* `windows_verbatim::Bool`: If `true` (defaults to `false`), then on Windows the `Cmd` will
  send a command-line string to the process with no quoting or escaping of arguments, even
  arguments containing spaces. (On Windows, arguments are sent to a program as a single
  "command-line" string, and programs are responsible for parsing it into arguments. By
  default, empty arguments and arguments with spaces or tabs are quoted with double quotes
  `"` in the command line, and `\\` or `"` are preceded by backslashes.
  `windows_verbatim=true` is useful for launching programs that parse their command line in
  nonstandard ways.) Has no effect on non-Windows systems.
* `windows_hide::Bool`: If `true` (defaults to `false`), then on Windows no new console
  window is displayed when the `Cmd` is executed. This has no effect if a console is
  already open or on non-Windows systems.
* `env`: Set environment variables to use when running the `Cmd`. `env` is either a
  dictionary mapping strings to strings, an array of strings of the form `"var=val"`, an
  array or tuple of `"var"=>val` pairs. In order to modify (rather than replace) the
  existing environment, initialize `env` with `copy(ENV)` and then set `env["var"]=val` as
  desired.  To add to an environment block within a `Cmd` object without replacing all
  elements, use [`addenv()`](@ref) which will return a `Cmd` object with the updated environment.
* `dir::AbstractString`: Specify a working directory for the command (instead
  of the current directory).

For any keywords that are not specified, the current settings from `cmd` are used. Normally,
to create a `Cmd` object in the first place, one uses backticks, e.g.

    Cmd(`echo "Hello world"`, ignorestatus=true, detach=false)
"""
Cmd

hash(x::Cmd, h::UInt) = hash(x.exec, hash(x.env, hash(x.ignorestatus, hash(x.dir, hash(x.flags, h)))))
==(x::Cmd, y::Cmd) = x.exec == y.exec && x.env == y.env && x.ignorestatus == y.ignorestatus &&
                     x.dir == y.dir && isequal(x.flags, y.flags)

struct OrCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    OrCmds(a::AbstractCmd, b::AbstractCmd) = new(a, b)
end

struct ErrOrCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    ErrOrCmds(a::AbstractCmd, b::AbstractCmd) = new(a, b)
end

struct AndCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    AndCmds(a::AbstractCmd, b::AbstractCmd) = new(a, b)
end

hash(x::AndCmds, h::UInt) = hash(x.a, hash(x.b, h))
==(x::AndCmds, y::AndCmds) = x.a == y.a && x.b == y.b

shell_escape(cmd::Cmd; special::AbstractString="") =
    shell_escape(cmd.exec..., special=special)
shell_escape_posixly(cmd::Cmd) =
    shell_escape_posixly(cmd.exec...)
shell_escape_csh(cmd::Cmd) =
    shell_escape_csh(cmd.exec...)
escape_microsoft_c_args(cmd::Cmd) =
    escape_microsoft_c_args(cmd.exec...)
escape_microsoft_c_args(io::IO, cmd::Cmd) =
    escape_microsoft_c_args(io::IO, cmd.exec...)

function show(io::IO, cmd::Cmd)
    print_env = cmd.env !== nothing
    print_dir = !isempty(cmd.dir)
    (print_env || print_dir) && print(io, "setenv(")
    print_cpus = cmd.cpus !== nothing
    print_cpus && print(io, "setcpuaffinity(")
    print(io, '`')
    join(io, map(cmd.exec) do arg
        replace(sprint(context=io) do io
            with_output_color(:underline, io) do io
                print_shell_word(io, arg, shell_special)
            end
        end, '`' => "\\`")
    end, ' ')
    print(io, '`')
    if print_cpus
        print(io, ", ")
        show(io, collect(Int, something(cmd.cpus)))
        print(io, ")")
    end
    print_env && (print(io, ","); show(io, cmd.env))
    print_dir && (print(io, "; dir="); show(io, cmd.dir))
    (print_dir || print_env) && print(io, ")")
    nothing
end

function show(io::IO, cmds::Union{OrCmds,ErrOrCmds})
    print(io, "pipeline(")
    show(io, cmds.a)
    print(io, ", ")
    print(io, isa(cmds, ErrOrCmds) ? "stderr=" : "stdout=")
    show(io, cmds.b)
    print(io, ")")
end

function show(io::IO, cmds::AndCmds)
    show(io, cmds.a)
    print(io, " & ")
    show(io, cmds.b)
end

const STDIN_NO  = 0
const STDOUT_NO = 1
const STDERR_NO = 2

struct FileRedirect
    filename::String
    append::Bool
    FileRedirect(filename::AbstractString, append::Bool) = FileRedirect(convert(String, filename), append)
    function FileRedirect(filename::String, append::Bool)
        if lowercase(filename) == (@static Sys.iswindows() ? "nul" : "/dev/null")
            @warn "For portability use devnull instead of a file redirect" maxlog=1
        end
        return new(filename, append)
    end
end

# setup_stdio ≈ cconvert
# rawhandle ≈ unsafe_convert
rawhandle(::DevNull) = C_NULL
rawhandle(x::OS_HANDLE) = x
if OS_HANDLE !== RawFD
    rawhandle(x::RawFD) = Libc._get_osfhandle(x)
end
setup_stdio(stdio::Union{DevNull,OS_HANDLE,RawFD}, ::Bool) = (stdio, false)

const Redirectable = Union{IO, FileRedirect, RawFD, OS_HANDLE}
const StdIOSet = NTuple{3, Redirectable}

struct CmdRedirect <: AbstractCmd
    cmd::AbstractCmd
    handle::Redirectable
    stream_no::Int
    readable::Bool
end
CmdRedirect(cmd, handle, stream_no) = CmdRedirect(cmd, handle, stream_no, stream_no == STDIN_NO)

function show(io::IO, cr::CmdRedirect)
    print(io, "pipeline(")
    show(io, cr.cmd)
    print(io, ", ")
    if cr.stream_no == STDOUT_NO
        print(io, "stdout")
    elseif cr.stream_no == STDERR_NO
        print(io, "stderr")
    elseif cr.stream_no == STDIN_NO
        print(io, "stdin")
    else
        print(io, cr.stream_no)
    end
    print(io, cr.readable ? "<" : ">")
    show(io, cr.handle)
    print(io, ")")
end

"""
    ignorestatus(command)

Mark a command object so that running it will not throw an error if the result code is non-zero.
"""
ignorestatus(cmd::Cmd) = Cmd(cmd, ignorestatus=true)
ignorestatus(cmd::Union{OrCmds,AndCmds}) =
    typeof(cmd)(ignorestatus(cmd.a), ignorestatus(cmd.b))

"""
    detach(command)

Mark a command object so that it will be run in a new process group, allowing it to outlive the julia process, and not have Ctrl-C interrupts passed to it.
"""
detach(cmd::Cmd) = Cmd(cmd; detach=true)

# like String(s), but throw an error if s contains NUL, since
# libuv requires NUL-terminated strings
function cstr(s)
    if Base.containsnul(s)
        throw(ArgumentError("strings containing NUL cannot be passed to spawned processes"))
    end
    return String(s)
end

# convert various env representations into an array of "key=val" strings
byteenv(env::AbstractArray{<:AbstractString}) =
    String[cstr(x) for x in env]
byteenv(env::AbstractDict) =
    String[cstr(string(k)*"="*string(v)) for (k,v) in env]
byteenv(env::Nothing) = nothing
byteenv(env::Union{AbstractVector{Pair{T,V}}, Tuple{Vararg{Pair{T,V}}}}) where {T<:AbstractString,V} =
    String[cstr(k*"="*string(v)) for (k,v) in env]

"""
    setenv(command::Cmd, env; dir)

Set environment variables to use when running the given `command`. `env` is either a
dictionary mapping strings to strings, an array of strings of the form `"var=val"`, or
zero or more `"var"=>val` pair arguments. In order to modify (rather than replace) the
existing environment, create `env` through `copy(ENV)` and then setting `env["var"]=val`
as desired, or use [`addenv`](@ref).

The `dir` keyword argument can be used to specify a working directory for the command.
`dir` defaults to the currently set `dir` for `command` (which is the current working
directory if not specified already).

See also [`Cmd`](@ref), [`addenv`](@ref), [`ENV`](@ref), [`pwd`](@ref).
"""
setenv(cmd::Cmd, env; dir=cmd.dir) = Cmd(cmd; env=byteenv(env), dir=dir)
setenv(cmd::Cmd, env::Pair{<:AbstractString}...; dir=cmd.dir) =
    setenv(cmd, env; dir=dir)
setenv(cmd::Cmd; dir=cmd.dir) = Cmd(cmd; dir=dir)

# split environment entry string into before and after first `=` (key and value)
function splitenv(e::String)
    i = findnext('=', e, 2)
    if i === nothing
        throw(ArgumentError("malformed environment entry"))
    end
    e[1:prevind(e, i)], e[nextind(e, i):end]
end

"""
    addenv(command::Cmd, env...; inherit::Bool = true)

Merge new environment mappings into the given [`Cmd`](@ref) object, returning a new `Cmd` object.
Duplicate keys are replaced.  If `command` does not contain any environment values set already,
it inherits the current environment at time of `addenv()` call if `inherit` is `true`.
Keys with value `nothing` are deleted from the env.

See also [`Cmd`](@ref), [`setenv`](@ref), [`ENV`](@ref).

!!! compat "Julia 1.6"
    This function requires Julia 1.6 or later.
"""
function addenv(cmd::Cmd, env::Dict; inherit::Bool = true)
    new_env = Dict{String,String}()
    if cmd.env === nothing
        if inherit
            merge!(new_env, ENV)
        end
    else
        for (k, v) in splitenv.(cmd.env)
            new_env[string(k)::String] = string(v)::String
        end
    end
    for (k, v) in env
        if v === nothing
            delete!(new_env, string(k)::String)
        else
            new_env[string(k)::String] = string(v)::String
        end
    end
    return setenv(cmd, new_env)
end

function addenv(cmd::Cmd, pairs::Pair{<:AbstractString}...; inherit::Bool = true)
    return addenv(cmd, Dict(k => v for (k, v) in pairs); inherit)
end

function addenv(cmd::Cmd, env::Vector{<:AbstractString}; inherit::Bool = true)
    return addenv(cmd, Dict(k => v for (k, v) in splitenv.(env)); inherit)
end

"""
    setcpuaffinity(original_command::Cmd, cpus) -> command::Cmd

Set the CPU affinity of the `command` by a list of CPU IDs (1-based) `cpus`.  Passing
`cpus = nothing` means to unset the CPU affinity if the `original_command` has any.

This function is supported only in Linux and Windows.  It is not supported in macOS because
libuv does not support affinity setting.

!!! compat "Julia 1.8"
    This function requires at least Julia 1.8.

# Examples

In Linux, the `taskset` command line program can be used to see how `setcpuaffinity` works.

```julia
julia> run(setcpuaffinity(`sh -c 'taskset -p \$\$'`, [1, 2, 5]));
pid 2273's current affinity mask: 13
```

Note that the mask value `13` reflects that the first, second, and the fifth bits (counting
from the least significant position) are turned on:

```julia
julia> 0b010011
0x13
```
"""
function setcpuaffinity end
setcpuaffinity(cmd::Cmd, ::Nothing) = Cmd(cmd; cpus = nothing)
setcpuaffinity(cmd::Cmd, cpus) = Cmd(cmd; cpus = collect(UInt16, cpus))

(&)(left::AbstractCmd, right::AbstractCmd) = AndCmds(left, right)
redir_out(src::AbstractCmd, dest::AbstractCmd) = OrCmds(src, dest)
redir_err(src::AbstractCmd, dest::AbstractCmd) = ErrOrCmds(src, dest)

# Stream Redirects
redir_out(dest::Redirectable, src::AbstractCmd) = CmdRedirect(src, dest, STDIN_NO)
redir_out(src::AbstractCmd, dest::Redirectable) = CmdRedirect(src, dest, STDOUT_NO)
redir_err(src::AbstractCmd, dest::Redirectable) = CmdRedirect(src, dest, STDERR_NO)

# File redirects
redir_out(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, false), STDOUT_NO)
redir_out(src::AbstractString, dest::AbstractCmd) = CmdRedirect(dest, FileRedirect(src, false), STDIN_NO)
redir_err(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, false), STDERR_NO)
redir_out_append(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, true), STDOUT_NO)
redir_err_append(src::AbstractCmd, dest::AbstractString) = CmdRedirect(src, FileRedirect(dest, true), STDERR_NO)

"""
    pipeline(command; stdin, stdout, stderr, append=false)

Redirect I/O to or from the given `command`. Keyword arguments specify which of the
command's streams should be redirected. `append` controls whether file output appends to the
file. This is a more general version of the 2-argument `pipeline` function.
`pipeline(from, to)` is equivalent to `pipeline(from, stdout=to)` when `from` is a command,
and to `pipeline(to, stdin=from)` when `from` is another kind of data source.

**Examples**:

```julia
run(pipeline(`dothings`, stdout="out.txt", stderr="errs.txt"))
run(pipeline(`update`, stdout="log.txt", append=true))
```
"""
function pipeline(cmd::AbstractCmd; stdin=nothing, stdout=nothing, stderr=nothing, append::Bool=false)
    if append && stdout === nothing && stderr === nothing
        throw(ArgumentError("append set to true, but no output redirections specified"))
    end
    if stdin !== nothing
        cmd = redir_out(stdin, cmd)
    end
    if stdout !== nothing
        cmd = append ? redir_out_append(cmd, stdout) : redir_out(cmd, stdout)
    end
    if stderr !== nothing
        cmd = append ? redir_err_append(cmd, stderr) : redir_err(cmd, stderr)
    end
    return cmd
end

pipeline(cmd::AbstractCmd, dest) = pipeline(cmd, stdout=dest)
pipeline(src::Union{Redirectable,AbstractString}, cmd::AbstractCmd) = pipeline(cmd, stdin=src)

"""
    pipeline(from, to, ...)

Create a pipeline from a data source to a destination. The source and destination can be
commands, I/O streams, strings, or results of other `pipeline` calls. At least one argument
must be a command. Strings refer to filenames. When called with more than two arguments,
they are chained together from left to right. For example, `pipeline(a,b,c)` is equivalent to
`pipeline(pipeline(a,b),c)`. This provides a more concise way to specify multi-stage
pipelines.

**Examples**:

```julia
run(pipeline(`ls`, `grep xyz`))
run(pipeline(`ls`, "out.txt"))
run(pipeline("out.txt", `grep xyz`))
```
"""
pipeline(a, b, c, d...) = pipeline(pipeline(a, b), c, d...)

## implementation of `cmd` syntax ##

cmd_interpolate(xs...) = cstr(string(map(cmd_interpolate1, xs)...))
cmd_interpolate1(x) = x
cmd_interpolate1(::Nothing) = throw(ArgumentError("`nothing` can not be interpolated into commands (`Cmd`)"))

arg_gen() = String[]
arg_gen(x::AbstractString) = String[cstr(x)]
function arg_gen(cmd::Cmd)
    if has_nondefault_cmd_flags(cmd)
        throw(ArgumentError("Non-default environment behavior is only permitted for the first interpolant."))
    end
    cmd.exec
end

function arg_gen(head)
    if isiterable(typeof(head))
        vals = String[]
        for x in head
            push!(vals, cmd_interpolate(x))
        end
        return vals
    else
        return String[cmd_interpolate(head)]
    end
end

function arg_gen(head, tail...)
    head = arg_gen(head)
    tail = arg_gen(tail...)
    vals = String[]
    for h = head, t = tail
        push!(vals, cmd_interpolate(h,t))
    end
    return vals
end

function cmd_gen(parsed)
    args = String[]
    if length(parsed) >= 1 && isa(parsed[1], Tuple{Cmd})
        cmd = (parsed[1]::Tuple{Cmd})[1]
        (ignorestatus, flags, env, dir) = (cmd.ignorestatus, cmd.flags, cmd.env, cmd.dir)
        append!(args, cmd.exec)
        for arg in tail(parsed)
            append!(args, arg_gen(arg...)::Vector{String})
        end
        return Cmd(Cmd(args), ignorestatus, flags, env, dir)
    else
        for arg in parsed
            append!(args, arg_gen(arg...)::Vector{String})
        end
        return Cmd(args)
    end
end

"""
    @cmd str

Similar to `cmd`, generate a `Cmd` from the `str` string which represents the shell command(s) to be executed.
The [`Cmd`](@ref) object can be run as a process and can outlive the spawning julia process (see `Cmd` for more).

# Examples
```jldoctest
julia> cm = @cmd " echo 1 "
`echo 1`

julia> run(cm)
1
Process(`echo 1`, ProcessExited(0))
```
"""
macro cmd(str)
    cmd_ex = shell_parse(str, special=shell_special, filename=String(__source__.file))[1]
    return :(cmd_gen($(esc(cmd_ex))))
end

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using Base.Threads

filesBuffer=Dict()
struct logfile
    fileLock::SpinLock
    fileObj
    filename::String
    function logfile( name::String )
        global filesBuffer
        fileLock = SpinLock()
        fileObj = open( name,"w")
        filename =name
        new(fileLock,fileObj,filename  )
    end
end

function writeLog(logObj::logfile,str::String )
    lock(logObj.fileLock)
        print(logObj.fileObj,str )
        print(logObj.fileObj,"\r\n" )
        flush(logObj.fileObj )
    unlock(logObj.fileLock)
      @warn  str , Base.Threads.threadid()
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function ub_smote{S<:Real}(X::Array{S, 2}, y, pct_over = 200, k = 5, pct_under = 200)

    dat = hcat(X, y)
    minority_indcs = find(y .== 1.0)
    n, p = size(dat)

    new_exs = smote_exs(dat[minority_indcs, :], p, pct_over, k)
    majority_indcs = setdiff(1:size(X, 1), minority_indcs)

    # Get the undersample of the "majority class" examples
    sel_majority = sample(majority_indcs,
                          round(Int, (pct_under/100) * size(new_exs, 1)),
                          replace = true)

    # Final dataset (the undersample + the rare cases + the smoted exs)
    newdata = vcat(dat[sel_majority, :], dat[minority_indcs, :], new_exs)
    n_new = size(newdata, 1)

    # Shuffle the order of instances
    newdata = newdata[sample(1:n_new, n_new, replace = false), :]

    X_new = newdata[:, 1:(end-1)]
    y_new = newdata[:, end]

    return (X_new, y_new)
end

# n_change(n, n1) = 0.5n - n1                        # give us number of new positive cases needed for balanced data

function smote(X, y, k::Int)
    pct_minority = pct_needed(y)
    # println("Percent oversampling: $pct_minority")
    n = length(y)
    n1 = count(x -> x == 1.0, y)
    needed = n - (n1 + n1*(pct_minority/100))
    # println("Needed: $needed")
    pct_majority = needed/(n1 * pct_minority/100) * 100

    # println("Percent undersampling: $pct_majority")
    X_new, y_new = ub_smote(X, y, pct_minority, k, pct_majority)
    return (X_new, y_new)
end

# n = 120
# n0 = 60
# X = randn(n, 10);
# y = vcat(zeros(n0), ones(n - n0));
# X2, y2 = smote(X, y, 5)
#
# length(y2)
# countmap(y2)


text:
## Let's see how finite difference approximations work.

using Printf
myfun = x -> sin(x)
myderiv = x -> cos(x) # so we can get error
finitediff(f,x,h) = (f(x+h) - f(x))/h
x = pi/6
truederiv = myderiv(x)
for h = 10.0.^-(1:16)
  fdderiv = finitediff(myfun,x,h)
  @printf("h=%.1e approx=%.8e err=%.8e\n",
    h, fdderiv, abs(fdderiv-truederiv))
end

## Now try a multidim version
function finitediff_i(f,x,h,i)
  fx = f(x)
  xi = x[i] # make a copy
  x[i] += h
  fxh = f(x) # eval f(x+h)
  x[i] = xi # restore
  return (fxh-fx)/h
end
finitediffgrad(f,x,h) = map( i -> finitediff_i(f,x,h,i), 1:length(x) )
## validate our work on gradient of quadratics
using LinearAlgebra
Q = [1 0.5;0.5 1]
s = [-1,-1.5]
quadratic2(A,b) = x -> x'*A*x + x'*b
myfun = quadratic2(Q,s)
myderiv = x -> 2*Q*x + s
x = [3.5;-1.5]
truederiv = myderiv(x)
for h = 10.0.^-(1:16)
  fdderiv = finitediffgrad(myfun,x,h)
  @printf("h=%.1e  err=%.8e\n",
    h, norm(fdderiv-truederiv))
end


text:
"""
    Parabolic{T <: Number} <: Construct{T}

A parabolic Cayley-Dickson construct as a pair of `Number`s.
"""
struct Parabolic{T <: Number} <: Construct{T}
    l::T
    r::T
end

function (==)(x::Parabolic, y::Parabolic)
    x.l == y.l && x.r == y.r
end

function real(z::Parabolic)
    real(z.l)
end

"""
    quadrance(z::CayleyDickson.Parabolic)

For parabolic Cayley-Dickson constructs, return the squared absolute value of the
left component of `z`. Depending on the type of the components of `z`, this can
either be non-negative, or positive, negative or zero.
"""
function quadrance(z::Parabolic)
    quadrance(z.l)
end

function iszerodivisor(z::Parabolic)
    iszerodivisor(z.l)
end

function zero(::Type{Parabolic{T}}) where T <: Number
    Parabolic{T}(zero(T), zero(T))
end

function zero(::Parabolic{T}) where T <: Number
    zero(Parabolic{T})
end

function one(::Type{Parabolic{T}}) where T <: Number
    Parabolic{T}(one(T), zero(T))
end

function one(::Parabolic{T}) where T <: Number
    one(Parabolic{T})
end

function (-)(z::Parabolic)
    Parabolic(-(z.l), -(z.r))
end

function conj(z::Parabolic)
    Parabolic(conj(z.l), -(z.r))
end

function (+)(x::Parabolic{T}, y::Parabolic{T}) where T <: Number
    Parabolic{T}(x.l + y.l, x.r + y.r)
end

function (+)(a::T, z::Parabolic{T}) where T <: Number
    Parabolic{T}(a + z.l, z.r)
end

function (+)(z::Parabolic{T}, a::T) where T <: Number
    Parabolic{T}(z.l + a, z.r)
end

function (+)(a::Real, z::Parabolic)
    l = a + z.l
    Parabolic(l, oftype(l, z.r))
end

function (+)(z::Parabolic, a::Real)
    l = z.l + a
    Parabolic(l, oftype(l, z.r))
end

function (-)(x::Parabolic{T}, y::Parabolic{T}) where T <: Number
    Parabolic{T}(x.l - y.l, x.r - y.r)
end

function (-)(a::T, z::Parabolic{T}) where T <: Number
    Parabolic{T}(a - z.l, -(z.r))
end

function (-)(z::Parabolic{T}, a::T) where T <: Number
    Parabolic{T}(z.l - a, z.r)
end

function (-)(a::Real, z::Parabolic)
    l = a - z.l
    Parabolic(l, -(oftype(l, z.r)))
end

function (-)(z::Parabolic, a::Real)
    l = z.l - a
    Parabolic(l, oftype(l, z.r))
end

function (*)(x::Parabolic{T}, y::Parabolic{T}) where T <: Number
    Parabolic{T}(
        (x.l * y.l),
        (y.r * x.l) + (x.r * conj(y.l))
    )
end

function (*)(a::T, z::Parabolic{T}) where T <: Number
    Parabolic{T}(a * z.l, z.r * a)
end

function (*)(z::Parabolic{T}, a::T) where T <: Number
    Parabolic{T}(z.l * a, z.r * conj(a))
end

function (*)(a::Real, z::Parabolic)
    Parabolic(a * z.l, a * z.r)
end

function (*)(z::Parabolic, a::Real)
    Parabolic(z.l * a, z.r * a)
end

function (/)(z::Parabolic, a::Real)
    if iszero(a)
        error(ZeroDenominator)
    end

    Parabolic(z.l / a, z.r / a)
end

function (\)(a::Real, z::Parabolic)
    if iszero(a)
        error(ZeroDenominator)
    end

    Parabolic(a \ z.l, a \ z.r)
end

function widen(::Type{Parabolic{T}}) where T
    Parabolic{widen(T)}
end

function float(::Type{Parabolic{T}}) where T <: AbstractFloat
    Parabolic{T}
end

function float(::Type{Parabolic{T}}) where T
    Parabolic{float(T)}
end

function cloak(z::Parabolic)
    Parabolic(cloak(z.l), dagger(z.r))
end

function dagger(z::Parabolic)
    Parabolic(dagger(z.l), cloak(z.r))
end

function star(z::Parabolic)
    Parabolic(dagger(star(z.r)), star(z.l))
end


text:
export Environment

# dependent package
#using DelimitedFiles

function Environment(option::Options,
                    source::Geometry,
                    receiver::Geometry,
                    ssp::Atmosphere,
                    terrain::Geometry,
                    reflection::Boundary)

    fn = option.filename
    filename = "$fn.env" # environment file

    freq = source.frequency
    nmedia = 1
    opt1 = option.option1
    opt2 = option.option2
    botdepth = maximum(ssp.sound_speed_profile.x)

    open(filename,"w") do io
        write(io,"'$fn'         ! title\n")
        write(io,"$freq         ! frequency (Hz)\n")
        write(io,"$nmedia       ! nmedia\n")
        write(io,"'$opt1'       ! C-linear, F-file\n") # option 1
        write(io,"51 0.0 $botdepth ! depth of bottom\n")

        ###### Sound speed profile
        for i=1:size(ssp.sound_speed_profile.x,1)
            zi = ssp.sound_speed_profile.x[i]
            vi = ssp.sound_speed_profile.z[i]
            #println(zi)
            write(io,"$zi $vi /\n")
        end
        ###### End read sound speed

        write(io,"'$opt2' 0.0 \n") # option 2

        ##---> Source information
        nsources = 1
        write(io,"$nsources ! number of sources \n")
        source_height = source.height[1]
        write(io,"$source_height / ! SD (1:NSD), m\n")

        ##----> Receiver information
        nrd = receiver.depth_point
        rd = receiver.depth
        nr = receiver.range_point
        r = receiver.range

        write(io,"$nrd      ! number of receivers in depth \n")
        if nrd==1
            rd1 = rd.x[1]
            write(io,"$rd1 / ! RD(1:RD) \n")
        else
            rdx = rd.x[1]
            rdz = rd.z[1]
            write(io,"$rdx $rdz / ! RD(1:RD) \n")
        end
        write(io,"$nr       ! number of receivers in range \n")
        rx = r.x[1]
        rz = r.z[1]
        write(io,"$rx $rz /  ! R(1:NR) \n")

        ## option 3
        opt3 = option.analyse
        nbeam = option.num_ray

        write(io,"'$opt3'   ! R/C/I/S\n") # option3
        write(io,"$nbeam    ! number of ray \n")
        neg_alpha = option.alpha.x[1]
        pos_alpha = option.alpha.z[1]
        write(io,"$neg_alpha $pos_alpha / ! alpha, deg \n")

        zmax = option.box.z[1] + 0.1*(option.box.z[1])
        rmax = option.box.x[1] + 0.1*(option.box.x[1])
        raystep = option.step
        write(io,"$raystep $zmax $rmax ! step (m), zbox(m), rbox(km) \n")

    end # end write environment file

    ### Write battery files
    filename2 = "$fn.bty" # battery file
    interp_tp = terrain.interp_type
    numpoint = size(terrain.profile.x,1)

    open(filename2,"w") do io
        write(io,"'$interp_tp' \n")
        write(io,"$numpoint   \n")
        for i=1:numpoint
            zi = terrain.profile.x[i]
            vi = terrain.profile.z[i]
            #println(zi)
            write(io,"$zi $vi \n")
        end
    end # end write battery file

    ### Write top and bottom reflection coefficient files
    filename3 = "$fn.brc" # bottom reflection
    num_bot = size(reflection.bottom_coeff.x,1)
    open(filename3,"w") do io
        write(io,"$num_bot   \n")
        for i=1:num_bot
            xi = reflection.bottom_coeff.x[i]
            yi = reflection.bottom_coeff.y[i]
            zi = reflection.bottom_coeff.z[i]
            #println(zi)
            write(io,"$xi $yi $zi \n")
        end
    end # end write bottom reflection coefficient file

    filename4 = "$fn.trc" # top reflection
    num_top = size(reflection.top_coeff.x,1)
    open(filename4,"w") do io
        write(io,"$num_top   \n")
        for i=1:num_top
            xi = reflection.top_coeff.x[i]
            yi = reflection.top_coeff.y[i]
            zi = reflection.top_coeff.z[i]
            #println(zi)
            write(io,"$xi $yi $zi \n")
        end
    end # end write top reflection coefficient file

end


text:
# vocoder.jl
# Compare Matlab and Sound versions of phase vocoder.

using Sound: phase_vocoder, hann, soundsc
using MATLAB
using Test: @testset, @test

ENV["MATLAB_ROOT"] = "/Applications/freeware/matlab"

@testset "vocoder" begin
end
    S = 44100
    N = S ÷ 5
    hopin = 121.
#   hopout = hopin * 2
    hopout = hopin * 1.5 # todo: large differences
# for hopin in (100., 121.), hopout in hopin .* (1, 1.5, 0.5)
@show hopin, hopout
    f = 540
    x = sin.(2π*f*(1:N)/S)
    yj = phase_vocoder(x, S; hopin, hopout)

    @show "start matlab"
    Nd = Float64(N)
    x2 = reshape(x, :, 1)
    mat"ym = pv($x2, $Nd, $hopin, $hopout)"
    @mget ym
    ym = Float32.(vec(ym))

    @show extrema(yj)
    @show extrema(ym)
    @show extrema(yj - ym)
    #=
    =#
#   @test isapprox(yj, ym; rtol = 0.01) # 1% seems OK for now?
# end

#=
=#
    # plot spectrograms to compare
    using Plots
    using MIRTjim: jim
    using DSP: spectrogram

    plot(yj - ym)
    xlims = (5000,5400)
p0 = plot(
    plot(yj, marker=:circle; xlims),
    plot(ym, marker=:circle; xlims),
    layout=(1,2),
); gui()
    sj = spectrogram(yj / maximum(abs, x); fs=S)
    sm = spectrogram(ym / maximum(abs, x); fs=S)
    dB = x -> max(log10.(x), -10)
    fun = (s) -> jim(s.time, s.freq, dB.(s.power)', aspect_ratio=:auto,
        ylims = (0,1800), color=:viridis, clim = (-10,0),
        yticks = 0:300:1800, yflip=false,
        xlabel="time", ylabel="freq",
        prompt = false,
    )
    sd = (; time=sj.time, freq=sj.freq, power=abs.(sj.power-sm.power))
    plot(fun(sj), fun(sm), fun(sd), p0), gui()
#   soundsc([x; ym; yj], S) # listen to compare


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

## numeric/object traits
# trait for objects that have an ordering
abstract TypeOrder
immutable HasOrder <: TypeOrder end
immutable Unordered <: TypeOrder end

(::Type{TypeOrder})(instance) = TypeOrder(typeof(instance))
(::Type{TypeOrder}){T<:Real}(::Type{T}) = HasOrder()
(::Type{TypeOrder}){T}(::Type{T}) = Unordered()

# trait for objects that support arithmetic
abstract TypeArithmetic
immutable ArithmeticRounds <: TypeArithmetic end     # least significant bits can be lost
immutable ArithmeticOverflows <: TypeArithmetic end  #  most significant bits can be lost
immutable ArithmeticUnknown <: TypeArithmetic end

(::Type{TypeArithmetic})(instance) = TypeArithmetic(typeof(instance))
(::Type{TypeArithmetic}){T<:AbstractFloat}(::Type{T}) = ArithmeticRounds()
(::Type{TypeArithmetic}){T<:Integer}(::Type{T}) = ArithmeticOverflows()
(::Type{TypeArithmetic}){T}(::Type{T}) = ArithmeticUnknown()


text:
using Base.Test

a_foo = 0

const foo_keymap = {
    'a' => :( global a_foo; a_foo += 1)
}

b_foo = 0

const foo2_keymap = {
    'b' => :( global b_foo; b_foo += 1)
}

a_bar = 0
b_bar = 0

const bar_keymap = {
    'a' => :( global a_bar; a_bar += 1),
    'b' => :( global b_bar; b_bar += 1)
}

@eval @Readline.keymap test1_func $foo_keymap

function run_test(f,buf)
    global a_foo, a_bar, b_bar
    a_foo = a_bar = b_bar = 0
    while !eof(buf)
        f(buf,nothing)
    end
end

run_test(test1_func,IOBuffer("aa"))
@test a_foo == 2

@eval @Readline.keymap test2_func $([foo2_keymap, foo_keymap])

run_test(test2_func,IOBuffer("aaabb"))
@test a_foo == 3
@test b_foo == 2

@eval @Readline.keymap test3_func $([bar_keymap, foo_keymap])

run_test(test3_func,IOBuffer("aab"))
@test a_bar == 2
@test b_bar == 1


text:
function obhod(r,side,path)
    count = 0
    count1 = 0
    f = true
    while isborder(r,side)==true
        if isborder(r,change_side(side))==true
            f = false
            break
        end
        move!(r,change_side(side))
        push!(path,rev(change_side(side)))
        count += 1
    end
    if f == true
        move!(r,side)
        push!(path,rev(side))
        while isborder(r,rev(change_side(side)))
            move!(r,side)
            push!(path,rev(side))
            count1 += 1
        end
        for i in 1:count
            move!(r,rev(change_side(side)))
            push!(path,change_side(side))
        end
    end
    return f, count1 - 1
end

function rev(side)
    if side == Nord
        return Sud
    end
    if side == Sud
        return Nord
    end
    if side == West
        return Ost
    end
    if side == Ost
        return West
    end
end

function move_to_wall_and_back(r, side)
    path = []
    count = 0
    while true
        if isborder(r, side) == false
            move!(r, side)
            push!(path, rev(side))
            count += 1
        else
            f, a = obhod(r, side, path)
            count += a
            if f == false
                break
            end
        end
    end
    way_back(r, path)
    return count
end

function move_to_wall(r, side)
    putmarker!(r)
    while isborder(r,side)
        move!(r, side)
        putmarker!(r)
    end
end

function move_count_and_back(r,side,count)
    path = []
    for i in 1:count
        if isborder(r, side) == false
            move!(r, side)
            putmarker!(r)
        else
            f, c = obhod(r, side, path)
            if f == false
                break
            end
            if i + c > count
                break
            end
        end
    end
    way_back(r, path)
end

function change_side(side)
    if side == Nord
        return West
    end
    if side == West
        return Sud
    end
    if side == Sud
        return Ost
    end
    if side == Ost
        return Nord
    end
end

function way_back(r,path)
    for i in 0:(length(path)-1)
        move!(r,path[length(path)-i])
    end
end

function corner(r)
    path = []
    while isborder(r, Sud) == false || isborder(r, Ost) == false
        if isborder(r, Sud) == false
            move!(r, Sud)
            push!(path, Nord)
        end
        if isborder(r, Ost) == false
            move!(r,Ost)
            push!(path, West)
        end
    end
    return path
end

function main(r)
    path = corner(r)
    count = move_to_wall_and_back(r, West)
    for i in 1:count
        move_count_and_back(r, West, count - i)
    end
    move_to_wall(r, Sud)
    way_back(r, path)
end


text:
using Test, Package

@test hello("Julia") == "Hello, Julia"
@test domath(2.0) ≈ 7.0
@test Package.compute_square(2.0) ≈ 4.0


text:
export discountFactor

struct ConstantBlackModel
    vol::Float64
    r::Float64
    q::Float64
end

struct TSBlackModel{S}
    surface::S
    r::Float64
    q::Float64
end

struct LocalVolatilityModel{S}
    surface::S
    r::Float64
    q::Float64
end

function discountFactor(model, t)
    exp(-model.r * t)
end

function logForward(model, lnspot, t)
    return lnspot + (model.r-model.q)*t
end

function forward(model, spot, t)
    return spot*exp((model.r-model.q)*t)
end

function varianceByLogmoneyness(model::ConstantBlackModel, y, t)
    return model.vol^2
end

function varianceByLogmoneyness(model, y, t)
    return varianceByLogmoneyness(model.surface, y, t)
end


text:
function open_pack!(
    pack_contents::Vector{Tuple{Int64,Int64}},
    pack_table,
    collection::Collection,
    decklist::DecklistArray,
    parameters::HearthstoneSimParameters,
    account_state::AccountState,
    sim_state::SimState,
    output::SimOutput
    )

    set = sim_state.set

    keep_going = true
    while keep_going
        StatsBase.sample!(pack_table[set][1], pack_table[set][2], pack_contents)

        if pack_okay(pack_contents, collection, account_state.pitytimer, set, pack_table)
            keep_going = false
        end
    end

    process_pack!(pack_contents, collection, decklist, account_state, set, parameters,sim_state, output)

    for r in 3:4
        if has_rarity(pack_contents, r)
            account_state.pitytimer[set][r] = 0
        else
            account_state.pitytimer[set][r] += 1
        end
    end
end

function process_pack!(pack_contents, collection, decklist, account_state, set, parameters, sim_state, output)
    for (r,c) in pack_contents
        collection_amount = collection[r][set][c] + 1
        necessary_amount = 0
        for deck in decklist[r][set]
            for (deck_card, deck_amt) in deck
                if (c == deck_card) && (deck_amt > necessary_amount)
                    necessary_amount = deck_amt
                end
            end
        end

        isgolden = rand() < parameters.golden_probs[r]

        if collection_amount > necessary_amount
            if isgolden
                account_state.dust += dust_golden[r]
            else
                account_state.dust += dust_regular[r]
            end
        else
            collection[r][set][c] += 1
            output.cards_from_packs[set][sim_state.rep, r] += 1
        end
        output.pack_card_rarity[sim_state.rep, r] += 1
    end
end

function pack_okay(pack_contents, collection, pity_timer, set, pack_table)
    ## at least one rare:
    has_minrarity(pack_contents, 2) || return false
    enforce_pitytimer!(pack_contents, pity_timer, set, pack_table)
    has_no_duplicate_leggos(pack_contents, collection, set) || return false
    has_no_playset_overflow(pack_contents) || return false

    return true

end

function has_minrarity(pack_contents, rarity)
    for (r,c) in pack_contents
        if r >= rarity
            return true
        end
    end
    return false
end

function has_rarity(pack_contents, rarity)
    for (r,c) in pack_contents
        if r == rarity
            return true
        end
    end
    return false
end

function has_no_pitytimer_violation(pack_contents, pity_timer, set)
    if pity_timer[set][4] >= 40 && !has_rarity(pack_contents, 4)
        return false
    end

    if pity_timer[set][3] >= 10 && !has_rarity(pack_contents, 3)
        return false
    end

    return true

end

function enforce_pitytimer!(pack_contents, pity_timer, set, pack_table)
    if pity_timer[set][4] >= 40 && !has_rarity(pack_contents, 4)
        for i in eachindex(pack_contents)
            r,c = pack_contents[i]
            if r >= 2
                pack_contents[i] = StatsBase.sample( filter(x -> x[1] == 4, pack_table[set][1]) )
                return nothing
            end
        end
    elseif pity_timer[set][3] >= 10 && !has_rarity(pack_contents,3)
        for i in eachindex(pack_contents)
            r,c = pack_contents[i]
            if r >= 2
                pack_contents[i] = StatsBase.sample( filter(x -> x[1] == 3, pack_table[set][1]) )
                return nothing
            end
        end
    end
end

function has_no_playset_overflow(pack_contents)
    for i in eachindex(pack_contents)
        count = 1
        for j in (i+1):length(pack_contents)
            if pack_contents[i] == pack_contents[j]
                count += 1
            end
        end
        if count > 2 || (count > 1 && pack_contents[i] == 4)
            return false
        end
    end

    return true
end

function has_no_duplicate_leggos(pack_contents, collection, set)
    all( collection[4][set] .>= 1 ) && return true

    for (r,c) in pack_contents
        if r == 4
            if collection[4][set][c] >= 1
                return false
            end
        end
    end
    return true
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
using DataDeps

register(DataDep(
    "zoo",
    """
    Dataset: zoo
    Author: Richard Forsyth
    Website: https://archive.ics.uci.edu/ml/datasets/zoo

    Observations: 101
    Features:     17 (16 numeric)
    Classes:      7

    Note: This version strips 'animal name' from attributes.
    """,
    "https://archive.ics.uci.edu/ml/machine-learning-databases/zoo/zoo.data",
    "cddc71c26ab9bc82795b8f4ff114cade41885d92720c6af29ffb69bcf73f0315",
    post_fetch_method = (path -> begin
        MVDA.process_dataset(path,
            header=false,
            target_index=18,
            feature_indices=2:17,
            ext=".csv",
        )
    end),
))


text:
using Test
using CorticalSpecrtralTemporalResponses
using AuditoryBistabilityLE
using SampledSignals
using TOML
using LinearAlgebra
using Statistics
using JLD2
using FileIO
using DataFrames
using Unitful

N = 50

pstream(p) = sum(p.counts[1][p.counts[2]]) / sum(p.counts[1])

@testset "Stremaing Bistability" begin
  @testset "Basic Streaming" begin
    params = load("object_test_params.jld2")["params"]

    # TODO: get this function running

    settings = TOML.parsefile("test_settings.toml")
    settings["stimulus"]["repeats"] = 10

    @test bistable_model(params, settings, interactive=true,
                   progressbar=false) != false
    # uncomment to profile
    # @time bistable_model(100, params, "test_settings.toml", interactive=true,
    #                      progressbar=true)

    # TODO: test with the always-fused stimulus from Elhilali et al 2009
    # (not at immediate concern, but should be present for the final
    # mod)
    params[:Δf] = 12
    settings["stimulus"]["repeats"] = 50
    df12 = bistable_model(params, "test_settings.toml",
                          interactive=true, progressbar=true)

    params[:Δf] = 3
    df3 = bistable_model(params, "test_settings.toml", interactive=true,
                         progressbar=true)

    @test pstream(df12.percepts) > 0.9
    @test pstream(df3.percepts) < 0.1
  end

  @testset "Track-level Bistability" begin

    params = Dict(
      :Δt         => 240ms, :Δf        => 3,
      :f          => 500Hz, :condition => :track,
      :τ_x        => 500ms, :c_x       => 3.0,
      :f_c_a => 0, :f_c_m => 0, :f_c_σ => 0,
      :s_c_a => 0, :s_c_m => 0, :s_c_σ => 0,
      :t_W_m_σ      => 15.0, #5.0
      :t_W_m_σ_t    => 7.0,   :t_W_m_σ_ϕ   => 7.0,
      :t_W_m_σ_N    => 3.0,   :t_W_m_c     => 6.0,
      :t_τ_m        => 350ms, :t_c_m       => 100,
      :t_τ_a        => 3s,    :t_c_a       => 6,
      :t_τ_σ        => 500ms, :t_c_σ       => 0.2
     )

    params[:Δf] = 12
    df12 = bistable_model(100, params, "test_settings.toml", interactive=true,
                          progressbar=true)

    params[:Δf] = 6
    df6 = bistable_model(100, params, "test_settings.toml", interactive=true,
                         progressbar=true)

    params[:Δf] = 3
    df3 = bistable_model(100, params, "test_settings.toml", interactive=true,
                         progressbar=true)

    # bare minimum tests that should pass if there's stimulus selective
    # bistability the results are somewhat random, so not all tests will pass
    # on all runs, but the should normally *mostly* pass.
    # once the parameters have been fit well, the goal will be to select
    # parameter values that consistently pass
    @test pstream(df12.percepts) > 0.8
    @test pstream(df3.percepts) < 0.2
    @test pstream(df12.percepts) > pstream(df6.percepts)
    @test pstream(df6.percepts) > pstream(df3.percepts)
  end

  @testset "Scale-level Bistability" begin

    params = Dict(
      :Δt         => 240ms, :Δf        => 3,
      :f          => 500Hz, :condition => :track,
      :τ_x        => 500ms, :c_x       => 3.0,
      :f_c_a => 0, :f_c_m => 0, :f_c_σ => 0,
      :s_W_m_σ      => 15.0,
      :s_W_m_c      => 6.0,
      :s_τ_m        => 350ms, :s_c_m       => 0,
      :s_τ_a        => 3s,    :s_c_a       => 25,
      :s_τ_σ        => 500ms, :s_c_σ       => 0.2,
      :t_c_a => 0, :t_c_m => 0, :t_c_σ => 0
     )

    params[:Δf] = 12
    df12 = bistable_model(100, params, "test_settings.toml", interactive=true,
                          progressbar=true)

    params[:Δf] = 6
    df6 = bistable_model(100, params, "test_settings.toml", interactive=true,
                         progressbar=true)

    params[:Δf] = 3
    df3 = bistable_model(100, params, "test_settings.toml", interactive=true,
                         progressbar=true)

    # bare minimum tests that should pass if there's stimulus
    # selective bistability
    @test pstream(df12.percepts) > 0.8
    @test pstream(df3.percepts) < 0.2
    @test pstream(df12.percepts) > pstream(df6.percepts)
    @test pstream(df6.percepts) > pstream(df3.percepts)
  end

  @testset "Freq-level Bistability" begin
    params = Dict(
      :Δt         => 240ms, :Δf        => 3,
      :f          => 500Hz, :condition => :track,
      :τ_x        => 500ms, :c_x       => 3.0,
      :f_W_m_σ      => 5.6,
      :f_W_m_c      => 6.0,
      :f_τ_m        => 350ms, :f_c_m       => 32,
      :f_τ_a        => 3s,    :f_c_a       => 5,
      :f_τ_σ        => 500ms, :f_c_σ       => 0.2,
      :s_c_a => 0, :s_c_m => 0, :s_c_σ => 0,
      :t_c_a => 0, :t_c_m => 0, :t_c_σ => 0
     )

    params[:Δf] = 12
    df12 = bistable_model(50, params, "test_settings.toml", interactive=true,
                          progressbar=true)

    params[:Δf] = 6
    df6 = bistable_model(50, params, "test_settings.toml", interactive=true,
                         progressbar=true)

    params[:Δf] = 3
    df3 = bistable_model(50, params, "test_settings.toml", interactive=true,
                          progressbar=true)

    # bare minimum tests that should pass if there's stimulus
    # selective bistability
    @test pstream(df12.percepts) > 0.8
    @test pstream(df3.percepts) < 0.2
    @test pstream(df12.percepts) > pstream(df6.percepts)
    @test pstream(df6.percepts) > pstream(df3.percepts)
  end
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
"""
euler = [ψ, θ, ϕ]  # yaw, pitch, roll
"""
function distribution_info(manoeuvre::Symbol)
    p_min, p_max = nothing, nothing
    v_min, v_max = nothing, nothing
    euler_min, euler_max = nothing, nothing
    ω_min, ω_max = nothing, nothing
    if manoeuvre == :hovering
        p_min, p_max = [-1, -1, -11.0], [1, 1, -9.0]
        v_min, v_max = [-1, -1, -1.0], [1, 1, 1.0]
        euler_min, euler_max = [deg2rad(-5), deg2rad(-5), deg2rad(-5)], [deg2rad(5), deg2rad(5), deg2rad(5)]
        ω_min, ω_max = [deg2rad(-5), deg2rad(-5), deg2rad(-5)], [deg2rad(5), deg2rad(5), deg2rad(5)]
    elseif manoeuvre == :forward
        p_min, p_max = [-1, -1, -11.0], [1, 1, -9.0]
        v_min, v_max = [3, -1, -1.0], [7, 1, 1.0]
        euler_min, euler_max = [deg2rad(-2), deg2rad(-10), deg2rad(-2)], [deg2rad(2), deg2rad(-5), deg2rad(2)]
        ω_min, ω_max = [deg2rad(-2), deg2rad(-2), deg2rad(-2)], [deg2rad(2), deg2rad(2), deg2rad(2)]
    elseif manoeuvre == :debug
        p_min, p_max = [0, 0, -10.0], [0, 0, -10.0]
        v_min, v_max = zeros(3), zeros(3)
        euler_min, euler_max = zeros(3), zeros(3)
        ω_min, ω_max = zeros(3), zeros(3)
    else
        error("Invalid manoeuvre")
    end
    min_nt = (;
              p=p_min,
              v=v_min,
              euler=euler_min,
              ω=ω_min,
             )
    max_nt = (;
              p=p_max,
              v=v_max,
              euler=euler_max,
              ω=ω_max,
             )
    min_nt, max_nt
end

"""
    run_multiple_sim()

A function for running multiple simulation.
This is used for 2nd-year report.

# Notes
- collector = collect (sequential computing)
- collector = Transducers.tcollect (parallel computing)
- manoeuvre = :hovering or :forward (:debug for debugging)
- `manoeuvres` is an array of manoeuvres.
"""
function run_multiple_sim(
        N::Int,
        manoeuvres::AbstractArray{Symbol},
        methods::AbstractArray{Symbol},
        _faults,
        θs_array,  # path parameter
        t0::Real,
        tf::Real,
        h_threshold::Union{Real, Nothing},  # m (nothing: no constraint)
        actual_time_limit::Union{Real, Nothing},  # s
        N_thread::Int,
        collector,
        will_plot::Bool,
        seed::Int,
    )
    println("# of simulation cases: $(N)")
    println("manoeuvres: $(manoeuvres)")
    println("methods: $(methods)")
    if collector == tcollect
        println("Parallel computing...")
        # GR has concurrency issues
        will_plot == true ? error("plotting figures not supported in tcollect") : nothing
    elseif collector == collect
        println("Sequential computing...")
    else
        error("Invalid collector")
    end
    Random.seed!(seed)
    _dir_log = joinpath("data", "N_$(N)_seed_$(seed)")
    _dir_figure = joinpath("figures", "N_$(N)_seed_$(seed)")
    multicopter = LeeHexacopter()  # dummy
    trajs_des = θs_array |> Map(θs -> Bezier(θs, t0, tf)) |> collect
    # run sim and save fig
    for manoeuvre in manoeuvres
        x0s = 1:N |> Map(i -> FTCTests.sample(multicopter, distribution_info(manoeuvre)...)) |> collect
        for method in methods
            dir_log = joinpath(joinpath(_dir_log, String(manoeuvre)), String(method))
            dir_figure = will_plot ? joinpath(joinpath(_dir_figure, String(manoeuvre)), String(method)) : nothing
            if method == :adaptive  # method `:adaptive` does not require FDI information; not affected by FDI delay time constant `τ`.
                τs = [0.0]
            elseif method == :adaptive2optim
                τs = [0.0, 0.1]
            else
                error("Invalid method")
            end
            for τ in τs
                case_numbers_partition = 1:N |> Partition(N_thread; flush=true) |> Map(copy) |> collect
                for case_numbers in case_numbers_partition
                    @time _ = zip(case_numbers,
                                x0s[case_numbers],
                                _faults[case_numbers],
                                trajs_des[case_numbers],
                               ) |>
                    MapSplat((i, x0, _fault, traj_des) -> FTCTests.run_sim(method, x0, multicopter,
                                                                           FaultSet(_fault...),
                                                                           DelayFDI(τ),
                                                                           traj_des,
                                                                           dir_log,
                                                                           i;
                                                                           dir_figure=dir_figure,
                                                                           t0=t0, tf=tf,
                                                                           h_threshold=h_threshold,
                                                                           actual_time_limit=actual_time_limit,
                                                                          )) |> collector
                end
            end
        end
    end
    nothing
end


text:
function nbytes_to_unsigned_type(n)
    if     n == 1  UInt8
    elseif n == 2  UInt16
    elseif n == 4  UInt32
    elseif n == 8  UInt64
    elseif n == 16 UInt128
    else "Length of n must be a power of two and maximum 16." |>
        ArgumentError |> throw
    end
end

# TODO: little endian only
# This one gets optimized out completely in 0.7+ :-)
"""
    BTCParser.to_unsigned(x::NTuple{N, UInt8})::Unsigned

Convert a tuple to an unsigned (`UInt8`, `UInt16`, `UInt32`, `UInt64`,
`UInt128`, `UInt256`) with the identical order of bits.

Results depend on the endianness of your CPU architecture.
"""
@generated function to_unsigned(x::NTuple{N, UInt8}) where N

    T = nbytes_to_unsigned_type(N)

    sb = ntuple(i -> T(1) << (8 * (i - 1)), N)

    ex = :($T(0))
    for i in 1:N
        ex = :($ex + $sb[$i] * x[$i])
    end

    return ex
end

# TODO: figure out a no-op way to do this
function to_unsigned(x::NTuple{32, UInt8})
    x |> Ref |>
        x -> Base.unsafe_convert(Ptr{NTuple{32, UInt8}}, x) |>
        x -> convert(Ptr{UInt256}, x) |>
        x -> unsafe_load(x)
end

# TODO: figure out a no-op way to do this
function to_byte_tuple(x::T) where T <: Union{UInt256, UInt128}
    x |>
        Ref |>
        x -> Base.unsafe_convert(Ptr{T}, x) |>
        x -> convert(Ptr{NTuple{sizeof(T), UInt8}}, x) |>
        x -> unsafe_load(x)
end

# TODO: Little endian only
# TODO: figure out a no-op way to do this
# NOTE: This is really slow for UInt128
"""
    BTCParser.to_byte_tuple(x::Unsigned)::NTuple{N, UInt8}

Convert an `Unsigned` (`UInt8`, `UInt16`, `UInt32`, `UInt64`, `UInt128`,
`UInt256`) into an `NTuple{N, UInt8}` with the identical order of bits.

Results depend on the endianness of your CPU architecture.
"""
function to_byte_tuple(x::T) where T <: Unsigned
    ntuple(i -> (x >> (8 * (i - 1))) % UInt8, sizeof(T))
end


text:
# A original flavor of the Conjugate Gradient Algorithm.
# * It performs the Conjugate Gradient and record all the residuals vectors for
#   computing the lanczos iterations.

mutable struct ConjGrad
    r::Vector{Float64}
    rnew::Vector{Float64}
    d::Vector{Float64}
    A::Function
    x::Vector{Float64}
    b::Vector{Float64}
    Ad::Vector{Float64}
    itr::Int64

    function ConjGrad(A::Function, b, x0=nothing)
        this = new()
        this.A = A
        this.x = x0 === nothing ? b .+ 0.1  : x0  # just to handle matrix A that has eigenvalue of exactly 1.
        this.r = b - A(this.x)
        this.rnew = similar(this.r)
        this.d = this.r
        this.itr = 0
        this.Ad = similar(this.d)

        return this
    end

    function ConjGrad(A::AbstractMatrix{Float64}, b::VecOrMat{Float64})
        return ConjGrad((x)->A*x, b)
    end

end

function (this::ConjGrad)()
    r = this.r
    if r == 0
        return 0 # The problem is solved already.
    end
    A = this.A
    d = this.d
    Ad = this.Ad
    Ad .= A(d)

    a = dot(r, r)/dot(d, Ad)
    if a < 0
        error("CG got a non-definite matrix")
    end

    this.x += a*d
    this.rnew .= r - a*Ad                # update rnew
    b = dot(this.rnew, this.rnew)/dot(r, r)
    # @assert abs(dot(rnew + β*d, Ad)) < 1e-8 "Not conjugate"
    @. this.d = this.rnew + b*d
    this.r = copy(this.rnew)                      # Override
    this.itr += 1
    return convert(Float64, norm(this.rnew))
end

function GetCurrentResidualNorm(this::ConjGrad)
    return norm(this.r)
end


text:
module MRIFiles

using Graphics: @mustimplement
using Reexport

using MRIBase
using Reexport
using FileIO
using LinearAlgebra
using HDF5
using LightXML
using NIfTI

abstract type MRIFile end
@mustimplement MRIBase.RawAcquisitionData(::MRIFile)

include("ISMRMRD/ISMRMRD.jl")
include("Bruker/Bruker.jl")
include("Siemens/Siemens.jl")
include("Philips/Philips.jl")
include("GE/GE.jl")
include("Others/Others.jl")

#include("Nifti.jl")

end # module


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "fiveDimSys_model.jl")
include(@relpath "fiveDimSys_specifications.jl")

S = fiveDimSys_model()
X0, options = fiveDimSys_specification()
options = Options(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# partition with 1D blocks
partition = [[i] for i in 1:5]

# reachability algorithms
algorithm_reach1 = BFFPSV18(:δ => 0.005, :vars => [1, 2],
                            :partition => partition)
algorithm_reach2 = GLGM06(:δ => 0.005, :max_order => 5)

# compute flowpipe
options[:mode] = "reach"
solution1 = solve(problem, options; op=algorithm_reach1)
solution2 = solve(problem, options; op=algorithm_reach2)

# project flowpipe to time and x₁
solution1.options[:plot_vars] = [1, 2]
solution_proj1 = project(solution1)
solution2.options[:plot_vars] = [1, 2]
solution_proj2 = project(solution2)

# plot projection
plot(solution_proj1)
savefig("fiveDimSys1")
plot(solution_proj2)
savefig("fiveDimSys2")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks