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

module Sort

import ..@__MODULE__, ..parentmodule
const Base = parentmodule(@__MODULE__)
using .Base.Order
using .Base: copymutable, LinearIndices, length, (:),
    eachindex, axes, first, last, similar, zip, OrdinalRange,
    AbstractVector, @inbounds, AbstractRange, @eval, @inline, Vector, @noinline,
    AbstractMatrix, AbstractUnitRange, isless, identity, eltype, >, <, <=, >=, |, +, -, *, !,
    extrema, sub_with_overflow, add_with_overflow, oneunit, div, getindex, setindex!,
    length, resize!, fill, Missing, require_one_based_indexing, keytype

using .Base: >>>, !==

import .Base:
    sort,
    sort!,
    issorted,
    sortperm,
    to_indices

export # also exported by Base
    # order-only:
    issorted,
    searchsorted,
    searchsortedfirst,
    searchsortedlast,
    # order & algorithm:
    sort,
    sort!,
    sortperm,
    sortperm!,
    partialsort,
    partialsort!,
    partialsortperm,
    partialsortperm!,
    # algorithms:
    InsertionSort,
    QuickSort,
    MergeSort,
    PartialQuickSort

export # not exported by Base
    Algorithm,
    DEFAULT_UNSTABLE,
    DEFAULT_STABLE,
    SMALL_ALGORITHM,
    SMALL_THRESHOLD

## functions requiring only ordering ##

function issorted(itr, order::Ordering)
    y = iterate(itr)
    y === nothing && return true
    prev, state = y
    y = iterate(itr, state)
    while y !== nothing
        this, state = y
        lt(order, this, prev) && return false
        prev = this
        y = iterate(itr, state)
    end
    return true
end

"""
    issorted(v, lt=isless, by=identity, rev:Bool=false, order::Ordering=Forward)

Test whether a vector is in sorted order. The `lt`, `by` and `rev` keywords modify what
order is considered to be sorted just as they do for [`sort`](@ref).

# Examples
```jldoctest
julia> issorted([1, 2, 3])
true

julia> issorted([(1, "b"), (2, "a")], by = x -> x[1])
true

julia> issorted([(1, "b"), (2, "a")], by = x -> x[2])
false

julia> issorted([(1, "b"), (2, "a")], by = x -> x[2], rev=true)
true
```
"""
issorted(itr;
    lt=isless, by=identity, rev::Union{Bool,Nothing}=nothing, order::Ordering=Forward) =
    issorted(itr, ord(lt,by,rev,order))

function partialsort!(v::AbstractVector, k::Union{Integer,OrdinalRange}, o::Ordering)
    inds = axes(v, 1)
    sort!(v, first(inds), last(inds), PartialQuickSort(k), o)
    maybeview(v, k)
end

maybeview(v, k) = view(v, k)
maybeview(v, k::Integer) = v[k]

"""
    partialsort!(v, k; by=<transform>, lt=<comparison>, rev=false)

Partially sort the vector `v` in place, according to the order specified by `by`, `lt` and
`rev` so that the value at index `k` (or range of adjacent values if `k` is a range) occurs
at the position where it would appear if the array were fully sorted via a non-stable
algorithm. If `k` is a single index, that value is returned; if `k` is a range, an array of
values at those indices is returned. Note that `partialsort!` does not fully sort the input
array.

# Examples
```jldoctest
julia> a = [1, 2, 4, 3, 4]
5-element Array{Int64,1}:
 1
 2
 4
 3
 4

julia> partialsort!(a, 4)
4

julia> a
5-element Array{Int64,1}:
 1
 2
 3
 4
 4

julia> a = [1, 2, 4, 3, 4]
5-element Array{Int64,1}:
 1
 2
 4
 3
 4

julia> partialsort!(a, 4, rev=true)
2

julia> a
5-element Array{Int64,1}:
 4
 4
 3
 2
 1
```
"""
partialsort!(v::AbstractVector, k::Union{Integer,OrdinalRange};
             lt=isless, by=identity, rev::Union{Bool,Nothing}=nothing, order::Ordering=Forward) =
    partialsort!(v, k, ord(lt,by,rev,order))

"""
    partialsort(v, k, by=<transform>, lt=<comparison>, rev=false)

Variant of [`partialsort!`](@ref) which copies `v` before partially sorting it, thereby returning the
same thing as `partialsort!` but leaving `v` unmodified.
"""
partialsort(v::AbstractVector, k::Union{Integer,OrdinalRange}; kws...) =
    partialsort!(copymutable(v), k; kws...)

# This implementation of `midpoint` is performance-optimized but safe
# only if `lo <= hi`.
midpoint(lo::T, hi::T) where T<:Integer = lo + ((hi - lo) >>> 0x01)
midpoint(lo::Integer, hi::Integer) = midpoint(promote(lo, hi)...)

# reference on sorted binary search:
#   http://www.tbray.org/ongoing/When/200x/2003/03/22/Binary

# index of the first value of vector a that is greater than or equal to x;
# returns length(v)+1 if x is greater than all values in v.
function searchsortedfirst(v::AbstractVector, x, lo::T, hi::T, o::Ordering)::keytype(v) where T<:Integer
    u = T(1)
    lo = lo - u
    hi = hi + u
    @inbounds while lo < hi - u
        m = midpoint(lo, hi)
        if lt(o, v[m], x)
            lo = m
        else
            hi = m
        end
    end
    return hi
end

# index of the last value of vector a that is less than or equal to x;
# returns 0 if x is less than all values of v.
function searchsortedlast(v::AbstractVector, x, lo::T, hi::T, o::Ordering)::keytype(v) where T<:Integer
    u = T(1)
    lo = lo - u
    hi = hi + u
    @inbounds while lo < hi - u
        m = midpoint(lo, hi)
        if lt(o, x, v[m])
            hi = m
        else
            lo = m
        end
    end
    return lo
end

# returns the range of indices of v equal to x
# if v does not contain x, returns a 0-length range
# indicating the insertion point of x
function searchsorted(v::AbstractVector, x, ilo::T, ihi::T, o::Ordering)::UnitRange{keytype(v)} where T<:Integer
    u = T(1)
    lo = ilo - u
    hi = ihi + u
    @inbounds while lo < hi - u
        m = midpoint(lo, hi)
        if lt(o, v[m], x)
            lo = m
        elseif lt(o, x, v[m])
            hi = m
        else
            a = searchsortedfirst(v, x, max(lo,ilo), m, o)
            b = searchsortedlast(v, x, m, min(hi,ihi), o)
            return a : b
        end
    end
    return (lo + 1) : (hi - 1)
end

function searchsortedlast(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering)::keytype(a)
    require_one_based_indexing(a)
    if step(a) == 0
        lt(o, x, first(a)) ? 0 : length(a)
    else
        n = round(Integer, clamp((x - first(a)) / step(a) + 1, 1, length(a)))
        lt(o, x, a[n]) ? n - 1 : n
    end
end

function searchsortedfirst(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering)::keytype(a)
    require_one_based_indexing(a)
    if step(a) == 0
        lt(o, first(a), x) ? length(a) + 1 : 1
    else
        n = round(Integer, clamp((x - first(a)) / step(a) + 1, 1, length(a)))
        lt(o, a[n] ,x) ? n + 1 : n
    end
end

function searchsortedlast(a::AbstractRange{<:Integer}, x::Real, o::DirectOrdering)::keytype(a)
    require_one_based_indexing(a)
    h = step(a)
    if h == 0
        lt(o, x, first(a)) ? 0 : length(a)
    elseif h > 0 && x < first(a)
        firstindex(a) - 1
    elseif h > 0 && x >= last(a)
        lastindex(a)
    elseif h < 0 && x > first(a)
        firstindex(a) - 1
    elseif h < 0 && x <= last(a)
        lastindex(a)
    else
        if o isa ForwardOrdering
            fld(floor(Integer, x) - first(a), h) + 1
        else
            fld(ceil(Integer, x) - first(a), h) + 1
        end
    end
end

function searchsortedfirst(a::AbstractRange{<:Integer}, x::Real, o::DirectOrdering)::keytype(a)
    require_one_based_indexing(a)
    h = step(a)
    if h == 0
        lt(o, first(a), x) ? length(a)+1 : 1
    elseif h > 0 && x <= first(a)
        firstindex(a)
    elseif h > 0 && x > last(a)
        lastindex(a) + 1
    elseif h < 0 && x >= first(a)
        firstindex(a)
    elseif h < 0 && x < last(a)
        lastindex(a) + 1
    else
        if o isa ForwardOrdering
            -fld(floor(Integer, -x) + Signed(first(a)), h) + 1
        else
            -fld(ceil(Integer, -x) + Signed(first(a)), h) + 1
        end
    end
end

function searchsortedfirst(a::AbstractRange{<:Integer}, x::Unsigned, o::DirectOrdering)::keytype(a)
    require_one_based_indexing(a)
    if lt(o, first(a), x)
        if step(a) == 0
            length(a) + 1
        else
            min(cld(x - first(a), step(a)), length(a)) + 1
        end
    else
        1
    end
end

function searchsortedlast(a::AbstractRange{<:Integer}, x::Unsigned, o::DirectOrdering)::keytype(a)
    require_one_based_indexing(a)
    if lt(o, x, first(a))
        0
    elseif step(a) == 0
        length(a)
    else
        min(fld(x - first(a), step(a)) + 1, length(a))
    end
end

searchsorted(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering) =
    searchsortedfirst(a, x, o) : searchsortedlast(a, x, o)

for s in [:searchsortedfirst, :searchsortedlast, :searchsorted]
    @eval begin
        $s(v::AbstractVector, x, o::Ordering) = (inds = axes(v, 1); $s(v,x,first(inds),last(inds),o))
        $s(v::AbstractVector, x;
           lt=isless, by=identity, rev::Union{Bool,Nothing}=nothing, order::Ordering=Forward) =
            $s(v,x,ord(lt,by,rev,order))
    end
end

"""
    searchsorted(a, x; by=<transform>, lt=<comparison>, rev=false)

Return the range of indices of `a` which compare as equal to `x` (using binary search)
according to the order specified by the `by`, `lt` and `rev` keywords, assuming that `a`
is already sorted in that order. Return an empty range located at the insertion point
if `a` does not contain values equal to `x`.

# Examples
```jldoctest
julia> searchsorted([1, 2, 4, 5, 5, 7], 4) # single match
3:3

julia> searchsorted([1, 2, 4, 5, 5, 7], 5) # multiple matches
4:5

julia> searchsorted([1, 2, 4, 5, 5, 7], 3) # no match, insert in the middle
3:2

julia> searchsorted([1, 2, 4, 5, 5, 7], 9) # no match, insert at end
7:6

julia> searchsorted([1, 2, 4, 5, 5, 7], 0) # no match, insert at start
1:0
```
""" searchsorted

"""
    searchsortedfirst(a, x; by=<transform>, lt=<comparison>, rev=false)

Return the index of the first value in `a` greater than or equal to `x`, according to the
specified order. Return `length(a) + 1` if `x` is greater than all values in `a`.
`a` is assumed to be sorted.

# Examples
```jldoctest
julia> searchsortedfirst([1, 2, 4, 5, 5, 7], 4) # single match
3

julia> searchsortedfirst([1, 2, 4, 5, 5, 7], 5) # multiple matches
4

julia> searchsortedfirst([1, 2, 4, 5, 5, 7], 3) # no match, insert in the middle
3

julia> searchsortedfirst([1, 2, 4, 5, 5, 7], 9) # no match, insert at end
7

julia> searchsortedfirst([1, 2, 4, 5, 5, 7], 0) # no match, insert at start
1
```
""" searchsortedfirst

"""
    searchsortedlast(a, x; by=<transform>, lt=<comparison>, rev=false)

Return the index of the last value in `a` less than or equal to `x`, according to the
specified order. Return `0` if `x` is less than all values in `a`. `a` is assumed to
be sorted.

# Examples
```jldoctest
julia> searchsortedlast([1, 2, 4, 5, 5, 7], 4) # single match
3

julia> searchsortedlast([1, 2, 4, 5, 5, 7], 5) # multiple matches
5

julia> searchsortedlast([1, 2, 4, 5, 5, 7], 3) # no match, insert in the middle
2

julia> searchsortedlast([1, 2, 4, 5, 5, 7], 9) # no match, insert at end
6

julia> searchsortedlast([1, 2, 4, 5, 5, 7], 0) # no match, insert at start
0
```
""" searchsortedlast

## sorting algorithms ##

abstract type Algorithm end

struct InsertionSortAlg <: Algorithm end
struct QuickSortAlg     <: Algorithm end
struct MergeSortAlg     <: Algorithm end

"""
    PartialQuickSort{T <: Union{Integer,OrdinalRange}}

Indicate that a sorting function should use the partial quick sort
algorithm. Partial quick sort returns the smallest `k` elements sorted from smallest
to largest, finding them and sorting them using [`QuickSort`](@ref).

Characteristics:
  * *not stable*: does not preserve the ordering of elements which
    compare equal (e.g. "a" and "A" in a sort of letters which
    ignores case).
  * *in-place* in memory.
  * *divide-and-conquer*: sort strategy similar to [`MergeSort`](@ref).
"""
struct PartialQuickSort{T <: Union{Integer,OrdinalRange}} <: Algorithm
    k::T
end

"""
    InsertionSort

Indicate that a sorting function should use the insertion sort
algorithm. Insertion sort traverses the collection one element
at a time, inserting each element into its correct, sorted position in
the output list.

Characteristics:
  * *stable*: preserves the ordering of elements which
    compare equal (e.g. "a" and "A" in a sort of letters
    which ignores case).
  * *in-place* in memory.
  * *quadratic performance* in the number of elements to be sorted:
    it is well-suited to small collections but should not be used for large ones.
"""
const InsertionSort = InsertionSortAlg()
"""
    QuickSort

Indicate that a sorting function should use the quick sort
algorithm, which is *not* stable.

Characteristics:
  * *not stable*: does not preserve the ordering of elements which
    compare equal (e.g. "a" and "A" in a sort of letters which
    ignores case).
  * *in-place* in memory.
  * *divide-and-conquer*: sort strategy similar to [`MergeSort`](@ref).
  * *good performance* for large collections.
"""
const QuickSort     = QuickSortAlg()
"""
    MergeSort

Indicate that a sorting function should use the merge sort
algorithm. Merge sort divides the collection into
subcollections and repeatedly merges them, sorting each
subcollection at each step, until the entire
collection has been recombined in sorted form.

Characteristics:
  * *stable*: preserves the ordering of elements which compare
    equal (e.g. "a" and "A" in a sort of letters which ignores
    case).
  * *not in-place* in memory.
  * *divide-and-conquer* sort strategy.
"""
const MergeSort     = MergeSortAlg()

const DEFAULT_UNSTABLE = QuickSort
const DEFAULT_STABLE   = MergeSort
const SMALL_ALGORITHM  = InsertionSort
const SMALL_THRESHOLD  = 20

function sort!(v::AbstractVector, lo::Integer, hi::Integer, ::InsertionSortAlg, o::Ordering)
    @inbounds for i = lo+1:hi
        j = i
        x = v[i]
        while j > lo
            if lt(o, x, v[j-1])
                v[j] = v[j-1]
                j -= 1
                continue
            end
            break
        end
        v[j] = x
    end
    return v
end

# selectpivot!
#
# Given 3 locations in an array (lo, mi, and hi), sort v[lo], v[mi], v[hi]) and
# choose the middle value as a pivot
#
# Upon return, the pivot is in v[lo], and v[hi] is guaranteed to be
# greater than the pivot

@inline function selectpivot!(v::AbstractVector, lo::Integer, hi::Integer, o::Ordering)
    @inbounds begin
        mi = midpoint(lo, hi)

        # sort v[mi] <= v[lo] <= v[hi] such that the pivot is immediately in place
        if lt(o, v[lo], v[mi])
            v[mi], v[lo] = v[lo], v[mi]
        end

        if lt(o, v[hi], v[lo])
            if lt(o, v[hi], v[mi])
                v[hi], v[lo], v[mi] = v[lo], v[mi], v[hi]
            else
                v[hi], v[lo] = v[lo], v[hi]
            end
        end

        # return the pivot
        return v[lo]
    end
end

# partition!
#
# select a pivot, and partition v according to the pivot

function partition!(v::AbstractVector, lo::Integer, hi::Integer, o::Ordering)
    pivot = selectpivot!(v, lo, hi, o)
    # pivot == v[lo], v[hi] > pivot
    i, j = lo, hi
    @inbounds while true
        i += 1; j -= 1
        while lt(o, v[i], pivot); i += 1; end;
        while lt(o, pivot, v[j]); j -= 1; end;
        i >= j && break
        v[i], v[j] = v[j], v[i]
    end
    v[j], v[lo] = pivot, v[j]

    # v[j] == pivot
    # v[k] >= pivot for k > j
    # v[i] <= pivot for i < j
    return j
end

function sort!(v::AbstractVector, lo::Integer, hi::Integer, a::QuickSortAlg, o::Ordering)
    @inbounds while lo < hi
        hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)
        j = partition!(v, lo, hi, o)
        if j-lo < hi-j
            # recurse on the smaller chunk
            # this is necessary to preserve O(log(n))
            # stack space in the worst case (rather than O(n))
            lo < (j-1) && sort!(v, lo, j-1, a, o)
            lo = j+1
        else
            j+1 < hi && sort!(v, j+1, hi, a, o)
            hi = j-1
        end
    end
    return v
end

function sort!(v::AbstractVector, lo::Integer, hi::Integer, a::MergeSortAlg, o::Ordering, t=similar(v,0))
    @inbounds if lo < hi
        hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)

        m = midpoint(lo, hi)
        (length(t) < m-lo+1) && resize!(t, m-lo+1)

        sort!(v, lo,  m,  a, o, t)
        sort!(v, m+1, hi, a, o, t)

        i, j = 1, lo
        while j <= m
            t[i] = v[j]
            i += 1
            j += 1
        end

        i, k = 1, lo
        while k < j <= hi
            if lt(o, v[j], t[i])
                v[k] = v[j]
                j += 1
            else
                v[k] = t[i]
                i += 1
            end
            k += 1
        end
        while k < j
            v[k] = t[i]
            k += 1
            i += 1
        end
    end

    return v
end

function sort!(v::AbstractVector, lo::Integer, hi::Integer, a::PartialQuickSort{<:Integer},
               o::Ordering)
    @inbounds while lo < hi
        hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)
        j = partition!(v, lo, hi, o)
        if j >= a.k
            # we don't need to sort anything bigger than j
            hi = j-1
        elseif j-lo < hi-j
            # recurse on the smaller chunk
            # this is necessary to preserve O(log(n))
            # stack space in the worst case (rather than O(n))
            lo < (j-1) && sort!(v, lo, j-1, a, o)
            lo = j+1
        else
            (j+1) < hi && sort!(v, j+1, hi, a, o)
            hi = j-1
        end
    end
    return v
end

function sort!(v::AbstractVector, lo::Integer, hi::Integer, a::PartialQuickSort{T},
               o::Ordering) where T<:OrdinalRange
    @inbounds while lo < hi
        hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)
        j = partition!(v, lo, hi, o)

        if j <= first(a.k)
            lo = j+1
        elseif j >= last(a.k)
            hi = j-1
        else
            if j-lo < hi-j
                lo < (j-1) && sort!(v, lo, j-1, a, o)
                lo = j+1
            else
                hi > (j+1) && sort!(v, j+1, hi, a, o)
                hi = j-1
            end
        end
    end
    return v
end

## generic sorting methods ##

defalg(v::AbstractArray) = DEFAULT_STABLE
defalg(v::AbstractArray{<:Union{Number, Missing}}) = DEFAULT_UNSTABLE

function sort!(v::AbstractVector, alg::Algorithm, order::Ordering)
    inds = axes(v,1)
    sort!(v,first(inds),last(inds),alg,order)
end

"""
    sort!(v; alg::Algorithm=defalg(v), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)

Sort the vector `v` in place. [`QuickSort`](@ref) is used by default for numeric arrays while
[`MergeSort`](@ref) is used for other arrays. You can specify an algorithm to use via the `alg`
keyword (see [Sorting Algorithms](@ref) for available algorithms). The `by` keyword lets you provide
a function that will be applied to each element before comparison; the `lt` keyword allows
providing a custom "less than" function; use `rev=true` to reverse the sorting order. These
options are independent and can be used together in all possible combinations: if both `by`
and `lt` are specified, the `lt` function is applied to the result of the `by` function;
`rev=true` reverses whatever ordering specified via the `by` and `lt` keywords.

# Examples
```jldoctest
julia> v = [3, 1, 2]; sort!(v); v
3-element Array{Int64,1}:
 1
 2
 3

julia> v = [3, 1, 2]; sort!(v, rev = true); v
3-element Array{Int64,1}:
 3
 2
 1

julia> v = [(1, "c"), (3, "a"), (2, "b")]; sort!(v, by = x -> x[1]); v
3-element Array{Tuple{Int64,String},1}:
 (1, "c")
 (2, "b")
 (3, "a")

julia> v = [(1, "c"), (3, "a"), (2, "b")]; sort!(v, by = x -> x[2]); v
3-element Array{Tuple{Int64,String},1}:
 (3, "a")
 (2, "b")
 (1, "c")
```
"""
function sort!(v::AbstractVector;
               alg::Algorithm=defalg(v),
               lt=isless,
               by=identity,
               rev::Union{Bool,Nothing}=nothing,
               order::Ordering=Forward)
    ordr = ord(lt,by,rev,order)
    if (ordr === Forward || ordr === Reverse) && eltype(v)<:Integer
        n = length(v)
        if n > 1
            min, max = extrema(v)
            (diff, o1) = sub_with_overflow(max, min)
            (rangelen, o2) = add_with_overflow(diff, oneunit(diff))
            if !o1 && !o2 && rangelen < div(n,2)
                return sort_int_range!(v, rangelen, min, ordr === Reverse ? reverse : identity)
            end
        end
    end
    sort!(v, alg, ordr)
end

# sort! for vectors of few unique integers
function sort_int_range!(x::AbstractVector{<:Integer}, rangelen, minval, maybereverse)
    offs = 1 - minval

    where = fill(0, rangelen)
    @inbounds for i = eachindex(x)
        where[x[i] + offs] += 1
    end

    idx = firstindex(x)
    @inbounds for i = maybereverse(1:rangelen)
        lastidx = idx + where[i] - 1
        val = i-offs
        for j = idx:lastidx
            x[j] = val
        end
        idx = lastidx + 1
    end

    return x
end

"""
    sort(v; alg::Algorithm=defalg(v), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)

Variant of [`sort!`](@ref) that returns a sorted copy of `v` leaving `v` itself unmodified.

# Examples
```jldoctest
julia> v = [3, 1, 2];

julia> sort(v)
3-element Array{Int64,1}:
 1
 2
 3

julia> v
3-element Array{Int64,1}:
 3
 1
 2
```
"""
sort(v::AbstractVector; kws...) = sort!(copymutable(v); kws...)

## partialsortperm: the permutation to sort the first k elements of an array ##

"""
    partialsortperm(v, k; by=<transform>, lt=<comparison>, rev=false)

Return a partial permutation `I` of the vector `v`, so that `v[I]` returns values of a fully
sorted version of `v` at index `k`. If `k` is a range, a vector of indices is returned; if
`k` is an integer, a single index is returned. The order is specified using the same
keywords as `sort!`. The permutation is stable, meaning that indices of equal elements
appear in ascending order.

Note that this function is equivalent to, but more efficient than, calling `sortperm(...)[k]`.

# Examples
```jldoctest
julia> v = [3, 1, 2, 1];

julia> v[partialsortperm(v, 1)]
1

julia> p = partialsortperm(v, 1:3)
3-element view(::Array{Int64,1}, 1:3) with eltype Int64:
 2
 4
 3

julia> v[p]
3-element Array{Int64,1}:
 1
 1
 2
```
"""
partialsortperm(v::AbstractVector, k::Union{Integer,OrdinalRange}; kwargs...) =
    partialsortperm!(similar(Vector{eltype(k)}, axes(v,1)), v, k; kwargs..., initialized=false)

"""
    partialsortperm!(ix, v, k; by=<transform>, lt=<comparison>, rev=false, initialized=false)

Like [`partialsortperm`](@ref), but accepts a preallocated index vector `ix` the same size as
`v`, which is used to store (a permutation of) the indices of `v`.

If the index vector `ix` is initialized with the indices of `v` (or a permutation thereof), `initialized` should be set to
`true`.

If `initialized` is `false` (the default), then `ix` is initialized to contain the indices of `v`.

If `initialized` is `true`, but `ix` does not contain (a permutation of) the indices of `v`, the behavior of
`partialsortperm!` is undefined.

(Typically, the indices of `v` will be `1:length(v)`, although if `v` has an alternative array type
with non-one-based indices, such as an `OffsetArray`, `ix` must also be an `OffsetArray` with the same
indices, and must contain as values (a permutation of) these same indices.)

Upon return, `ix` is guaranteed to have the indices `k` in their sorted positions, such that

```julia
partialsortperm!(ix, v, k);
v[ix[k]] == partialsort(v, k)
```

The return value is the `k`th element of `ix` if `k` is an integer, or view into `ix` if `k` is
a range.

# Examples
```jldoctest
julia> v = [3, 1, 2, 1];

julia> ix = Vector{Int}(undef, 4);

julia> partialsortperm!(ix, v, 1)
2

julia> ix = [1:4;];

julia> partialsortperm!(ix, v, 2:3, initialized=true)
2-element view(::Array{Int64,1}, 2:3) with eltype Int64:
 4
 3
```
 """
function partialsortperm!(ix::AbstractVector{<:Integer}, v::AbstractVector,
                          k::Union{Integer, OrdinalRange};
                          lt::Function=isless,
                          by::Function=identity,
                          rev::Union{Bool,Nothing}=nothing,
                          order::Ordering=Forward,
                          initialized::Bool=false)
    if axes(ix,1) != axes(v,1)
        throw(ArgumentError("The index vector is used as a workspace and must have the " *
                            "same length/indices as the source vector, $(axes(ix,1)) != $(axes(v,1))"))
    end
    if !initialized
        @inbounds for i = axes(ix,1)
            ix[i] = i
        end
    end

    # do partial quicksort
    sort!(ix, PartialQuickSort(k), Perm(ord(lt, by, rev, order), v))

    maybeview(ix, k)
end

## sortperm: the permutation to sort an array ##

"""
    sortperm(v; alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)

Return a permutation vector `I` that puts `v[I]` in sorted order. The order is specified
using the same keywords as [`sort!`](@ref). The permutation is guaranteed to be stable even
if the sorting algorithm is unstable, meaning that indices of equal elements appear in
ascending order.

See also [`sortperm!`](@ref).

# Examples
```jldoctest
julia> v = [3, 1, 2];

julia> p = sortperm(v)
3-element Array{Int64,1}:
 2
 3
 1

julia> v[p]
3-element Array{Int64,1}:
 1
 2
 3
```
"""
function sortperm(v::AbstractVector;
                  alg::Algorithm=DEFAULT_UNSTABLE,
                  lt=isless,
                  by=identity,
                  rev::Union{Bool,Nothing}=nothing,
                  order::Ordering=Forward)
    ordr = ord(lt,by,rev,order)
    if ordr === Forward && isa(v,Vector) && eltype(v)<:Integer
        n = length(v)
        if n > 1
            min, max = extrema(v)
            (diff, o1) = sub_with_overflow(max, min)
            (rangelen, o2) = add_with_overflow(diff, oneunit(diff))
            if !o1 && !o2 && rangelen < div(n,2)
                return sortperm_int_range(v, rangelen, min)
            end
        end
    end
    ax = axes(v, 1)
    p = similar(Vector{eltype(ax)}, ax)
    for (i,ind) in zip(eachindex(p), ax)
        p[i] = ind
    end
    sort!(p, alg, Perm(ordr,v))
end

"""
    sortperm!(ix, v; alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward, initialized::Bool=false)

Like [`sortperm`](@ref), but accepts a preallocated index vector `ix`.  If `initialized` is `false`
(the default), `ix` is initialized to contain the values `1:length(v)`.

# Examples
```jldoctest
julia> v = [3, 1, 2]; p = zeros(Int, 3);

julia> sortperm!(p, v); p
3-element Array{Int64,1}:
 2
 3
 1

julia> v[p]
3-element Array{Int64,1}:
 1
 2
 3
```
"""
function sortperm!(x::AbstractVector{<:Integer}, v::AbstractVector;
                   alg::Algorithm=DEFAULT_UNSTABLE,
                   lt=isless,
                   by=identity,
                   rev::Union{Bool,Nothing}=nothing,
                   order::Ordering=Forward,
                   initialized::Bool=false)
    if axes(x,1) != axes(v,1)
        throw(ArgumentError("index vector must have the same length/indices as the source vector, $(axes(x,1)) != $(axes(v,1))"))
    end
    if !initialized
        @inbounds for i = axes(v,1)
            x[i] = i
        end
    end
    sort!(x, alg, Perm(ord(lt,by,rev,order),v))
end

# sortperm for vectors of few unique integers
function sortperm_int_range(x::Vector{<:Integer}, rangelen, minval)
    offs = 1 - minval
    n = length(x)

    where = fill(0, rangelen+1)
    where[1] = 1
    @inbounds for i = 1:n
        where[x[i] + offs + 1] += 1
    end

    #cumsum!(where, where)
    @inbounds for i = 2:length(where)
        where[i] += where[i-1]
    end

    P = Vector{Int}(undef, n)
    @inbounds for i = 1:n
        label = x[i] + offs
        P[where[label]] = i
        where[label] += 1
    end

    return P
end

## sorting multi-dimensional arrays ##

"""
    sort(A; dims::Integer, alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)

Sort a multidimensional array `A` along the given dimension.
See [`sort!`](@ref) for a description of possible
keyword arguments.

To sort slices of an array, refer to [`sortslices`](@ref).

# Examples
```jldoctest
julia> A = [4 3; 1 2]
2×2 Array{Int64,2}:
 4  3
 1  2

julia> sort(A, dims = 1)
2×2 Array{Int64,2}:
 1  2
 4  3

julia> sort(A, dims = 2)
2×2 Array{Int64,2}:
 3  4
 1  2
```
"""
function sort(A::AbstractArray;
              dims::Integer,
              alg::Algorithm=DEFAULT_UNSTABLE,
              lt=isless,
              by=identity,
              rev::Union{Bool,Nothing}=nothing,
              order::Ordering=Forward)
    dim = dims
    order = ord(lt,by,rev,order)
    n = length(axes(A, dim))
    if dim != 1
        pdims = (dim, setdiff(1:ndims(A), dim)...)  # put the selected dimension first
        Ap = permutedims(A, pdims)
        Av = vec(Ap)
        sort_chunks!(Av, n, alg, order)
        permutedims(Ap, invperm(pdims))
    else
        Av = A[:]
        sort_chunks!(Av, n, alg, order)
        reshape(Av, axes(A))
    end
end

@noinline function sort_chunks!(Av, n, alg, order)
    inds = LinearIndices(Av)
    for s = first(inds):n:last(inds)
        sort!(Av, s, s+n-1, alg, order)
    end
    Av
end

"""
    sort!(A; dims::Integer, alg::Algorithm=defalg(A), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)

Sort the multidimensional array `A` along dimension `dims`.
See [`sort!`](@ref) for a description of possible keyword arguments.

To sort slices of an array, refer to [`sortslices`](@ref).

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

# Examples
```jldoctest
julia> A = [4 3; 1 2]
2×2 Array{Int64,2}:
 4  3
 1  2

julia> sort!(A, dims = 1); A
2×2 Array{Int64,2}:
 1  2
 4  3

julia> sort!(A, dims = 2); A
2×2 Array{Int64,2}:
 1  2
 3  4
```
"""
function sort!(A::AbstractArray;
               dims::Integer,
               alg::Algorithm=defalg(A),
               lt=isless,
               by=identity,
               rev::Union{Bool,Nothing}=nothing,
               order::Ordering=Forward)
    ordr = ord(lt, by, rev, order)
    nd = ndims(A)
    k = dims

    1 <= k <= nd || throw(ArgumentError("dimension out of range"))

    remdims = ntuple(i -> i == k ? 1 : size(A, i), nd)
    for idx in CartesianIndices(remdims)
        Av = view(A, ntuple(i -> i == k ? Colon() : idx[i], nd)...)
        sort!(Av, alg, ordr)
    end
    A
end

## fast clever sorting for floats ##

module Float
using ..Sort
using ...Order
using ..Base: @inbounds, AbstractVector, Vector, last, axes

import Core.Intrinsics: slt_int
import ..Sort: sort!
import ...Order: lt, DirectOrdering

const Floats = Union{Float32,Float64}

struct Left <: Ordering end
struct Right <: Ordering end

left(::DirectOrdering) = Left()
right(::DirectOrdering) = Right()

left(o::Perm) = Perm(left(o.order), o.data)
right(o::Perm) = Perm(right(o.order), o.data)

lt(::Left, x::T, y::T) where {T<:Floats} = slt_int(y, x)
lt(::Right, x::T, y::T) where {T<:Floats} = slt_int(x, y)

isnan(o::DirectOrdering, x::Floats) = (x!=x)
isnan(o::Perm, i::Integer) = isnan(o.order,o.data[i])

function nans2left!(v::AbstractVector, o::Ordering, lo::Integer=first(axes(v,1)), hi::Integer=last(axes(v,1)))
    i = lo
    @inbounds while i <= hi && isnan(o,v[i])
        i += 1
    end
    j = i + 1
    @inbounds while j <= hi
        if isnan(o,v[j])
            v[i], v[j] = v[j], v[i]
            i += 1
        end
        j += 1
    end
    return i, hi
end
function nans2right!(v::AbstractVector, o::Ordering, lo::Integer=first(axes(v,1)), hi::Integer=last(axes(v,1)))
    i = hi
    @inbounds while lo <= i && isnan(o,v[i])
        i -= 1
    end
    j = i - 1
    @inbounds while lo <= j
        if isnan(o,v[j])
            v[i], v[j] = v[j], v[i]
            i -= 1
        end
        j -= 1
    end
    return lo, i
end

nans2end!(v::AbstractVector, o::ForwardOrdering) = nans2right!(v,o)
nans2end!(v::AbstractVector, o::ReverseOrdering) = nans2left!(v,o)
nans2end!(v::AbstractVector{<:Integer}, o::Perm{<:ForwardOrdering}) = nans2right!(v,o)
nans2end!(v::AbstractVector{<:Integer}, o::Perm{<:ReverseOrdering}) = nans2left!(v,o)

issignleft(o::ForwardOrdering, x::Floats) = lt(o, x, zero(x))
issignleft(o::ReverseOrdering, x::Floats) = lt(o, x, -zero(x))
issignleft(o::Perm, i::Integer) = issignleft(o.order, o.data[i])

function fpsort!(v::AbstractVector, a::Algorithm, o::Ordering)
    i, j = lo, hi = nans2end!(v,o)
    @inbounds while true
        while i <= j &&  issignleft(o,v[i]); i += 1; end
        while i <= j && !issignleft(o,v[j]); j -= 1; end
        i <= j || break
        v[i], v[j] = v[j], v[i]
        i += 1; j -= 1
    end
    sort!(v, lo, j,  a, left(o))
    sort!(v, i,  hi, a, right(o))
    return v
end

fpsort!(v::AbstractVector, a::Sort.PartialQuickSort, o::Ordering) =
    sort!(v, first(axes(v,1)), last(axes(v,1)), a, o)

sort!(v::AbstractVector{<:Floats}, a::Algorithm, o::DirectOrdering) = fpsort!(v,a,o)
sort!(v::Vector{Int}, a::Algorithm, o::Perm{<:DirectOrdering,<:Vector{<:Floats}}) = fpsort!(v,a,o)

end # module Sort.Float

end # module Sort

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using Test
using StaticArrays, IntervalSets, LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators

@testset "sphere divergence" begin
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(0.0),
        Geometry.ZPoint{FT}(1.0);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = 10)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.SphereDomain(30.0)
    horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 4)
    horztopology = Topologies.Topology2D(horzmesh)
    quad = Spaces.Quadratures.GLL{3 + 1}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)

    coords = Fields.coordinate_field(hv_center_space)
    x = Geometry.UVWVector.(cosd.(coords.lat), 0.0, 0.0)
    div = Operators.Divergence()
    @test norm(div.(x)) < 2e-2
end

function hvspace_3D(
    xlim = (-π, π),
    ylim = (-π, π),
    zlim = (0, 4π),
    xelem = 4,
    yelem = 4,
    zelem = 16,
    npoly = 7,
)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_tags = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = zelem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(xlim[1]) .. Geometry.XPoint{FT}(xlim[2]),
        Geometry.YPoint{FT}(ylim[1]) .. Geometry.YPoint{FT}(ylim[2]),
        x1periodic = true,
        x2periodic = true,
    )
    horzmesh = Meshes.RectilinearMesh(horzdomain, xelem, yelem)
    horztopology = Topologies.Topology2D(horzmesh)

    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

@testset "2D SE, 1D FD Extruded Domain level extraction" begin
    hv_center_space, hv_face_space = hvspace_3D()

    coord = Fields.coordinate_field(hv_face_space)
    @test parent(Fields.field_values(level(coord.x, half))) == parent(
        Fields.field_values(
            Fields.coordinate_field(hv_face_space.horizontal_space).x,
        ),
    )
    @test parent(Fields.field_values(level(coord.z, half))) ==
          parent(
        Fields.field_values(
            Fields.coordinate_field(hv_face_space.horizontal_space).x,
        ),
    ) .* 0
end

@testset "2D SE, 1D FV Extruded Domain ∇ ODE Solve vertical" begin

    hv_center_space, hv_face_space = hvspace_3D()
    V =
        Geometry.UVWVector.(
            zeros(Float64, hv_face_space),
            zeros(Float64, hv_face_space),
            ones(Float64, hv_face_space),
        )

    function rhs!(dudt, u, _, t)
        A = Operators.AdvectionC2C(
            bottom = Operators.SetValue(sin(-t)),
            top = Operators.Extrapolate(),
        )
        return @. dudt = -A(V, u)
    end

    U = sin.(Fields.coordinate_field(hv_center_space).z)
    dudt = zeros(eltype(U), hv_center_space)
    rhs!(dudt, U, nothing, 0.0)

    using OrdinaryDiffEq
    Δt = 0.01
    prob = ODEProblem(rhs!, U, (0.0, 2π))
    sol = solve(prob, SSPRK33(), dt = Δt)

    htopo = Spaces.topology(hv_center_space)
    for h in 1:Topologies.nlocalelems(htopo)
        sol_column_field = ClimaCore.column(sol.u[end], 1, 1, h)
        ref_column_field = ClimaCore.column(U, 1, 1, h)
        @test norm(sol_column_field .- ref_column_field) ≤ 5e-2
    end
end

@testset "2D SE, 1D FV Extruded Domain ∇ ODE Solve horzizontal" begin

    # Advection Equation
    # ∂_t f + c ∂_x f  = 0
    # the solution translates to the right at speed c,
    # so if you you have a periodic domain of size [-π, π]
    # at time t, the solution is f(x - c * t, y)
    # here c == 1, integrate t == 2π or one full period

    function rhs!(dudt, u, _, t)
        # horizontal divergence operator applied to all levels
        hdiv = Operators.Divergence()
        @. dudt = -hdiv(u * Geometry.UVVector(1.0, 1.0))
        Spaces.weighted_dss!(dudt)
        return dudt
    end

    hv_center_space, _ = hvspace_3D()
    U = sin.(Fields.coordinate_field(hv_center_space).x)
    dudt = zeros(eltype(U), hv_center_space)
    rhs!(dudt, U, nothing, 0.0)

    using OrdinaryDiffEq
    Δt = 0.01
    prob = ODEProblem(rhs!, U, (0.0, 2π))
    sol = solve(prob, SSPRK33(), dt = Δt)

    @test norm(U .- sol.u[end]) ≤ 5e-5
end

@testset "2D SE, 1D FV Extruded Domain ∇ ODE Solve diagonally" begin

    # Advection equation in Cartesian domain with
    # uₕ = (cₕ, 0), uᵥ = (0, cᵥ)
    # ∂ₜf + ∇ₕ⋅(uₕ * f) + ∇ᵥ⋅(uᵥ * f)  = 0
    # the solution translates diagonally to the top right and
    # at time t, the solution is f(x - cₕ * t, y - cᵥ * t)
    # here cₕ == cᵥ == 1, integrate t == 2π or one full period.
    # This is only correct if the solution is localized in the vertical
    # as we don't use periodic boundary conditions in the vertical.
    #
    # NOTE: the equation setup is only correct for Cartesian domains!

    hv_center_space, hv_face_space =
        hvspace_3D((-1.0, 1.0), (-1.0, 1.0), (-1.0, 1.0))

    function rhs!(dudt, u, _, t)
        h = u.h
        dh = dudt.h

        # vertical advection no inflow at bottom
        # and outflow at top
        Ic2f = Operators.InterpolateC2F(top = Operators.Extrapolate())
        divf2c = Operators.DivergenceF2C(
            bottom = Operators.SetValue(Geometry.UVWVector(0.0, 0.0, 0.0)),
        )
        # only upward advection
        @. dh = -divf2c(Ic2f(h) * Geometry.UVWVector(0.0, 0.0, 1.0))

        # only horizontal advection
        hdiv = Operators.Divergence()
        @. dh += -hdiv(h * Geometry.UVVector(1.0, 1.0))
        Spaces.weighted_dss!(dh)

        return dudt
    end

    # initial conditions
    function h_init(x_init, y_init, z_init, A = 1.0, σ = 0.2)
        coords = Fields.coordinate_field(hv_center_space)
        h = map(coords) do coord
            A * exp(
                -(
                    (coord.x + x_init)^2 +
                    (coord.y + y_init)^2 +
                    (coord.z + z_init)^2
                ) / (2 * σ^2),
            )
        end
        return h
    end

    using OrdinaryDiffEq
    U = Fields.FieldVector(h = h_init(0.5, 0.5, 0.5))
    Δt = 0.01
    t_end = 1.0
    prob = ODEProblem(rhs!, U, (0.0, t_end))
    sol = solve(prob, SSPRK33(), dt = Δt)

    h_end = h_init(0.5 - t_end, 0.5 - t_end, 0.5 - t_end)
    @test norm(h_end .- sol.u[end].h) ≤ 8e2
end

ext: jl


text:
using Surrogates
using LinearAlgebra
using Zygote
using Test
#using Zygote: @nograd
#=
#FORWARD
###### 1D ######
lb = 0.0
ub = 10.0
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x^2
y = f.(x)

#Radials
my_rad = RadialBasis(x,y,lb,ub,x->norm(x),2)
g = x -> ForwardDiff.derivative(my_rad,x)
g(5.0)

#Kriging
p = 1.5
my_krig = Kriging(x,y,p)
g = x -> ForwardDiff.derivative(my_krig,x)
g(5.0)

#Linear Surrogate
my_linear = LinearSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.derivative(my_linear,x)
g(5.0)

#Inverse distance
p = 1.4
my_inverse = InverseDistanceSurrogate(x,y,p,lb,ub)
g = x -> ForwardDiff.derivative(my_inverse,x)
g(5.0)

#Lobachevsky
n = 4
α = 2.4
my_loba = LobachevskySurrogate(x,y,α,n,lb,ub)
g = x -> ForwardDiff.derivative(my_loba,x)
g(5.0)

#Second order polynomial
my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.derivative(my_second,x)
g(5.0)

###### ND ######
lb = [0.0,0.0]
ub = [10.0,10.0]
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]
y = f.(x)

#Radials
my_rad = RadialBasis(x,y,[lb,ub],z->norm(z),2)
g = x -> ForwardDiff.gradient(my_rad,x)
g([2.0,5.0])

#Kriging
theta = [2.0,2.0]
p = [1.9,1.9]
my_krig = Kriging(x,y,p,theta)
g = x -> ForwardDiff.gradient(my_krig,x)
g([2.0,5.0])

#Linear Surrogate
my_linear = LinearSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.gradient(my_linear,x)
g([2.0,5.0])

#Inverse Distance
p = 1.4
my_inverse = InverseDistanceSurrogate(x,y,p,lb,ub)
g = x -> ForwardDiff.gradient(my_inverse,x)
g([2.0,5.0])

#Lobachevsky
alpha = [1.4,1.4]
n = 4
my_loba_ND = LobachevskySurrogate(x,y,alpha,n,lb,ub)
g = x -> ForwardDiff.gradient(my_loba_ND,x)
g([2.0,5.0])

#Second order polynomial
my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.gradient(my_second,x)
g([2.0,5.0])
=#

##############
### ZYGOTE ###
##############

############
#### 1D ####
############
lb = 0.0
ub = 3.0
n = 10
x = sample(n,lb,ub,SobolSample())
f = x -> x^2
y = f.(x)

#Radials
@testset "Radials 1D" begin
    my_rad = RadialBasis(x,y,lb,ub,rad = linearRadial())
    g = x -> my_rad'(x)
    g(5.0)
end

# #Kriging
@testset "Kriging 1D" begin
    my_p = 1.5
    my_krig = Kriging(x,y,lb,ub,p=my_p)
    g = x -> my_krig'(x)
    g(5.0)
end

# #Linear Surrogate
@testset "Linear Surrogate" begin
    my_linear = LinearSurrogate(x,y,lb,ub)
    g = x -> my_linear'(x)
    g(5.0)
end

#Inverse distance
@testset "Inverse Distance" begin
    my_p = 1.4
    my_inverse = InverseDistanceSurrogate(x,y,lb,ub,p=my_p)
    g = x -> my_inverse'(x)
    g(5.0)
end

#Second order polynomial
@testset "Second Order Polynomial" begin
    my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
    g = x -> my_second'(x)
    g(5.0)
end

#Lobachevsky
@testset "Lobachevsky" begin
    n = 4
    α = 2.4
    my_loba = LobachevskySurrogate(x,y,lb,ub, alpha = α, n = 4)
    g = x -> my_loba'(x)
    g(0.0)
end

#Wendland
@testset "Wendland" begin
    my_wend = Wendland(x,y,lb,ub)
    g = x -> my_wend'(x)
    g(3.0)
end

# MOE and VariableFidelity for free because they are Linear combinations
# of differentiable surrogates

# #Gek
@testset "GEK" begin
    n = 10
    lb = 0.0
    ub = 5.0
    x = sample(n,lb,ub,SobolSample())
    f = x-> x^2
    y1 = f.(x)
    der = x->2*x
    y2 = der.(x)
    y = vcat(y1,y2)
    my_gek = GEK(x,y,lb,ub)
    g = x-> my_gek'(x)
    g(3.0)
end

################
###### ND ######
################

lb = [0.0,0.0]
ub = [10.0,10.0]
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]
y = f.(x)

#Radials
@testset "Radials ND" begin
    my_rad = RadialBasis(x,y,lb,ub,rad = linearRadial(), scale_factor = 2.1)
    g = x -> Zygote.gradient(my_rad,x)
    g((2.0,5.0))
end

# Kriging
@testset "Kriging ND" begin
    my_theta = [2.0,2.0]
    my_p = [1.9,1.9]
    my_krig = Kriging(x,y,lb,ub,p=my_p,theta=my_theta)
    g = x -> Zygote.gradient(my_krig,x)
    g((2.0,5.0))
end

# Linear Surrogate
@testset "Linear Surrogate ND" begin
    my_linear = LinearSurrogate(x,y,lb,ub)
    g = x -> Zygote.gradient(my_linear,x)
    g((2.0,5.0))
end

#Inverse Distance
@testset "Inverse Distance ND" begin
    my_p = 1.4
    my_inverse = InverseDistanceSurrogate(x,y,lb,ub,p=my_p)
    g = x -> Zygote.gradient(my_inverse,x)
    g((2.0,5.0))
end

#Lobachevsky not working yet weird issue with Zygote @nograd
#=
Zygote.refresh()
alpha = [1.4,1.4]
n = 4
my_loba_ND = LobachevskySurrogate(x,y,alpha,n,lb,ub)
g = x -> Zygote.gradient(my_loba_ND,x)
g((2.0,5.0))
=#

# Second order polynomial mutating arrays
@testset "SecondOrderPolynomialSurrogate ND" begin
    my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
    g = x -> Zygote.gradient(my_second,x)
    g((2.0,5.0))
end

#wendland
@testset "Wendland ND" begin
    my_wend_ND = Wendland(x,y,lb,ub)
    g = x -> Zygote.gradient(my_wend_ND,x)
    g((2.0,5.0))
end
# #MOE and VariableFidelity for free because they are Linear combinations
# #of differentiable surrogates


text:
module LIKWID
import Base.Threads
using Libdl
using Unitful

# liblikwid
const liblikwid = "liblikwid"
include("LibLikwid.jl")

# Julia types
include("types.jl")

# Note: underscore prefix -> liblikwid API
const topo_initialized = Ref{Bool}(false)
const numa_initialized = Ref{Bool}(false)
const affinity_initialized = Ref{Bool}(false)
const timer_initialized = Ref{Bool}(false)
const power_initialized = Ref{Bool}(false)
const config_initialized = Ref{Bool}(false)
const access_initialized = Ref{Bool}(false)
const perfmon_initialized = Ref{Bool}(false)
const gputopo_initialized = Ref{Bool}(false)
const nvmon_initialized = Ref{Bool}(false)
const _cputopo = Ref{Union{LibLikwid.CpuTopology,Nothing}}(nothing) # (Julia) API struct
const cputopo = Ref{Union{CpuTopology,Nothing}}(nothing) # Julia struct
const _cpuinfo = Ref{Union{LibLikwid.CpuInfo,Nothing}}(nothing) # (Julia) API struct
const cpuinfo = Ref{Union{CpuInfo,Nothing}}(nothing) # Julia struct
const _numainfo = Ref{Union{LibLikwid.NumaTopology,Nothing}}(nothing) # (Julia) API struct
const numainfo = Ref{Union{NumaTopology,Nothing}}(nothing) # Julia struct
const _affinity = Ref{Union{LibLikwid.AffinityDomains,Nothing}}(nothing) # (Julia) API struct
const affinity = Ref{Union{AffinityDomains,Nothing}}(nothing) # Julia struct
const _powerinfo = Ref{Union{LibLikwid.PowerInfo,Nothing}}(nothing) # (Julia) API struct
const powerinfo = Ref{Union{PowerInfo,Nothing}}(nothing) # Julia struct
const _config = Ref{Union{LibLikwid.Likwid_Configuration,Nothing}}(nothing) # (Julia) API struct
const config = Ref{Union{Likwid_Configuration,Nothing}}(nothing) # Julia struct
const _gputopo = Ref{Union{LibLikwid.GpuTopology,Nothing}}(nothing) # (Julia) API struct
const gputopo = Ref{Union{GpuTopology,Nothing}}(nothing) # Julia struct

const likwid_gpusupport = Ref{Union{Nothing,Bool}}(nothing)

# functions
include("topology.jl")
include("numa.jl")
include("configuration.jl")
include("affinity.jl")
include("timer.jl")
include("thermal.jl")
include("power.jl")
include("access.jl")
include("perfmon.jl")
include("misc.jl")
include("marker.jl")
import .Marker
import .Marker: region, @region
export Marker, region, @region
include("markerfile.jl")
include("topology_gpu.jl")
include("nvmon.jl")
include("marker_gpu.jl")
import .GPUMarker
import .GPUMarker: gpuregion, @gpuregion
export GPUMarker, gpuregion, @gpuregion
include("frequency.jl")

function init(; gpu=false)
    Marker.init()
    init_topology()
    init_numa()
    init_affinity()
    PerfMon.init()
    HPM.init()
    Timer.init()
    Power.init()
    Freq.init()
    if gpu && gpusupport()
        GPUMarker.init()
        init_topology_gpu()
        NvMon.init()
    end
    return nothing
end

function finalize(; gpu=true)
    Marker.close()
    finalize_topology()
    finalize_numa()
    finalize_affinity()
    PerfMon.finalize()
    HPM.finalize()
    Timer.finalize()
    Power.finalize()
    Freq.finalize()
    if gpu && gpusupport()
        GPUMarker.close()
        finalize_topology_gpu()
        NvMon.finalize()
    end
    return nothing
end

end


text:
"""
    ArrayIndex{N}

Subtypes of `ArrayIndex` represent series of transformations for a provided index to some
buffer which is typically accomplished with square brackets (e.g., `buffer[index[inds...]]`).
The only behavior that is required of a subtype of `ArrayIndex` is the ability to transform
individual index elements (i.e. not collections). This does not guarantee bounds checking or
the ability to iterate (although additional functionallity may be provided for specific
types).
"""
abstract type ArrayIndex{N} end

const MatrixIndex = ArrayIndex{2}

const VectorIndex = ArrayIndex{1}

Base.ndims(::Type{<:ArrayIndex{N}}) where {N} = N

struct BidiagonalIndex <: MatrixIndex
    count::Int
    isup::Bool
end

struct TridiagonalIndex <: MatrixIndex
    count::Int# count==nsize+nsize-1+nsize-1
    nsize::Int
    isrow::Bool
end

struct BandedMatrixIndex <: MatrixIndex
    count::Int
    rowsize::Int
    colsize::Int
    bandinds::Array{Int,1}
    bandsizes::Array{Int,1}
    isrow::Bool
end

function _bandsize(bandind, rowsize, colsize)
    -(rowsize - 1) <= bandind <= colsize - 1 || throw(ErrorException("Invalid Bandind"))
    if (bandind * (colsize - rowsize) > 0) & (abs(bandind) <= abs(colsize - rowsize))
        return min(rowsize, colsize)
    elseif bandind * (colsize - rowsize) <= 0
        return min(rowsize, colsize) - abs(bandind)
    else
        return min(rowsize, colsize) - abs(bandind) + abs(colsize - rowsize)
    end
end

function BandedMatrixIndex(rowsize, colsize, lowerbandwidth, upperbandwidth, isrow)
    upperbandwidth > -lowerbandwidth || throw(ErrorException("Invalid Bandwidths"))
    bandinds = upperbandwidth:-1:-lowerbandwidth
    bandsizes = [_bandsize(band, rowsize, colsize) for band in bandinds]
    BandedMatrixIndex(sum(bandsizes), rowsize, colsize, bandinds, bandsizes, isrow)
end

struct BlockBandedMatrixIndex <: MatrixIndex
    count::Int
    refinds::Array{Int,1}
    refcoords::Array{Int,1}# storing col or row inds at ref points
    isrow::Bool
end

function BlockBandedMatrixIndex(nrowblock, ncolblock, rowsizes, colsizes, l, u)
    blockrowind = BandedMatrixIndex(nrowblock, ncolblock, l, u, true)
    blockcolind = BandedMatrixIndex(nrowblock, ncolblock, l, u, false)
    sortedinds = sort(
        [(blockrowind[i], blockcolind[i]) for i = 1:length(blockrowind)],
        by = x -> x[1],
    )
    sort!(sortedinds, by = x -> x[2], alg = InsertionSort)# stable sort keeps the second index in order
    refinds = Array{Int,1}()
    refrowcoords = Array{Int,1}()
    refcolcoords = Array{Int,1}()
    rowheights = pushfirst!(copy(rowsizes), 1)
    cumsum!(rowheights, rowheights)
    blockheight = 0
    blockrow = 1
    blockcol = 1
    currenti = 1
    lastrowind = sortedinds[1][1] - 1
    lastcolind = sortedinds[1][2]
    for ind in sortedinds
        rowind, colind = ind
        if colind == lastcolind
            if rowind > lastrowind
                blockheight += rowsizes[rowind]
            end
        else
            for j = blockcol:blockcol+colsizes[lastcolind]-1
                push!(refinds, currenti)
                push!(refrowcoords, blockrow)
                push!(refcolcoords, j)
                currenti += blockheight
            end
            blockcol += colsizes[lastcolind]
            blockrow = rowheights[rowind]
            blockheight = rowsizes[rowind]
        end
        lastcolind = colind
        lastrowind = rowind
    end
    for j = blockcol:blockcol+colsizes[lastcolind]-1
        push!(refinds, currenti)
        push!(refrowcoords, blockrow)
        push!(refcolcoords, j)
        currenti += blockheight
    end
    push!(refinds, currenti)# guard
    push!(refrowcoords, -1)
    push!(refcolcoords, -1)
    rowindobj = BlockBandedMatrixIndex(currenti - 1, refinds, refrowcoords, true)
    colindobj = BlockBandedMatrixIndex(currenti - 1, refinds, refcolcoords, false)
    rowindobj, colindobj
end

struct BandedBlockBandedMatrixIndex <: MatrixIndex
    count::Int
    refinds::Array{Int,1}
    refcoords::Array{Int,1}# storing col or row inds at ref points
    reflocalinds::Array{BandedMatrixIndex,1}
    isrow::Bool
end

function BandedBlockBandedMatrixIndex(
    nrowblock,
    ncolblock,
    rowsizes,
    colsizes,
    l,
    u,
    lambda,
    mu,
)
    blockrowind = BandedMatrixIndex(nrowblock, ncolblock, l, u, true)
    blockcolind = BandedMatrixIndex(nrowblock, ncolblock, l, u, false)
    sortedinds = sort(
        [(blockrowind[i], blockcolind[i]) for i = 1:length(blockrowind)],
        by = x -> x[1],
    )
    sort!(sortedinds, by = x -> x[2], alg = InsertionSort)# stable sort keeps the second index in order
    rowheights = pushfirst!(copy(rowsizes), 1)
    cumsum!(rowheights, rowheights)
    colwidths = pushfirst!(copy(colsizes), 1)
    cumsum!(colwidths, colwidths)
    currenti = 1
    refinds = Array{Int,1}()
    refrowcoords = Array{Int,1}()
    refcolcoords = Array{Int,1}()
    reflocalrowinds = Array{BandedMatrixIndex,1}()
    reflocalcolinds = Array{BandedMatrixIndex,1}()
    for ind in sortedinds
        rowind, colind = ind
        localrowind =
            BandedMatrixIndex(rowsizes[rowind], colsizes[colind], lambda, mu, true)
        localcolind =
            BandedMatrixIndex(rowsizes[rowind], colsizes[colind], lambda, mu, false)
        push!(refinds, currenti)
        push!(refrowcoords, rowheights[rowind])
        push!(refcolcoords, colwidths[colind])
        push!(reflocalrowinds, localrowind)
        push!(reflocalcolinds, localcolind)
        currenti += localrowind.count
    end
    push!(refinds, currenti)
    push!(refrowcoords, -1)
    push!(refcolcoords, -1)
    rowindobj = BandedBlockBandedMatrixIndex(
        currenti - 1,
        refinds,
        refrowcoords,
        reflocalrowinds,
        true,
    )
    colindobj = BandedBlockBandedMatrixIndex(
        currenti - 1,
        refinds,
        refcolcoords,
        reflocalcolinds,
        false,
    )
    rowindobj, colindobj
end

struct StrideIndex{N,R,C,S,O,O1} <: ArrayIndex{N}
    strides::S
    offsets::O
    offset1::O1

    function StrideIndex{N,R,C}(s::S, o::O, o1::O1) where {N,R,C,S<:Tuple{Vararg{<:CanonicalInt,N}},O<:Tuple{Vararg{<:CanonicalInt,N}},O1}
        return new{N,R::NTuple{N,Int},C::Int,S,O,O1}(s, o, o1)
    end
    function StrideIndex{N,R,C}(a::A) where {N,R,C,A}
        return StrideIndex{N,R,C}(strides(a),offsets(a),offset1(a))
    end
    function StrideIndex(a::A) where {A}
        return StrideIndex{ndims(A),known(stride_rank(A)), known(contiguous_axis(A))}(a)
    end
end

Base.firstindex(i::Union{TridiagonalIndex,BandedBlockBandedMatrixIndex,BandedMatrixIndex,BidiagonalIndex,BlockBandedMatrixIndex}) = 1
Base.lastindex(i::Union{TridiagonalIndex,BandedBlockBandedMatrixIndex,BandedMatrixIndex,BidiagonalIndex,BlockBandedMatrixIndex}) = i.count
Base.length(i::Union{TridiagonalIndex,BandedBlockBandedMatrixIndex,BandedMatrixIndex,BidiagonalIndex,BlockBandedMatrixIndex}) = i.count
function Base.getindex(ind::BidiagonalIndex, i::Int)
    1 <= i <= ind.count || throw(BoundsError(ind, i))
    if ind.isup
        ii = i + 1
    else
        ii = i + 1 + 1
    end
    convert(Int, floor(ii / 2))
end

function Base.getindex(ind::TridiagonalIndex, i::Int)
    1 <= i <= ind.count || throw(BoundsError(ind, i))
    offsetu = ind.isrow ? 0 : 1
    offsetl = ind.isrow ? 1 : 0
    if 1 <= i <= ind.nsize
        return i
    elseif ind.nsize < i <= ind.nsize + ind.nsize - 1
        return i - ind.nsize + offsetu
    else
        return i - (ind.nsize + ind.nsize - 1) + offsetl
    end
end

function Base.getindex(ind::BandedMatrixIndex, i::Int)
    1 <= i <= ind.count || throw(BoundsError(ind, i))
    _i = i
    p = 1
    while _i - ind.bandsizes[p] > 0
        _i -= ind.bandsizes[p]
        p += 1
    end
    bandind = ind.bandinds[p]
    startfromone = !xor(ind.isrow, (bandind > 0))
    if startfromone
        return _i
    else
        return _i + abs(bandind)
    end
end

function Base.getindex(ind::BlockBandedMatrixIndex, i::Int)
    1 <= i <= ind.count || throw(BoundsError(ind, i))
    p = 1
    while i - ind.refinds[p] >= 0
        p += 1
    end
    p -= 1
    _i = i - ind.refinds[p]
    if ind.isrow
        return ind.refcoords[p] + _i
    else
        return ind.refcoords[p]
    end
end

function Base.getindex(ind::BandedBlockBandedMatrixIndex, i::Int)
    1 <= i <= ind.count || throw(BoundsError(ind, i))
    p = 1
    while i - ind.refinds[p] >= 0
        p += 1
    end
    p -= 1
    _i = i - ind.refinds[p] + 1
    ind.reflocalinds[p][_i] + ind.refcoords[p] - 1
end

@inline function Base.getindex(x::StrideIndex{N}, i::CanonicalInt...) where {N}
    return _strides2int(x.offsets, x.strides, Tuple(i)) + x.offset1
end
@inline function _strides2int(o::Tuple, s::Tuple, i::Tuple)
    return ((first(i) - first(o)) * first(s)) + _strides2int(tail(o), tail(s), tail(i))
end
_strides2int(::Tuple{}, ::Tuple{}, ::Tuple{}) = static(0)


text:
####
#### OpTerm and OpSum backed by SparseVector
####

# using ..AbstractPaulis: PauliPhaseCounts

SparseArraysN.neutral(T::Type{<:AbstractOp}) = one(T)
SparseArraysN.neutral(x::AbstractOp) = one(x)
SparseArraysN.isneutral(x::AbstractOp) = isone(x)
SparseArraysN.neutrals(T::Type{<:AbstractOp}, args...) = ones(T, args...)
SparseArraysN.neutrals(T::Type{<:AbstractOp}, args::Integer...) = ones(T, args...)

const SparseOpTerm = OpTerm{T, V} where {T<:AbstractOp, V<:SparseArraysN.SparseVector{T}}

#function Base.:*(t1::OpTerm{T, V}, t2::OpTerm{T, V}) where {T<:AbstractOp, V<:SparseArraysN.SparseVector{T}}
function Base.:*(t1::SparseOpTerm, t2::SparseOpTerm)
    (t_out, phase) = mul(t1.ops, t2.ops)
    return OpTerm(t_out, phase * t1.coeff * t2.coeff)
end

"""
    mul(v1::SparseArraysN.SparseVector{T}, v2::SparseArraysN.SparseVector{T}) where {T <: AbstractOp}

Multiply `v1` and `v2`, factorwise, returning a tuple of a `SparseVector` and the phase accumulated
from each multiplication.
"""
function mul(v1::SparseArraysN.SparseVector{T}, v2::SparseArraysN.SparseVector{T}) where {T <: AbstractOp}
    vals = T[]
    inds = empty(v1.nzind)
    phase_data = AbstractOps.phase_data(T)
    i1 = 1  # index into terms in t1
    i2 = 1
    while i1 <= weight(v1) && i2 <= weight(v2)
        ind1 = v1.nzind[i1]
        ind2 = v2.nzind[i2]
        val1 = v1.nzval[i1]
        val2 = v2.nzval[i2]
        ## ind1 corresponds to non-structural identity in v2
        ## Multiplying v1[ind1] by identity gives v1[ind1], so we just push
        ## this value to the result.
        if ind1 < ind2
            push!(inds, ind1)
            push!(vals, val1)
            i1 += 1
        ## Symmetric to the case above.
        elseif ind1 > ind2
            push!(inds, ind2)
            push!(vals, val2)
            i2 += 1
        else  # tt1.op and tt2.op operate on the same index (degree of freedom)
            ## Only here do we multiply non-trivial factors.
            new_op = val1 * val2  # product up to a phase
            if iszero(new_op)  # if any factor vanishes, the term vanishes.
                _phase = 0
                return (SparseArraysN.SparseVector(0, empty!(inds), empty!(vals)), _phase)
            end
            i1 += 1
            i2 += 1
            if isone(new_op)  # Identity is not stored in sparse representation
                continue
            end
#            phase_data = AbstractOps.accumulate_phase(phase_data, val1, val2)
            AbstractOps.accumulate_phase!(phase_data, val1, val2)
            push!(inds, ind1)
            push!(vals, new_op)
        end
    end
    ## Include remaining factors from the longer string.
    if i1 <= weight(v1)
        for i in i1:weight(v1)
            push!(inds, v1.nzind[i])
            push!(vals, v1.nzval[i])
        end
    elseif i2 <= weight(v2)
        for i in i2:weight(v2)
            push!(inds, v2.nzind[i])
            push!(vals, v2.nzval[i])
        end
    end
    n = isempty(inds) ? 0 : inds[end]
    _phase = AbstractOps.compute_phase(phase_data, v1, v2)
    return (SparseArraysN.SparseVector(n, inds, vals), _phase)
end

function _show_sparse_term(io, term)
    ops = term.ops
    xnnz = length(ops.nzind)
    for i in eachindex(ops.nzind)
        _show_op_plain(io, ops.nzval[i])
        print(io, ops.nzind[i])
        if i < length(ops.nzind)
            print(io, " ")
        end
    end
    if xnnz != 0
        print(io, " * ")
    end
    if term.coeff isa Real
        print(io, term.coeff)
    else
        print(io, "(", term.coeff, ")")
    end
end

function Base.show(io::IO, term::OpTerm{T, <:SparseArraysN.SparseVector}) where {T}
    ops = term.ops
    xnnz = length(ops.nzind)
    print(io, length(term), "-element ", typeof(term), " with ", xnnz,
          " stored ", xnnz == 1 ? "entry" : "entries")
    print(io, ":\n")
    _show_sparse_term(io, term)
end

function Base.show(io::IO, opsum::OpSum{T, V}) where {T<:AbstractOp, V <: Vector{<:SparseArraysN.SparseVector}}
    (m, n) = size(opsum)
    print(io, m, "x", n, " ", typeof(opsum), ":\n")
    for i in eachindex(opsum)
        _show_sparse_term(io, opsum[i])
        if i != lastindex(opsum)
            print(io, "\n")
        end
    end
end

## TODO: should we remove this? I think the remaining method would then return a copy
sparse_op(term::OpTerm{T, V}) where {T<:AbstractOp, V<:SparseArraysN.SparseVector{T}} =
    term

"""
    sparse_op(x::OpTerm)
    sparse_op(x::OpSum)

Convert `x` to a sparse representation. Note that the function `sparse` instead may convert `x` to
a sparse matrix.
"""
function sparse_op(term::OpTerm)
    return OpTerm(SparseArraysN.sparse(term.ops), term.coeff)
end

"""
    dense_op(term::OpTerm)
    dense_op(_sum::OpSum)

Convert `term` or `_sum` to a dense representation.
"""
function dense_op(term::OpTerm)
    return OpTerm(Vector(term.ops), term.coeff)
end

function _convert_op(_sum::OpSum, convert_func)
    sum_out = OpSum([convert_func(first(_sum))]; already_sorted=true)
    @inbounds for i in 2:length(_sum)
        push!(sum_out, convert_func(_sum[i]))
    end
    return sum_out
end

dense_op(_sum::OpSum) = _convert_op(_sum, dense_op)
sparse_op(_sum::OpSum) = _convert_op(_sum, sparse_op)


text:
Materials = (
    tissueplastic = Material(
        "A-150 Tissue-Equivalent Plastic",
        0.54903,
        65.1eV,
        1.127g / cm^3,
        Dict(
            1 => 0.10133,
            6 => 0.775498,
            7 => 0.035057,
            8 => 0.052315,
            9 => 0.017423,
            20 => 0.018377,
        ),
    ),
    adipose = Material(
        "Adipose Tissue (ICRU-44)",
        0.55579,
        64.8eV,
        0.95g / cm^3,
        Dict(
            1 => 0.114,
            6 => 0.598,
            7 => 0.007,
            8 => 0.278,
            11 => 0.001,
            16 => 0.001,
            17 => 0.001,
        ),
    ),
    air = Material(
        "Air, Dry (near sea level)",
        0.49919,
        85.7eV,
        0.001205g / cm^3,
        Dict(6 => 0.000124, 7 => 0.755268, 8 => 0.231781, 18 => 0.012827),
    ),
    alanine = Material(
        "Alanine",
        0.53876,
        71.9eV,
        1.424g / cm^3,
        Dict(1 => 0.079192, 6 => 0.404437, 7 => 0.157213, 8 => 0.359157),
    ),
    boneplastic = Material(
        "B-100 Bone-Equivalent Plastic",
        0.5274,
        85.9eV,
        1.45g / cm^3,
        Dict(
            1 => 0.065473,
            6 => 0.536942,
            7 => 0.0215,
            8 => 0.032084,
            9 => 0.167415,
            20 => 0.176585,
        ),
    ),
    bakelite = Material(
        "Bakelite",
        0.52792,
        72.4eV,
        1.25g / cm^3,
        Dict(1 => 0.057444, 6 => 0.774589, 8 => 0.167968),
    ),
    wholeblood = Material(
        "Blood, Whole (ICRU-44)",
        0.54999,
        75.2eV,
        1.06g / cm^3,
        Dict(
            1 => 0.102,
            6 => 0.11,
            7 => 0.033,
            8 => 0.745,
            11 => 0.001,
            15 => 0.001,
            16 => 0.002,
            17 => 0.003,
            19 => 0.002,
            26 => 0.001,
        ),
    ),
    corticalbone = Material(
        "Bone, Cortical (ICRU-44)",
        0.51478,
        112.0eV,
        1.92g / cm^3,
        Dict(
            1 => 0.034,
            6 => 0.155,
            7 => 0.042,
            8 => 0.435,
            11 => 0.001,
            12 => 0.002,
            15 => 0.103,
            16 => 0.003,
            20 => 0.225,
        ),
    ),
    brain = Material(
        "Brain, Grey/White Matter (ICRU-44)",
        0.55239,
        73.9eV,
        1.04g / cm^3,
        Dict(
            1 => 0.107,
            6 => 0.145,
            7 => 0.022,
            8 => 0.712,
            11 => 0.002,
            15 => 0.004,
            16 => 0.002,
            17 => 0.003,
            19 => 0.003,
        ),
    ),
    breasttissue = Material(
        "Breast Tissue (ICRU-44)",
        0.55196,
        70.3eV,
        1.02g / cm^3,
        Dict(
            1 => 0.106,
            6 => 0.332,
            7 => 0.03,
            8 => 0.527,
            11 => 0.001,
            15 => 0.001,
            16 => 0.002,
            17 => 0.001,
        ),
    ),
    airplastic = Material(
        "C-552 Air-equivalent Plastic",
        0.49969,
        86.8eV,
        1.76g / cm^3,
        Dict(
            1 => 0.024681,
            6 => 0.50161,
            8 => 0.004527,
            9 => 0.465209,
            14 => 0.003973,
        ),
    ),
    CdTe = Material(
        "Cadmium Telluride",
        0.41665,
        539.3eV,
        6.2g / cm^3,
        Dict(48 => 0.468358, 52 => 0.531642),
    ),
    CaF = Material(
        "Calcium Fluoride",
        0.48671,
        166.0eV,
        3.18g / cm^3,
        Dict(9 => 0.486672, 20 => 0.513328),
    ),
    CaSO4 = Material(
        "Calcium Sulfate",
        0.49948,
        152.3eV,
        2.96g / cm^3,
        Dict(8 => 0.470081, 16 => 0.235534, 20 => 0.294385),
    ),
    ammoniumsulfate = Material(
        "15 mmol L-1 Ceric Ammonium Sulfate Solution",
        0.55282,
        76.7eV,
        1.030E+00g / cm^3,
        Dict(
            1 => 0.107694,
            7 => 0.000816,
            8 => 0.875172,
            16 => 0.01427958,
            58 => 0.00204,
        ),
    ),
    CsI = Material(
        "Cesium Iodide",
        0.41569,
        553.1eV,
        4.51g / cm^3,
        Dict(53 => 0.488451, 55 => 0.511549),
    ),
    concrete = Material(
        "Concrete, Ordinary",
        0.50932,
        124.5eV,
        2.3g / cm^3,
        Dict(
            1 => 0.0221,
            6 => 0.002484,
            8 => 0.57493,
            11 => 0.015208,
            12 => 0.001266,
            13 => 0.019953,
            14 => 0.304627,
            19 => 0.010045,
            20 => 0.042951,
            26 => 0.006435,
        ),
    ),
    concretebarite = Material(
        "Concrete, Barite (TYPE BA)",
        0.45714,
        248.2eV,
        3.35g / cm^3,
        Dict(
            1 => 0.003585,
            8 => 0.311622,
            12 => 0.001195,
            13 => 0.004183,
            14 => 0.010457,
            16 => 0.107858,
            20 => 0.050194,
            26 => 0.047505,
            56 => 0.4634,
        ),
    ),
    eyelens = Material(
        "Eye Lens (ICRU-44)",
        0.54709,
        74.3eV,
        1.07g / cm^3,
        Dict(
            1 => 0.096,
            6 => 0.195,
            7 => 0.057,
            8 => 0.646,
            11 => 0.001,
            15 => 0.001,
            16 => 0.003,
            17 => 0.001,
        ),
    ),
    ferroussulfate = Material(
        "Ferrous Sulfate Standard Fricke",
        0.55334,
        76.3eV,
        1.024g / cm^3,
        Dict(
            1 => 0.108376,
            8 => 0.878959,
            11 => 2.2e-05,
            16 => 0.012553,
            17 => 3.5e-05,
            26 => 5.5e-05,
        ),
    ),
    GOS = Material(
        "Gadolinium Oxysulfide",
        0.42265,
        493.3eV,
        7.44g / cm^3,
        Dict(8 => 0.084527, 16 => 0.084704, 64 => 0.830769),
    ),
    gafchromicsensor = Material(
        "Gafchromic Sensor",
        0.54384,
        67.2eV,
        1.3g / cm^3,
        Dict(1 => 0.0897, 6 => 0.6058, 7 => 0.1122, 8 => 0.1923),
    ),
    GaAs = Material(
        "Gallium Arsenide",
        0.44246,
        384.9eV,
        5.31g / cm^3,
        Dict(31 => 0.48203, 33 => 0.51797),
    ),
    pyrex = Material(
        "Glass, Borosilicate (Pyrex)",
        0.49707,
        134.0eV,
        2.23g / cm^3,
        Dict(
            5 => 0.040066,
            8 => 0.539559,
            11 => 0.028191,
            13 => 0.011644,
            14 => 0.37722,
            19 => 0.003321,
        ),
    ),
    leadglass = Material(
        "Glass, Lead",
        0.42101,
        526.4eV,
        6.22g / cm^3,
        Dict(
            8 => 0.156453,
            14 => 0.080866,
            22 => 0.008092,
            33 => 0.002651,
            82 => 0.751938,
        ),
    ),
    LiF = Material(
        "Lithium Fluride",
        0.46262,
        94.0eV,
        2.635g / cm^3,
        Dict(3 => 0.267585, 9 => 0.732415),
    ),
    LiB4O7 = Material(
        "Lithium Tetraborate",
        0.48485,
        94.6eV,
        2.44g / cm^3,
        Dict(3 => 0.082081, 5 => 0.255715, 8 => 0.662204),
    ),
    lung = Material(
        "Lung Tissue (ICRU-44)",
        0.55048,
        75.2eV,
        1.05g / cm^3,
        Dict(
            1 => 0.103,
            6 => 0.105,
            7 => 0.031,
            8 => 0.749,
            11 => 0.002,
            15 => 0.002,
            16 => 0.003,
            17 => 0.003,
            19 => 0.002,
        ),
    ),
    MgB4O7 = Material(
        "Magnesium Tetroborate",
        0.49012,
        108.3eV,
        2.53g / cm^3,
        Dict(5 => 0.24087, 8 => 0.623762, 12 => 0.135367),
    ),
    MgI2 = Material(
        "Mercuric Iodide",
        0.40933,
        684.5eV,
        6.36g / cm^3,
        Dict(53 => 0.55856, 80 => 0.44144),
    ),
    muscle = Material(
        "Muscle, Skeletal (ICRU-44)",
        0.55,
        74.6eV,
        1.05g / cm^3,
        Dict(
            1 => 0.102,
            6 => 0.143,
            7 => 0.034,
            8 => 0.71,
            11 => 0.001,
            15 => 0.002,
            16 => 0.003,
            17 => 0.001,
            19 => 0.004,
        ),
    ),
    ovary = Material(
        "Ovary (ICRU-44)",
        0.55149,
        75.0eV,
        1.05g / cm^3,
        Dict(
            1 => 0.105,
            6 => 0.093,
            7 => 0.024,
            8 => 0.768,
            11 => 0.002,
            15 => 0.002,
            16 => 0.002,
            17 => 0.002,
            19 => 0.002,
        ),
    ),
    emulsionkodak = Material(
        "Photographic Emulsion (Kodak Type AA)",
        0.48176,
        179.0eV,
        2.2g / cm^3,
        Dict(
            1 => 0.0305,
            6 => 0.2107,
            7 => 0.0721,
            8 => 0.1632,
            35 => 0.2228,
            47 => 0.3007,
        ),
    ),
    emulsionstandard = Material(
        "Photographic Emulsion (Standard Nuclear)",
        0.45453,
        331.0eV,
        3.815g / cm^3,
        Dict(
            1 => 0.0141,
            6 => 0.072261,
            7 => 0.01932,
            8 => 0.066101,
            16 => 0.00189,
            35 => 0.349104,
            47 => 0.474105,
            53 => 0.00312,
        ),
    ),
    vinyltoluene = Material(
        "Plastic Scintillator, Vinyltoluene",
        0.54141,
        64.7eV,
        1.032g / cm^3,
        Dict(1 => 0.085, 6 => 0.915),
    ),
    polyethylene = Material(
        "Polyethylene",
        0.57033,
        57.4eV,
        0.93g / cm^3,
        Dict(1 => 0.143716, 6 => 0.856284),
    ),
    mylar = Material(
        "Polyethylene Terephthalate, (Mylar)",
        0.52037,
        78.7eV,
        1.38g / cm^3,
        Dict(1 => 0.04196, 6 => 0.625016, 8 => 0.333024),
    ),
    PMMA = Material(
        "Polymethyl Methacrylate",
        0.53937,
        74.0eV,
        1.19g / cm^3,
        Dict(1 => 0.080541, 6 => 0.599846, 8 => 0.319613),
    ),
    polystyrene = Material(
        "Polystyrene",
        0.53768,
        68.7eV,
        1.06g / cm^3,
        Dict(1 => 0.077421, 6 => 0.922579),
    ),
    teflon = Material(
        "Polytetrafluoroethylene, (Teflon)",
        0.47993,
        99.1eV,
        2.25g / cm^3,
        Dict(6 => 0.240183, 9 => 0.759818),
    ),
    PVC = Material(
        "Polyvinyl Chloride",
        0.51201,
        108.2eV,
        1.406g / cm^3,
        Dict(1 => 0.048382, 6 => 0.384361, 17 => 0.567257),
    ),
    nylonfilm = Material(
        "Radiochromic Dye Film, Nylon Base",
        0.54987,
        64.5eV,
        1.08g / cm^3,
        Dict(1 => 0.101996, 6 => 0.654396, 7 => 0.098915, 8 => 0.144693),
    ),
    testis = Material(
        "Testis (ICRU-44)",
        0.552,
        74.7eV,
        1.04g / cm^3,
        Dict(
            1 => 0.106,
            6 => 0.099,
            7 => 0.02,
            8 => 0.766,
            11 => 0.002,
            15 => 0.001,
            16 => 0.002,
            17 => 0.002,
            19 => 0.002,
        ),
    ),
    softtissue = Material(
        "Tissue, Soft (ICRU-44)",
        0.54996,
        74.7eV,
        1.06g / cm^3,
        Dict(
            1 => 0.102,
            6 => 0.143,
            7 => 0.034,
            8 => 0.708,
            11 => 0.002,
            15 => 0.003,
            16 => 0.003,
            17 => 0.002,
            19 => 0.003,
        ),
    ),
    softtissue4 = Material(
        "Tissue, Soft (ICRU Four-Component)",
        0.54975,
        74.9eV,
        1.0g / cm^3,
        Dict(1 => 0.101174, 6 => 0.111, 7 => 0.026, 8 => 0.761826),
    ),
    methane = Material(
        "Tissue-Equivalent Gas, Methane Based",
        0.54992,
        61.2eV,
        0.001064g / cm^3,
        Dict(1 => 0.101873, 6 => 0.456177, 7 => 0.035172, 8 => 0.406778),
    ),
    propane = Material(
        "Tissue-Equivalent Gas, Propane Based",
        0.55027,
        59.5eV,
        0.001826g / cm^3,
        Dict(1 => 0.102676, 6 => 0.568937, 7 => 0.035022, 8 => 0.293365),
    ),
    water = Material(
        "Water, Liquid",
        0.55508,
        75.0eV,
        1.0g / cm^3,
        Dict(1 => 0.111898, 8 => 0.888102),
    ),
)


text:
include("../utility.jl")
using Distributions
using Turing
using Test

xnoreturn = [1.5 2.0]

@model noreturn(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt.(s))
  for i in 1:length(x)
    x[i] ~ Normal(m, sqrt.(s))
  end
end

chain = sample(noreturn(xnoreturn), HMC(3000, 0.15, 6))

check_numerical(chain, [:s, :m], [49/24, 7/6])


text:
module MyMatching
include("deferred_acceptance.jl")
export my_deferred_acceptance
end


text:
using Yao, YaoBlocks.Optimise
using YaoTensorNetwork
using Test
using YaoExtensions
using BitBasis

@testset "generate" begin
    @test seperatecontrol() == false
    for sep in [true, false]
        seperatecontrol(sep)
        gb = GraphBuilder(4)
        attach_gate!(gb, put(4,3=>X))
        @test gb.lines == [1,2,5,4]
        attach_gate!(gb, control(4,2,3=>X))
        @test gb.lines == [1,6,7+sep,4]
        finish!(gb, [0,0,1,0])
        @test gb.lines == [-1,-1,-1,-1]
        @test ntensors(gb) == 10 + sep
        @test contract(generate_eingraph(gb)) ≈ 1
    end
end

@testset "circuit contract" begin
    nbit = 2
    # single qubit gates
    c = chain(put(2, 1=>Rx(0.5)), put(2,2=>Rx(1.5)))
    gb = circuit2tn(c)
    @test contract(gb) ≈ select(zero_state(nbit) |> c, bit"00").state[]

    # add cz
    c = chain(put(2, 1=>Rx(0.5)), put(2,2=>Rx(1.5)), control(2, 2, 1=>Z))
    for config in basis(BitStr64{2})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(nbit) |> c, config).state[]
    end

    # add cnot
    c = chain(put(2, 1=>Rx(0.5)), put(2,2=>Rx(1.5)), control(2, 2, 1=>X))
    for config in basis(BitStr64{2})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(nbit) |> c, config).state[]
    end

    # add swap
    c = chain(put(3, 1=>Rx(0.5)), put(3,2=>Rx(1.5)), put(3,3=>Rx(2.5)), put(3, (1,3)=>SWAP))
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # put rot kron
    c = chain(put(3, (2,3)=>rot(kron(X,Y), 0.5)), put(3,2=>Rx(1.5)))
    seperatecontrol(true)
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test maximum(ndims, gb.tensors) == 3
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # add control swap
    c = chain(put(3, 1=>Rx(0.5)), put(3,2=>Rx(1.5)), put(3,3=>Rx(2.5)), control(3, 2, (1,3)=>SWAP))
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # add i-cnot
    c = chain(put(3, 1=>Rx(0.5)), put(3,2=>Rx(1.5)), put(3,3=>Rx(2.5)), control(3, -2, (1,3)=>SWAP))
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # add multi-control
    c = chain(put(3, 1=>Rx(0.5)), put(3,2=>Rx(1.5)), put(3,3=>Rx(2.5)), control(3, (1,2), 3=>X))
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # add multi-control with inverse
    c = chain(put(3, 1=>Rx(0.5)), put(3,2=>Rx(1.5)), put(3,3=>Rx(2.5)), control(3, (1,-2), 3=>X))
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # add two qubit matblock
    u1 = rand_unitary(4)
    u2 = rand_unitary(4)
    c = chain(3,put((1,2)=>matblock(u1)), put((2,3)=>matblock(u2)))
    for config in basis(BitStr64{3})
        gb = circuit2tn(c; final_config=config)
        @test contract(gb) ≈ select(zero_state(3) |> c, config).state[]
    end

    # random circuit
    nbit = 2
    c = dispatch!(variational_circuit(nbit, 1, [1=>2]), :random)
    for config in basis(BitStr64{nbit})
        gb = circuit2tn(c, final_config=config)
        @test contract(gb) ≈ select(zero_state(nbit) |> c, config).state[]
    end

    # qft circuit
    nbit = 4
    c = qft_circuit(4)
    for sep in [true, false]
        seperatecontrol(sep)
        for config in basis(BitStr64{nbit})
            gb = circuit2tn(c, final_config=config)
            @test contract(gb) ≈ select(zero_state(nbit) |> c, config).state[]
        end
    end

    # supremacy 2d circuit
    nx, ny = 2, 2
    c = rand_supremacy2d(nx, ny, 6)
    c = simplify_blocktypes(c)
    deleteat!(c[2].blocks, 1)
    for sep in [true, false]
        seperatecontrol(sep)
        for config in basis(BitStr64{nbit})
            gb = circuit2tn(c, final_config=config)
            l = contract(gb)
            r = select(zero_state(nbit) |> c, config).state[]
            @test isapprox(l, r; atol=1e-8)
        end
    end

    nx, ny = 2, 2
    nbit = nx*ny
    c = variational_circuit(nbit, 1, pair_square(nx, ny; periodic=false))
    dispatch!(c, :random)
    for sep in [true, false]
        seperatecontrol(sep)
        for config in basis(BitStr64{nbit})
            gb = circuit2tn(c, final_config=config)
            #@show contract(gb)
            @test contract(gb) ≈ select(zero_state(nbit) |> c, config).state[]
        end
    end
end

@testset "saveload" begin
    nbit = 2
    c = dispatch!(variational_circuit(nbit, 1, [1=>2]), :random)
    gb = circuit2tn(c, final_config=bit"01")
    dump_graph("_test", gb)
    gb2 = load_graph(eltype(gb),"_test")
    @test gb ≈ gb2
end


text:
"""
Define abstract type for the distrubance, each subtype should define the method
genEvent(d::Event,Nlon,Nlat,Ntime)
"""
abstract Event
export Event
outevent!(ev::Event,x)=x

"""
This function generates a 3D Array of Ntime, Nlat, Nlon. It should `= 0` for all points that are not affected by an event and `> 0`
for affected grid cells. In the simplest case the resulting array is binary (either 0 or 1)
"""
genEvent(d::Event,Ntime,Nlat,Nlon)=error("The method genEvent is not defined for $(typeof(d))")

"""
n simple cubic disturbances of relative sizes sx,sy,sz (0..1) for each event with the centers at px,py,pz (0..1) for each event
The positions and sizes are relative to the cube size
"""
type CubeEvent <: Event
    n::Int
    sx::Vector{Float64}
    sy::Vector{Float64}
    sz::Vector{Float64}
    px::Vector{Float64}
    py::Vector{Float64}
    pz::Vector{Float64}
end
export CubeEvent
export genEvent
function genEvent(d::CubeEvent,Ntime,Nlat,Nlon)
    sx=d.sx*Nlon/2-0.5;px=d.px*(Nlon+1)
    sy=d.sy*Nlat/2-0.5;py=d.py*(Nlat+1)
    sz=d.sz*Ntime/2-0.5;pz=d.pz*(Ntime+1)
    # First check that disturbance fits into the cube
    a = zeros(Float64,Ntime,Nlat,Nlon)
    for i=1:d.n
        a[max(1,round(Int,pz[i]-sz[i])):min(Ntime,round(Int,pz[i]+sz[i])),max(1,round(Int,py[i]-sy[i])):min(Nlat,round(Int,py[i]+sy[i])),max(1,round(Int,px[i]-sx[i])):min(Nlon,round(Int,px[i]+sx[i]))]+=1
    end
    a
end
# Some convenience Constructors to create centered disturbances of a certain size
CubeEvent(s::AbstractFloat) = CubeEvent(1,[s],[s],[s],[0.5],[0.5],[0.5])

function CubeEvent(n::Integer,sx::Number,sy::Number,sz::Number)
    sx=fill(convert(Float64,sx),n)
    sy=fill(convert(Float64,sy),n)
    sz=fill(convert(Float64,sz),n)
    px=Float64[i/(n+1) for i=1:n]
    py=Float64[i/(n+1) for i=1:n]
    pz=Float64[i/(n+1) for i=1:n]
    CubeEvent(n,sx,sy,sz,px,py,pz)
end

"""
simple local distrubance that covers a single longitude-latitude point (0..1,0..1) over the time span s (0..1) starting at time t (0..1)
"""
type LocalEvent <: Event
    xlon::Float64
    xlat::Float64
    s   ::Float64
    t   ::Float64
end
export LocalEvent
function genEvent(d::LocalEvent,Ntime,Nlat,Nlon)
    sx = round(Int,d.xlon*Nlon)
    sy = round(Int,d.xlat*Nlat)
    tstart = round(Int,d.s*Ntime+0.5)
    tend   = tstart+round(Int,d.t*Ntime+0.5)-1
    a = zeros(Float64,Ntime,Nlat,Nlon)
    a[tstart:tend,sx,sy]=1
    a
end

"Type for an empty Disturbance"
type EmptyEvent <: Event
end
export EmptyEvent
genEvent(d::EmptyEvent,Ntime,Nlat,Nlon)=zeros(Int,Ntime,Nlat,Nlon)

"Gaussian disturbance, no sharp edegs. This is constructed by supplying a center px,py,pz and bandwidth sx,sy,sz"
type GaussianEvent
    sx::Float64
    sy::Float64
    sz::Float64
    px::Float64
    py::Float64
    pz::Float64
end
export GaussianEvent
GaussianEvent(s::Number)=GaussianEvent(s,s,s,0.5,0.5,0.5)
function genEvent(d::GaussianEvent,Ntime,Nlat,Nlon)
    sx=d.sx*Nlon/2-0.5;px=d.px*(Nlon+1)
    sy=d.sy*Nlat/2-0.5;py=d.py*(Nlat+1)
    sz=d.sz*Ntime/2-0.5;pz=d.pz*(Ntime+1)
    a=[exp((i-px)^2/sx^2 + (j-py)^2/sy^2) + (k-pz)^2/sz^2 for k=1:Ntime,j=1:Nlat,i=1:Nlon]
end

"""
This is a wrapper type that generates trend-like events in time. It wraps a given event `ev` which is usually a step
function and integrates it in the time domain to obtain slowly changing behaviour.
"""
type TrendEvent{T<:Event} <: Event
    ev::T
    persist::Bool
end
export TrendEvent
function genEvent(d::TrendEvent,Ntime,Nlat,Nlon)
    #First generate innner event
    a = genEvent(d.ev,Ntime,Nlat,Nlon)
    #Integrate over time
    d.persist ? cumsum!(a,a,1) : cumsumcontinous!(a)
    #If trend increase does not persist, set event to 0 after extreme
    #And rescale
    scale!(a,1./maximum(a))
    a
end
function outevent!(ev::TrendEvent,x)
    for ilon=1:size(x,3), ilat=1:size(x,2), itime=size(x,1):-1:2
        x[itime,ilat,ilon]=x[itime,ilat,ilon]-x[itime-1,ilat,ilon]
        x[itime,ilat,ilon]<0.0 && (x[itime,ilat,ilon]=0.0)
    end
    scale!(x,1/maximum(x))
    x
end

function cumsumcontinous!(a)
    for x in 1:size(a,3), y in 1:size(a,2), z in 2:size(a,1)
        if a[z,y,x]>0.0
            a[z,y,x]+=a[z-1,y,x]
        end
    end
end

"This generates an event of spatial size sx,sy at px,py that starts at time os and lasts until the end of the time series"
type OnsetEvent <: Event
    sx::Float64
    sy::Float64
    px::Float64
    py::Float64
    os::Float64
end
export OnsetEvent
function genEvent(d::OnsetEvent,Ntime,Nlat,Nlon)
    sx=d.sx*Nlon/2-0.5;px=d.px*(Nlon+1)
    sy=d.sy*Nlat/2-0.5;py=d.py*(Nlat+1)
    pz=d.os*(Ntime+1)
    # FIrst check that disturbance fits into the cube
    a = zeros(Float64,Ntime,Nlat,Nlon)
    a[round(Int,pz):Ntime,round(Int,py-sy):round(Int,py+sy),round(Int,px-sx):round(Int,px+sx)]=1
    a
end

"RandomWalkEvent, generates N events. The event is generated by doing a Random Walk in 3D crossing starting at px,py,pz. with step probabilities px,py,pz in the respective directions"
type RandomWalkEvent <: Event
    n::Int
    sx::Vector{Float64}
    sy::Vector{Float64}
    sz::Vector{Float64}
    px::Vector{Float64}
    py::Vector{Float64}
    pz::Vector{Float64}
end
export RandomWalkEvent
function RandomWalkEvent(n::Integer,sx::Number,sy::Number,sz::Number)
    sx=fill(convert(Float64,sx),n);sy=fill(convert(Float64,sy),n);sz=fill(convert(Float64,sz),n);px=Float64[i/(n+1) for i=1:n];py=Float64[i/(n+1) for i=1:n];pz=Float64[i/(n+1) for i=1:n];
    RandomWalkEvent(n,sx,sy,sz,px,py,pz)
end
function genEvent(d::RandomWalkEvent,Ntime,Nlat,Nlon)
    sx=d.sx*Nlon;px=round(Int,d.px*(Nlon+1))
    sy=d.sy*Nlat;py=round(Int,d.py*(Nlat+1))
    sz=d.sz*Ntime;pz=round(Int,d.pz*(Ntime+1))
    a = zeros(Float64,Ntime,Nlat,Nlon)
    for ievent=1:d.n
        walkweights=[sx[ievent]^2,sx[ievent]^2,sy[ievent]^2,sy[ievent]^2,sz[ievent]^2,sz[ievent]^2]
        scale!(walkweights,1/sum(walkweights))
        distri=Categorical(walkweights)
        Ntot=round(Int,sx[ievent]*sy[ievent]*sz[ievent])
        x=px[ievent];y=py[ievent];z=pz[ievent]
        a[z,y,x]=1 #Set starting point to 1
        itot=1
        while itot<Ntot
            direc = rand(distri)
            (direc == 1) && (x=max(1,x-1))
            (direc == 2) && (x=min(Nlon,x+1))
            (direc == 3) && (y=max(1,y-1))
            (direc == 4) && (y=min(Nlat,y+1))
            (direc == 5) && (z=max(1,z-1))
            (direc == 6) && (z=min(Ntime,z+1))
            if a[z,y,x]==0
                a[z,y,x]=1
                itot=itot+1
            end
        end
    end
    a
end

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
export NKqLandscape

# Construct an NKq landscape, see:
# Newman, M., & Engelhardt, R. (1998). Effects of selective neutrality on the
# evolution of molecular species. Proceedings of the Royal Society of London.
# Series B: Biological Sciences, 265(1403)
type NKqLandscape <: Landscape
  n::Int64
  k::Int64
  q::Int64
  a::Int64
  links::AlleleLinks
  contribs::AlleleContribs
end

function NKqLandscape(n::Int64, k::Int64, q::Int64; a::Int64=2, near::Bool=false)
  if a != 2
    # TODO:  implement a values other than 2
    error("only a == 2 is implemented at this time")
  end
  links = makelinks(n, k, near)
  contribs = newAlleleContribs(n)
  return NKqLandscape(n, k, q, a, links, contribs)
end

show(io::Base.IO, l::NKqLandscape) = print(io, "NKqLandscape($(l.n), $(l.k), $(l.q))")


text:
import Compat.Pkg

@testset "MathProgBase tests" begin
    mathprogbase_test = joinpath(Pkg.dir("MathProgBase"), "test")

    include(joinpath(mathprogbase_test, "linprog.jl"))
    linprogtest(GLPKSolverLP())

    include(joinpath(mathprogbase_test, "linproginterface.jl"))
    linprogsolvertest(GLPKSolverLP())

    include(joinpath(mathprogbase_test, "mixintprog.jl"))
    mixintprogtest(GLPKSolverMIP())

    include(joinpath(mathprogbase_test,"conicinterface.jl"))
    coniclineartest(GLPKSolverLP())

    solver = GLPKSolverMIP(msg_lev=GLPK.MSG_ALL)

    # Silent should override GLPK.MSG_ALL
    MathProgBase.setparameters!(solver, Silent=true, TimeLimit=100.0)
    coniclineartest(solver)

end


text:
using LinearAlgebra, Quaternions, ReferenceFrameRotations
# Check if a given vector is three-dimensional
function check_vector_3(v)
    if size(v) != (1,3) && size(v) != (3, 1) && size(v) != (3,)
        println(size(v))
        println("\nIncorrect vector3. You must introduce an R3 vector.")
        return false
    end
end
# Format vector3 to be vertically oriented
function format_vector_3(v)
    if size(v) == (1,3)
        v = v'
    end
    return v
end

function to_2d(vector3)
    # Check if the given vector is actually on three dimensions
    if check_vector_3(vector3) == false return end
    # If the vector is given horizontal, we transpose it
    vector3 = format_vector_3(vector3)

    # Defined matrix to convert R3 vector to R2 vector
    # (we find this matrix on the .pdf attached file)
    axonometry_matrix = [-0.5*cosd(42)    cosd(7)     0;
                       -0.5*sind(42)    -sind(7)    1]

    #=axonometry_matrix = [-cosd(30)    cosd(30)     0;
                        -sind(30)    -sind(30)    1]=#

    # Return the result of the multiplication
    return axonometry_matrix * vector3
end

function rotate_phi_z(X, phi)
    # Check if the given vector is actually on three dimensions
    if check_vector_3(X) == false
        return
    end
    # If the vector is given horizontal, we transpose it
    X = format_vector_3(X)

    # Define rotation matrix on Z axis
    R = [ 0 -1 0;
          1 0 0;
          0 0 0 ]

    # Apply the Rodrigues formula
    result = I + sind(phi) * R + (1 - cosd(phi)) * R^2

    # Calculate result using the rotation matrix
    T = result * X
    return T

end

function rotate_phi(X, phi, axis)
    # Check if the given vector is actually on three dimensions
    if check_vector_3(X) == false
        return
    end
    # If the vector is given horizontal, we transpose it
    X = format_vector_3(X)

    axis = normalize(axis)

    # Define rotation matrix on Z axis
    R = [ 0      -axis[3]    axis[2];
         axis[3]    0      -axis[1];
         -axis[2]  axis[1]      0 ]

    # Apply the Rodrigues formula
    result = I + sind(phi) * R + (1 - cosd(phi)) * R^2

    # Calculate result using the rotation matrix
    T = result * X
    return T

end

# To check if the given vector is of dimension 3
function check_vector_3(v)
    if size(v) != (1,3) && size(v) != (3, 1) && size(v) != (3,)
        print("Your vector is size ", size(v))
        println("\nIncorrect vector3. You must introduce an R3 vector.")
        return false
    end
end

# To check if the given object is a Quaternion
function check_quaternion(q)
    if typeof(q) == Quaternions.Quaternion{Float64} || typeof(q) == Quaternions.Quaternion{Int64}
        return true
    end
    println("You must introduce a valid Quaternion! Using Quaternion()")
    return false
end

# To check if the given matrix is of dimensions 3x3
function check_matrix_3(R)
    if size(R) != (3,3)
        print("Your matrix is size ", size(R))
        println("\nIncorrect matrix. You must introduce a 3x3 matrix.")
        return false
    end
end

function check_euler_angleaxis(E)
    if typeof(E) == EulerAngleAxis{Float64} || typeof(E) == EulerAngleAxis{Int64}
        return true
    end
    println("Incorrect EulerAngleAxis. You must introduce an EulerAngleAxis()")
    return false
end
# Axis/Angle to Quaternion
function axis_angle_to_quat(E)

    if check_euler_angleaxis(E) == false
        return
    end

    # We get the axis from our EulerAngleAxis
    V = Vector{Float64}(E.v)
    # We get the angle from our EulerAngleAxis
    phi = E.a

    # Use qrotation function to create a Quaternion with the given axis and angle
    qR = qrotation(V, phi)

    return qR

end

function axis_angle_to_mat(E)

    if check_euler_angleaxis(E) == false
        return
    end
    # We get the axis from our EulerAngleAxis
    V = normalize(E.v)
    # We get the angle from our EulerAngleAxis
    phi = E.a

    # We use Rodrigues Formula to find the Rotation matrix
    # Define rotation matrix on Z axis
    R = [ 0     -V[3]   V[2];
          V[3]   0      -V[1];
         -V[2]   V[1]      0 ]

    # Apply the Rodrigues formula
    result = I + sin(phi) * R + (1 - cos(phi)) * R^2

    return result
end

function mat_to_axis_angle(R)
    if check_matrix_3(R) == false
        return
    end

# Find angle
phi = acos((tr(R) -1) / 2)

# Find axis matrix
V = (R - R') / (2 * sin(phi))

# Get axis values from the axis matrix into a vector
axis = [V[3,2]; V[1,3]; V[2,1]]

return EulerAngleAxis(phi, axis)
end

function quat_to_axis_angle(q)
    if check_quaternion(q) == false
        return
    end

    # Get the angle
    angle = 2 * atan(sqrt(q.v1^2 + q.v2^2 + q.v3^2) / q.s)

    # Normalize the Quaternion
    q = normalize(q)

    # We do the same process as in the
    s = sin(angle/2)

    # If the s is not zero, we divide. If it is zero, we can't divide.
    # Therefore, we set axis as the [1 0 0] vector
    s != 0 ? axis = [q.v1; q.v2; q.v3] / s : axis = [1.0; 0.0; 0.0]

    return EulerAngleAxis(angle, axis)
end

function quat_to_mat(q)
    if check_quaternion(q) == false
        return
    end
    # We apply the matrix to Quaternion transformation
    R = [q.s^2 + q.v1^2 - q.v2^2 - q.v3^2    2*q.v1*q.v2 - 2*q.s*q.v3            2*q.v1*q.v3 + 2*q.s*q.v2;
     2*q.v1*q.v2 + 2*q.s*q.v3                q.s^2 - q.v1^2 + q.v2^2 - q.v3^2    2*q.v2*q.v3 - 2*q.s*q.v1;
     2*q.v1*q.v3 - 2*q.s*q.v2                2*q.v2*q.v3 + 2*q.s*q.v1            q.s^2 - q.v1^2 - q.v2^2 + q.v3^2]

     return R

end

function mat_to_quat(R)
    if check_matrix_3(R) == false
        return
    end

    E = mat_to_axis_angle(R)

    q = axis_angle_to_quat(E)

    return q
end

function rescale_transalte_vector(V, pos, scale)
    # Get Vector without considering origin
    # Change Vector position using pos
    # Position is on world coordinates (always positive)

    #Check vector and position are on the correct format
    V = format_vector_3(V)
    pos = format_vector_3(pos)

    # Rescale vector
    V = V*scale

    # Change vector position
    V = V+pos

    return V

end


text:
# for geom in (:Point, :MultiPoint, :LineString, :MultiLineString, :LinearRing, :Polygon, :MultiPolygon)
#     @eval Base.copy(obj::$geom) = ($geom)(GEOSGeom_clone(obj.ptr))
# end

GeoInterface.coordinates(obj::Point) =
    isEmpty(obj.ptr) ? Float64[] : getCoordinates(getCoordSeq(obj.ptr), 1)
GeoInterface.coordinates(obj::LineString) = getCoordinates(getCoordSeq(obj.ptr))
GeoInterface.coordinates(obj::LinearRing) = getCoordinates(getCoordSeq(obj.ptr))

function GeoInterface.coordinates(polygon::Polygon)
    exterior = getCoordinates(getCoordSeq(exteriorRing(polygon.ptr)))
    interiors = [getCoordinates(getCoordSeq(ring)) for ring in interiorRings(polygon.ptr)]
    if length(interiors) == 0
        return Vector{Vector{Float64}}[exterior]
    else
        return [Vector{Vector{Float64}}[exterior]; interiors]
    end
end

GeoInterface.coordinates(multipoint::MultiPoint) = Vector{Float64}[
    getCoordinates(getCoordSeq(geom), 1) for geom in getGeometries(multipoint.ptr)
]
GeoInterface.coordinates(multiline::MultiLineString) = Vector{Vector{Float64}}[
    getCoordinates(getCoordSeq(geom)) for geom in getGeometries(multiline.ptr)
]
function GeoInterface.coordinates(multipolygon::MultiPolygon)
    geometries = getGeometries(multipolygon.ptr)
    coords = Array{Vector{Vector{Vector{Float64}}}}(undef, length(geometries))
    for (i, geom) in enumerate(getGeometries(multipolygon.ptr))
        exterior = getCoordinates(getCoordSeq(exteriorRing(geom)))
        interiors = [getCoordinates(getCoordSeq(ring)) for ring in interiorRings(geom)]
        coords[i] = [Vector{Vector{Float64}}[exterior]; interiors]
    end
    coords
end

function GeoInterface.geometries(obj::GeometryCollection)
    collection = GeoInterface.AbstractGeometry[]
    sizehint!(collection, numGeometries(obj.ptr))
    for geom in getGeometries(obj.ptr)
        if geomTypeId(geom) == GEOS_POINT
            push!(collection, Point(cloneGeom(geom)))
        elseif geomTypeId(geom) == GEOS_LINESTRING
            push!(collection, LineString(cloneGeom(geom)))
        elseif geomTypeId(geom) == GEOS_LINEARRING
            push!(collection, LinearRing(cloneGeom(geom)))
        elseif geomTypeId(geom) == GEOS_POLYGON
            push!(collection, Polygon(cloneGeom(geom)))
        elseif geomTypeId(geom) == GEOS_MULTIPOINT
            push!(collection, MultiPoint(cloneGeom(geom)))
        elseif geomTypeId(geom) == GEOS_MULTILINESTRING
            push!(collection, MultiLineString(cloneGeom(geom)))
        elseif geomTypeId(geom) == GEOS_MULTIPOLYGON
            push!(collection, MultiPolygon(cloneGeom(geom)))
        else
            @assert geomTypeId(geom) == GEOS_GEOMETRYCOLLECTION
            push!(collection, GeometryCollection(cloneGeom(geom)))
        end
    end
    collection
end


text:
#!/usr/bin/env julia

using Luxor, Colors

fname = "/tmp/Tiler-test1.pdf"
pagewidth, pageheight = 600, 900
Drawing(pagewidth, pageheight, fname)
origin() # move 0/0 to center
background("ivory")
setopacity(0.9)
setline(0.3)

# width, height, nrows, ncols, margin
pagetiles = Tiler(pagewidth, pageheight, 8, 2, margin=20)

setopacity(0.5)

for (pos, n) in pagetiles
  randomhue()
  ellipse(pos, pagetiles.tilewidth, pagetiles.tileheight, :fill)
end

pagetiles = Tiler(200, 300, 4, 5, margin=20)
for (pos, n) in pagetiles
  randomhue()
  ellipse(pos, pagetiles.tilewidth, pagetiles.tileheight, :fill)
end

finish()
preview()