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

module Iterators

import Base: start, done, next, isempty, length, size, eltype, iteratorsize, iteratoreltype, indices, ndims

using Base: tuple_type_cons, SizeUnknown, HasLength, HasShape, IsInfinite, EltypeUnknown, HasEltype, OneTo, @propagate_inbounds

export enumerate, zip, rest, countfrom, take, drop, cycle, repeated, product, flatten, partition

_min_length(a, b, ::IsInfinite, ::IsInfinite) = min(length(a),length(b)) # inherit behaviour, error
_min_length(a, b, A, ::IsInfinite) = length(a)
_min_length(a, b, ::IsInfinite, B) = length(b)
_min_length(a, b, A, B) = min(length(a),length(b))

_diff_length(a, b, A, ::IsInfinite) = 0
_diff_length(a, b, ::IsInfinite, ::IsInfinite) = 0
_diff_length(a, b, ::IsInfinite, B) = length(a) # inherit behaviour, error
_diff_length(a, b, A, B) = max(length(a)-length(b), 0)

and_iteratorsize(isz::T, ::T) where {T} = isz
and_iteratorsize(::HasLength, ::HasShape) = HasLength()
and_iteratorsize(::HasShape, ::HasLength) = HasLength()
and_iteratorsize(a, b) = SizeUnknown()

and_iteratoreltype(iel::T, ::T) where {T} = iel
and_iteratoreltype(a, b) = EltypeUnknown()

# enumerate

struct Enumerate{I}
    itr::I
end

"""
    enumerate(iter)

An iterator that yields `(i, x)` where `i` is a counter starting at 1,
and `x` is the `i`th value from the given iterator. It's useful when
you need not only the values `x` over which you are iterating, but
also the number of iterations so far. Note that `i` may not be valid
for indexing `iter`; it's also possible that `x != iter[i]`, if `iter`
has indices that do not start at 1. See the `enumerate(IndexLinear(),
iter)` method if you want to ensure that `i` is an index.

# Examples
```jldoctest
julia> a = ["a", "b", "c"];

julia> for (index, value) in enumerate(a)
           println("\$index \$value")
       end
1 a
2 b
3 c
```
"""
enumerate(iter) = Enumerate(iter)

length(e::Enumerate) = length(e.itr)
size(e::Enumerate) = size(e.itr)
@inline start(e::Enumerate) = (1, start(e.itr))
@inline function next(e::Enumerate, state)
    n = next(e.itr,state[2])
    (state[1],n[1]), (state[1]+1,n[2])
end
@inline done(e::Enumerate, state) = done(e.itr, state[2])

eltype(::Type{Enumerate{I}}) where {I} = Tuple{Int, eltype(I)}

iteratorsize(::Type{Enumerate{I}}) where {I} = iteratorsize(I)
iteratoreltype(::Type{Enumerate{I}}) where {I} = iteratoreltype(I)

struct IndexValue{I,A<:AbstractArray}
    data::A
    itr::I
end

"""
    enumerate(IndexLinear(), A)
    enumerate(IndexCartesian(), A)
    enumerate(IndexStyle(A), A)

An iterator that accesses each element of the array `A`, returning
`(i, x)`, where `i` is the index for the element and `x = A[i]`.  This
is similar to `enumerate(A)`, except `i` will always be a valid index
for `A`.

Specifying `IndexLinear()` ensures that `i` will be an integer;
specifying `IndexCartesian()` ensures that `i` will be a
`CartesianIndex`; specifying `IndexStyle(A)` chooses whichever has
been defined as the native indexing style for array `A`.

# Examples
```jldoctest
julia> A = ["a" "d"; "b" "e"; "c" "f"];

julia> for (index, value) in enumerate(IndexStyle(A), A)
           println("\$index \$value")
       end
1 a
2 b
3 c
4 d
5 e
6 f

julia> S = view(A, 1:2, :);

julia> for (index, value) in enumerate(IndexStyle(S), S)
           println("\$index \$value")
       end
CartesianIndex{2}((1, 1)) a
CartesianIndex{2}((2, 1)) b
CartesianIndex{2}((1, 2)) d
CartesianIndex{2}((2, 2)) e
```

Note that `enumerate(A)` returns `i` as a *counter* (always starting
at 1), whereas `enumerate(IndexLinear(), A)` returns `i` as an *index*
(starting at the first linear index of `A`, which may or may not be
1).

See also: [`IndexStyle`](@ref), [`indices`](@ref).
"""
enumerate(::IndexLinear,    A::AbstractArray) = IndexValue(A, linearindices(A))
enumerate(::IndexCartesian, A::AbstractArray) = IndexValue(A, CartesianRange(indices(A)))

length(v::IndexValue)  = length(v.itr)
indices(v::IndexValue) = indices(v.itr)
size(v::IndexValue)    = size(v.itr)
@inline start(v::IndexValue) = start(v.itr)
@propagate_inbounds function next(v::IndexValue, state)
    indx, n = next(v.itr, state)
    item = v.data[indx]
    (indx, item), n
end
@inline done(v::IndexValue, state) = done(v.itr, state)

eltype(::Type{IndexValue{I,A}}) where {I,A} = Tuple{eltype(I), eltype(A)}

iteratorsize(::Type{IndexValue{I}}) where {I} = iteratorsize(I)
iteratoreltype(::Type{IndexValue{I}}) where {I} = iteratoreltype(I)

# zip

abstract type AbstractZipIterator end

zip_iteratorsize(a, b) = and_iteratorsize(a,b) # as `and_iteratorsize` but inherit `Union{HasLength,IsInfinite}` of the shorter iterator
zip_iteratorsize(::HasLength, ::IsInfinite) = HasLength()
zip_iteratorsize(::HasShape, ::IsInfinite) = HasLength()
zip_iteratorsize(a::IsInfinite, b) = zip_iteratorsize(b,a)
zip_iteratorsize(a::IsInfinite, b::IsInfinite) = IsInfinite()

struct Zip1{I} <: AbstractZipIterator
    a::I
end
zip(a) = Zip1(a)
length(z::Zip1) = length(z.a)
size(z::Zip1) = size(z.a)
indices(z::Zip1) = indices(z.a)
eltype(::Type{Zip1{I}}) where {I} = Tuple{eltype(I)}
@inline start(z::Zip1) = start(z.a)
@inline function next(z::Zip1, st)
    n = next(z.a,st)
    return ((n[1],), n[2])
end
@inline done(z::Zip1, st) = done(z.a,st)

iteratorsize(::Type{Zip1{I}}) where {I} = iteratorsize(I)
iteratoreltype(::Type{Zip1{I}}) where {I} = iteratoreltype(I)

struct Zip2{I1, I2} <: AbstractZipIterator
    a::I1
    b::I2
end
zip(a, b) = Zip2(a, b)
length(z::Zip2) = _min_length(z.a, z.b, iteratorsize(z.a), iteratorsize(z.b))
size(z::Zip2) = promote_shape(size(z.a), size(z.b))
indices(z::Zip2) = promote_shape(indices(z.a), indices(z.b))
eltype(::Type{Zip2{I1,I2}}) where {I1,I2} = Tuple{eltype(I1), eltype(I2)}
@inline start(z::Zip2) = (start(z.a), start(z.b))
@inline function next(z::Zip2, st)
    n1 = next(z.a,st[1])
    n2 = next(z.b,st[2])
    return ((n1[1], n2[1]), (n1[2], n2[2]))
end
@inline done(z::Zip2, st) = done(z.a,st[1]) | done(z.b,st[2])

iteratorsize(::Type{Zip2{I1,I2}}) where {I1,I2} = zip_iteratorsize(iteratorsize(I1),iteratorsize(I2))
iteratoreltype(::Type{Zip2{I1,I2}}) where {I1,I2} = and_iteratoreltype(iteratoreltype(I1),iteratoreltype(I2))

struct Zip{I, Z<:AbstractZipIterator} <: AbstractZipIterator
    a::I
    z::Z
end

"""
    zip(iters...)

For a set of iterable objects, returns an iterable of tuples, where the `i`th tuple contains
the `i`th component of each input iterable.

Note that [`zip`](@ref) is its own inverse: `collect(zip(zip(a...)...)) == collect(a)`.

# Examples
```jldoctest
julia> a = 1:5
1:5

julia> b = ["e","d","b","c","a"]
5-element Array{String,1}:
 "e"
 "d"
 "b"
 "c"
 "a"

julia> c = zip(a,b)
Base.Iterators.Zip2{UnitRange{Int64},Array{String,1}}(1:5, String["e", "d", "b", "c", "a"])

julia> length(c)
5

julia> first(c)
(1, "e")
```
"""
zip(a, b, c...) = Zip(a, zip(b, c...))
length(z::Zip) = _min_length(z.a, z.z, iteratorsize(z.a), iteratorsize(z.z))
size(z::Zip) = promote_shape(size(z.a), size(z.z))
indices(z::Zip) = promote_shape(indices(z.a), indices(z.z))
eltype(::Type{Zip{I,Z}}) where {I,Z} = tuple_type_cons(eltype(I), eltype(Z))
@inline start(z::Zip) = tuple(start(z.a), start(z.z))
@inline function next(z::Zip, st)
    n1 = next(z.a, st[1])
    n2 = next(z.z, st[2])
    (tuple(n1[1], n2[1]...), (n1[2], n2[2]))
end
@inline done(z::Zip, st) = done(z.a,st[1]) | done(z.z,st[2])

iteratorsize(::Type{Zip{I1,I2}}) where {I1,I2} = zip_iteratorsize(iteratorsize(I1),iteratorsize(I2))
iteratoreltype(::Type{Zip{I1,I2}}) where {I1,I2} = and_iteratoreltype(iteratoreltype(I1),iteratoreltype(I2))

# filter

struct Filter{F,I}
    flt::F
    itr::I
end

"""
    Iterators.filter(flt, itr)

Given a predicate function `flt` and an iterable object `itr`, return an
iterable object which upon iteration yields the elements `x` of `itr` that
satisfy `flt(x)`. The order of the original iterator is preserved.

This function is *lazy*; that is, it is guaranteed to return in ``Θ(1)`` time
and use ``Θ(1)`` additional space, and `flt` will not be called by an
invocation of `filter`. Calls to `flt` will be made when iterating over the
returned iterable object. These calls are not cached and repeated calls will be
made when reiterating.

See [`Base.filter`](@ref) for an eager implementation of filtering for arrays.

# Examples
```jldoctest
julia> f = Iterators.filter(isodd, [1, 2, 3, 4, 5])
Base.Iterators.Filter{Base.#isodd,Array{Int64,1}}(isodd, [1, 2, 3, 4, 5])

julia> foreach(println, f)
1
3
5
```
"""
filter(flt, itr) = Filter(flt, itr)

start(f::Filter) = start_filter(f.flt, f.itr)
function start_filter(pred, itr)
    s = start(itr)
    while !done(itr,s)
        v,t = next(itr,s)
        if pred(v)
            return (false, v, t)
        end
        s=t
    end
    (true,)
end

next(f::Filter, s) = advance_filter(f.flt, f.itr, s)
function advance_filter(pred, itr, st)
    _, v, s = st
    while !done(itr,s)
        w,t = next(itr,s)
        if pred(w)
            return v, (false, w, t)
        end
        s=t
    end
    v, (true, v, s)
end

done(f::Filter, s) = s[1]

eltype(::Type{Filter{F,I}}) where {F,I} = eltype(I)
iteratoreltype(::Type{Filter{F,I}}) where {F,I} = iteratoreltype(I)
iteratorsize(::Type{<:Filter}) = SizeUnknown()

# Rest -- iterate starting at the given state

struct Rest{I,S}
    itr::I
    st::S
end

"""
    rest(iter, state)

An iterator that yields the same elements as `iter`, but starting at the given `state`.

# Examples
```jldoctest
julia> collect(Iterators.rest([1,2,3,4], 2))
3-element Array{Any,1}:
 2
 3
 4
```
"""
rest(itr,state) = Rest(itr,state)

start(i::Rest) = i.st
next(i::Rest, st) = next(i.itr, st)
done(i::Rest, st) = done(i.itr, st)

eltype(::Type{Rest{I}}) where {I} = eltype(I)
iteratoreltype(::Type{Rest{I,S}}) where {I,S} = iteratoreltype(I)
rest_iteratorsize(a) = SizeUnknown()
rest_iteratorsize(::IsInfinite) = IsInfinite()
iteratorsize(::Type{Rest{I,S}}) where {I,S} = rest_iteratorsize(iteratorsize(I))

# Count -- infinite counting

struct Count{S<:Number}
    start::S
    step::S
end

"""
    countfrom(start=1, step=1)

An iterator that counts forever, starting at `start` and incrementing by `step`.

# Examples
```jldoctest
julia> for v in Iterators.countfrom(5, 2)
           v > 10 && break
           println(v)
       end
5
7
9
```
"""
countfrom(start::Number, step::Number) = Count(promote(start, step)...)
countfrom(start::Number)               = Count(start, oneunit(start))
countfrom()                            = Count(1, 1)

eltype(::Type{Count{S}}) where {S} = S

start(it::Count) = it.start
next(it::Count, state) = (state, state + it.step)
done(it::Count, state) = false

iteratorsize(::Type{<:Count}) = IsInfinite()

# Take -- iterate through the first n elements

struct Take{I}
    xs::I
    n::Int
end

"""
    take(iter, n)

An iterator that generates at most the first `n` elements of `iter`.

# Examples
```jldoctest
julia> a = 1:2:11
1:2:11

julia> collect(a)
6-element Array{Int64,1}:
  1
  3
  5
  7
  9
 11

julia> collect(Iterators.take(a,3))
3-element Array{Int64,1}:
 1
 3
 5
```
"""
take(xs, n::Integer) = Take(xs, Int(n))
take(xs::Take, n::Integer) = Take(xs.xs, min(Int(n), xs.n))

eltype(::Type{Take{I}}) where {I} = eltype(I)
iteratoreltype(::Type{Take{I}}) where {I} = iteratoreltype(I)
take_iteratorsize(a) = HasLength()
take_iteratorsize(::SizeUnknown) = SizeUnknown()
iteratorsize(::Type{Take{I}}) where {I} = take_iteratorsize(iteratorsize(I))
length(t::Take) = _min_length(t.xs, 1:t.n, iteratorsize(t.xs), HasLength())

start(it::Take) = (it.n, start(it.xs))

function next(it::Take, state)
    n, xs_state = state
    v, xs_state = next(it.xs, xs_state)
    return v, (n - 1, xs_state)
end

function done(it::Take, state)
    n, xs_state = state
    return n <= 0 || done(it.xs, xs_state)
end

# Drop -- iterator through all but the first n elements

struct Drop{I}
    xs::I
    n::Int
end

"""
    drop(iter, n)

An iterator that generates all but the first `n` elements of `iter`.

# Examples
```jldoctest
julia> a = 1:2:11
1:2:11

julia> collect(a)
6-element Array{Int64,1}:
  1
  3
  5
  7
  9
 11

julia> collect(Iterators.drop(a,4))
2-element Array{Int64,1}:
  9
 11
```
"""
drop(xs, n::Integer) = Drop(xs, Int(n))
drop(xs::Take, n::Integer) = Take(drop(xs.xs, Int(n)), max(0, xs.n - Int(n)))
drop(xs::Drop, n::Integer) = Drop(xs.xs, Int(n) + xs.n)

eltype(::Type{Drop{I}}) where {I} = eltype(I)
iteratoreltype(::Type{Drop{I}}) where {I} = iteratoreltype(I)
drop_iteratorsize(::SizeUnknown) = SizeUnknown()
drop_iteratorsize(::Union{HasShape, HasLength}) = HasLength()
drop_iteratorsize(::IsInfinite) = IsInfinite()
iteratorsize(::Type{Drop{I}}) where {I} = drop_iteratorsize(iteratorsize(I))
length(d::Drop) = _diff_length(d.xs, 1:d.n, iteratorsize(d.xs), HasLength())

function start(it::Drop)
    xs_state = start(it.xs)
    for i in 1:it.n
        if done(it.xs, xs_state)
            break
        end

        _, xs_state = next(it.xs, xs_state)
    end
    xs_state
end

next(it::Drop, state) = next(it.xs, state)
done(it::Drop, state) = done(it.xs, state)

# Cycle an iterator forever

struct Cycle{I}
    xs::I
end

"""
    cycle(iter)

An iterator that cycles through `iter` forever.

# Examples
```jldoctest
julia> for (i, v) in enumerate(Iterators.cycle("hello"))
           print(v)
           i > 10 && break
       end
hellohelloh
```
"""
cycle(xs) = Cycle(xs)

eltype(::Type{Cycle{I}}) where {I} = eltype(I)
iteratoreltype(::Type{Cycle{I}}) where {I} = iteratoreltype(I)
iteratorsize(::Type{Cycle{I}}) where {I} = IsInfinite()

function start(it::Cycle)
    s = start(it.xs)
    return s, done(it.xs, s)
end

function next(it::Cycle, state)
    s, d = state
    if done(it.xs, s)
        s = start(it.xs)
    end
    v, s = next(it.xs, s)
    return v, (s, false)
end

done(it::Cycle, state) = state[2]

# Repeated - repeat an object infinitely many times

struct Repeated{O}
    x::O
end
repeated(x) = Repeated(x)

"""
    repeated(x[, n::Int])

An iterator that generates the value `x` forever. If `n` is specified, generates `x` that
many times (equivalent to `take(repeated(x), n)`).

# Examples
```jldoctest
julia> a = Iterators.repeated([1 2], 4);

julia> collect(a)
4-element Array{Array{Int64,2},1}:
 [1 2]
 [1 2]
 [1 2]
 [1 2]
```
"""
repeated(x, n::Integer) = take(repeated(x), Int(n))

eltype(::Type{Repeated{O}}) where {O} = O

start(it::Repeated) = nothing
next(it::Repeated, state) = (it.x, nothing)
done(it::Repeated, state) = false

iteratorsize(::Type{<:Repeated}) = IsInfinite()
iteratoreltype(::Type{<:Repeated}) = HasEltype()

# Product -- cartesian product of iterators

abstract type AbstractProdIterator end

length(p::AbstractProdIterator) = prod(size(p))
_length(p::AbstractProdIterator) = prod(map(unsafe_length, indices(p)))
size(p::AbstractProdIterator) = _prod_size(p.a, p.b, iteratorsize(p.a), iteratorsize(p.b))
indices(p::AbstractProdIterator) = _prod_indices(p.a, p.b, iteratorsize(p.a), iteratorsize(p.b))
ndims(p::AbstractProdIterator) = length(indices(p))

# generic methods to handle size of Prod* types
_prod_size(a, ::HasShape)  = size(a)
_prod_size(a, ::HasLength) = (length(a),)
_prod_size(a, A) =
    throw(ArgumentError("Cannot compute size for object of type $(typeof(a))"))
_prod_size(a, b, ::HasLength, ::HasLength)  = (length(a),  length(b))
_prod_size(a, b, ::HasLength, ::HasShape)   = (length(a),  size(b)...)
_prod_size(a, b, ::HasShape,  ::HasLength)  = (size(a)..., length(b))
_prod_size(a, b, ::HasShape,  ::HasShape)   = (size(a)..., size(b)...)
_prod_size(a, b, A, B) =
    throw(ArgumentError("Cannot construct size for objects of types $(typeof(a)) and $(typeof(b))"))

_prod_indices(a, ::HasShape)  = indices(a)
_prod_indices(a, ::HasLength) = (OneTo(length(a)),)
_prod_indices(a, A) =
    throw(ArgumentError("Cannot compute indices for object of type $(typeof(a))"))
_prod_indices(a, b, ::HasLength, ::HasLength)  = (OneTo(length(a)),  OneTo(length(b)))
_prod_indices(a, b, ::HasLength, ::HasShape)   = (OneTo(length(a)),  indices(b)...)
_prod_indices(a, b, ::HasShape,  ::HasLength)  = (indices(a)..., OneTo(length(b)))
_prod_indices(a, b, ::HasShape,  ::HasShape)   = (indices(a)..., indices(b)...)
_prod_indices(a, b, A, B) =
    throw(ArgumentError("Cannot construct indices for objects of types $(typeof(a)) and $(typeof(b))"))

# one iterator
struct Prod1{I} <: AbstractProdIterator
    a::I
end
product(a) = Prod1(a)

eltype(::Type{Prod1{I}}) where {I} = Tuple{eltype(I)}
size(p::Prod1) = _prod_size(p.a, iteratorsize(p.a))
indices(p::Prod1) = _prod_indices(p.a, iteratorsize(p.a))

@inline start(p::Prod1) = start(p.a)
@inline function next(p::Prod1, st)
    n, st = next(p.a, st)
    (n,), st
end
@inline done(p::Prod1, st) = done(p.a, st)

iteratoreltype(::Type{Prod1{I}}) where {I} = iteratoreltype(I)
iteratorsize(::Type{Prod1{I}}) where {I} = iteratorsize(I)

# two iterators
struct Prod2{I1, I2} <: AbstractProdIterator
    a::I1
    b::I2
end

"""
    product(iters...)

Returns an iterator over the product of several iterators. Each generated element is
a tuple whose `i`th element comes from the `i`th argument iterator. The first iterator
changes the fastest.

# Examples
```jldoctest
julia> collect(Iterators.product(1:2,3:5))
2×3 Array{Tuple{Int64,Int64},2}:
 (1, 3)  (1, 4)  (1, 5)
 (2, 3)  (2, 4)  (2, 5)
```
"""
product(a, b) = Prod2(a, b)

eltype(::Type{Prod2{I1,I2}}) where {I1,I2} = Tuple{eltype(I1), eltype(I2)}

iteratoreltype(::Type{Prod2{I1,I2}}) where {I1,I2} = and_iteratoreltype(iteratoreltype(I1),iteratoreltype(I2))
iteratorsize(::Type{Prod2{I1,I2}}) where {I1,I2} = prod_iteratorsize(iteratorsize(I1),iteratorsize(I2))

function start(p::AbstractProdIterator)
    s1, s2 = start(p.a), start(p.b)
    s1, s2, Nullable{eltype(p.b)}(), (done(p.a,s1) || done(p.b,s2))
end

@inline function prod_next(p, st)
    s1, s2 = st[1], st[2]
    v1, s1 = next(p.a, s1)

    nv2 = st[3]
    if isnull(nv2)
        v2, s2 = next(p.b, s2)
    else
        v2 = nv2.value
    end

    if done(p.a, s1)
        return (v1,v2), (start(p.a), s2, Nullable{eltype(nv2)}(), done(p.b,s2))
    end
    return (v1,v2), (s1, s2, Nullable(v2), false)
end

@inline next(p::Prod2, st) = prod_next(p, st)
@inline done(p::AbstractProdIterator, st) = st[4]

# n iterators
struct Prod{I1, I2<:AbstractProdIterator} <: AbstractProdIterator
    a::I1
    b::I2
end
product(a, b, c...) = Prod(a, product(b, c...))

eltype(::Type{Prod{I1,I2}}) where {I1,I2} = tuple_type_cons(eltype(I1), eltype(I2))

iteratoreltype(::Type{Prod{I1,I2}}) where {I1,I2} = and_iteratoreltype(iteratoreltype(I1),iteratoreltype(I2))
iteratorsize(::Type{Prod{I1,I2}}) where {I1,I2} = prod_iteratorsize(iteratorsize(I1),iteratorsize(I2))

@inline function next(p::Prod, st)
    x = prod_next(p, st)
    ((x[1][1],x[1][2]...), x[2])
end

prod_iteratorsize(::Union{HasLength,HasShape}, ::Union{HasLength,HasShape}) = HasShape()
# products can have an infinite iterator
prod_iteratorsize(::IsInfinite, ::IsInfinite) = IsInfinite()
prod_iteratorsize(a, ::IsInfinite) = IsInfinite()
prod_iteratorsize(::IsInfinite, b) = IsInfinite()
prod_iteratorsize(a, b) = SizeUnknown()

# flatten an iterator of iterators

struct Flatten{I}
    it::I
end

"""
    flatten(iter)

Given an iterator that yields iterators, return an iterator that yields the
elements of those iterators.
Put differently, the elements of the argument iterator are concatenated.

# Examples
```jldoctest
julia> collect(Iterators.flatten((1:2, 8:9)))
4-element Array{Int64,1}:
 1
 2
 8
 9
```
"""
flatten(itr) = Flatten(itr)

eltype(::Type{Flatten{I}}) where {I} = eltype(eltype(I))
iteratoreltype(::Type{Flatten{I}}) where {I} = _flatteneltype(I, iteratoreltype(I))
_flatteneltype(I, ::HasEltype) = iteratoreltype(eltype(I))
_flatteneltype(I, et) = EltypeUnknown()

flatten_iteratorsize(::Union{HasShape, HasLength}, b::Type{<:Tuple}) = isleaftype(b) ? HasLength() : SizeUnknown()
flatten_iteratorsize(::Union{HasShape, HasLength}, b::Type{<:Number}) = HasLength()
flatten_iteratorsize(a, b) = SizeUnknown()

iteratorsize(::Type{Flatten{I}}) where {I} = flatten_iteratorsize(iteratorsize(I), eltype(I))

function flatten_length(f, ::Type{T}) where {T<:Tuple}
    if !isleaftype(T)
        throw(ArgumentError(
            "Cannot compute length of a tuple-type which is not a leaf-type"))
    end
    fieldcount(T)*length(f.it)
end
flatten_length(f, ::Type{<:Number}) = length(f.it)
flatten_length(f, T) = throw(ArgumentError(
    "Iterates of the argument to Flatten are not known to have constant length"))
length(f::Flatten{I}) where {I} = flatten_length(f, eltype(I))

function start(f::Flatten)
    local inner, s2
    s = start(f.it)
    d = done(f.it, s)
    # this is a simple way to make this function type stable
    d && throw(ArgumentError("argument to Flatten must contain at least one iterator"))
    while !d
        inner, s = next(f.it, s)
        s2 = start(inner)
        !done(inner, s2) && break
        d = done(f.it, s)
    end
    return s, inner, s2
end

@inline function next(f::Flatten, state)
    s, inner, s2 = state
    val, s2 = next(inner, s2)
    while done(inner, s2) && !done(f.it, s)
        inner, s = next(f.it, s)
        s2 = start(inner)
    end
    return val, (s, inner, s2)
end

@inline function done(f::Flatten, state)
    s, inner, s2 = state
    return done(f.it, s) && done(inner, s2)
end

"""
    partition(collection, n)

Iterate over a collection `n` elements at a time.

# Examples
```jldoctest
julia> collect(Iterators.partition([1,2,3,4,5], 2))
3-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
 [5]
```
"""
partition(c::T, n::Integer) where {T} = PartitionIterator{T}(c, Int(n))

mutable struct PartitionIterator{T}
    c::T
    n::Int
end

eltype(::Type{PartitionIterator{T}}) where {T} = Vector{eltype(T)}
partition_iteratorsize(::HasShape) = HasLength()
partition_iteratorsize(isz) = isz
function iteratorsize(::Type{PartitionIterator{T}}) where {T}
    partition_iteratorsize(iteratorsize(T))
end

function length(itr::PartitionIterator)
    l = length(itr.c)
    return div(l, itr.n) + ((mod(l, itr.n) > 0) ? 1 : 0)
end

start(itr::PartitionIterator) = start(itr.c)

done(itr::PartitionIterator, state) = done(itr.c, state)

function next(itr::PartitionIterator{<:Vector}, state)
    l = state
    r = min(state + itr.n-1, length(itr.c))
    return view(itr.c, l:r), r + 1
end

function next(itr::PartitionIterator, state)
    v = Vector{eltype(itr.c)}(itr.n)
    i = 0
    while !done(itr.c, state) && i < itr.n
        i += 1
        v[i], state = next(itr.c, state)
    end
    return resize!(v, i), state
end

end

max_stars_repo_name: jdlangs/julia

max_issues_repo_name: jdlangs/julia

max_forks_repo_name: jdlangs/julia


text:
#==============================================================================#
# AWSCore/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSCore
using AWSCore: Services
using Dates
using HTTP
using HTTP: Headers, URI
using IniFile
using JSON
using Mocking
using Retry
using SymDict
using Test
using UUIDs
using XMLDict
using .SignatureV4

Mocking.activate()
aws = aws_config()

@testset "AWSCore" begin
    include("aws4.jl")
    include("arn.jl")
    include("credentials.jl")
    include("exceptions.jl")
    include("endpoints.jl")
    include("localhost.jl")
    include("signaturev4.jl")
    include("xml.jl")
end

ext: jl


text:
function __init__cora()
    DEPNAME = "Cora"
    # LINK = "https://github.com/shchur/gnn-benchmark/raw/master/data/npz"
    # LINK = "https://github.com/abojchevski/graph2gauss/raw/master/data/"
    LINK = "https://github.com/kimiyoung/planetoid/raw/master/data"
    DOCS = "https://github.com/kimiyoung/planetoid"
    DATA = "ind.cora." .* ["x", "y", "tx", "allx", "ty", "ally", "graph", "test.index"]

    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $DOCS
        """,
        map(x -> "$LINK/$x", DATA),
        "81de017067dc045ebdb8ffd5c0e69a209973ffdb1fe2d5b434e94d3614f3f5c7",  # if checksum omitted, will be generated by DataDeps
        # post_fetch_method = unpack
    ))
end

"""
    Cora()

The Cora citation network dataset from Ref. [1].
Nodes represent documents and edges represent citation links.
Each node has a predefined feature with 1433 dimensions.
The dataset is designed for the node classification task.
The task is to predict the category of certain paper.
The dataset is retrieved from Ref. [2].

# Statistics

- Nodes: 2708
- Edges: 10556
- Number of Classes: 7
- Label split:
    - Train:  140
    - Val:    500
    - Test:  1000

The split is the one used in the original paper [1] and
doesn't consider all nodes.

# References

[1]: [Deep Gaussian Embedding of Graphs: Unsupervised Inductive Learning via Ranking](https://arxiv.org/abs/1707.03815)
[2]: [Planetoid](https://github.com/kimiyoung/planetoid
"""
struct Cora <: AbstractDataset
    metadata::Dict{String, Any}
    graphs::Vector{Graph}
end

function Cora(; dir=nothing, reverse_edges=true)
    metadata, g = read_planetoid_data("Cora", dir=dir, reverse_edges=reverse_edges)
    return Cora(metadata, [g])
end

Base.length(d::Cora) = length(d.graphs)
Base.getindex(d::Cora) = d.graphs
Base.getindex(d::Cora, i) = getindex(d.graphs, i)


text:
########################
#### Energy balance ####
########################

#' Biochemical Energy
#'
#' Radiant energy absorbed in photosynthesis or heat release by respiration calculated
#'              from net ecosystem exchange of CO2 (NEE).
#'
#' - NEE     Net ecosystem exchange (umol CO2 m-2 s-1)
#' - alpha   Energy taken up/released by photosynthesis/respiration per mol CO2 fixed/respired (J umol-1)
#'
#' # Details
 The following sign convention is employed: NEE is negative when carbon is taken up by the ecosystem.
#'          Positive values of the resulting biochemical energy mean that energy (heat) is taken up by the ecosystem,
#'          negative ones that heat is released.
#'          The value of alpha is taken from Nobel 1974 (see Meyers & Hollinger 2004), but other values
#'          have been used (e.g. Blanken et al., 1997)
#'
#' # Value
 - Sp -: biochemical energy (W m-2)
#'
#' #References
#' Meyers, T_P., Hollinger, S_E. 2004: An assessment of storage terms in the surface energy
#'             balance of maize and soybean. Agricultural and Forest Meteorology 125, 105-115.
#'
#'             Nobel, P_S., 1974: Introduction to Biophysical Plant Physiology.
#'             Freeman, New York.
#'
#'             Blanken, P_D. et al., 1997: Energy balance and canopy conductance of a boreal aspen
#'             forest: Partitioning overstory and understory components.
#'             Journal of Geophysical Research 102, 28915-28927.
#'
#' ```@example; output = false
#' ```
#' # Calculate biochemical energy taken up by the ecosystem with
#' # a measured NEE of -30umol CO2 m-2 s-1
#' biochemical_energy(NEE=-30)
#'
#' @export
function biochemical_energy(NEE,alpha=0.422)
  Sp = alpha*-NEE
  return(Sp)
end

#' Energy-Use Efficiency (EUE)
#'
#' Fraction of net radiation fixed by primary productivity.
#'
#' - GPP     Gross primary productivity exchange (umol CO2 m-2 s-1)
#' - alpha   Energy taken up/released by photosynthesis/respiration (J umol-1)
#' - Rn      Net radiation (W m-2)
#'
#' # Details
 Energy use efficiency is calculated as:
#'
#'            ``EUE = sum(GPP)/sum(Rn)``
#'
#'          where the sums are calculated for complete cases of GPP and Rn over
#'          the entire time period.
#'
#' # Value
 - EUE -: Energy use efficiency (-)
#'
#' #See also
#' [`light_use_efficiency`](@ref)
#'
#' ```@example; output = false
#' ```
#' energy_use_efficiency(GPP=20,Rn=500)
#'
#' @export
function energy_use_efficiency(GPP,alpha=0.422,Rn)

  Sp = biochemical_energy(-GPP,alpha)

  comp  = complete_cases(Sp,Rn)

  Sp_sum = sum(Sp[comp],na_rm=T)
  Rn_sum = sum(Rn[comp],na_rm=T)

  EUE = Sp_sum/Rn_sum

  return(c("EUE"=EUE))
end

#' Energy Balance Closure
#'
#' Calculates the degree of the energy balance non-closure for the entire time span
#'              based on the ratio of two sums (energy balance ratio), and ordinary least squares (OLS).
#'
#' - data  DataFrame or matrix containing all required variables.
#' - Rn    Net radiation (W m-2)
#' - G     Ground heat flux (W m-2); optional
#' - S     Sum of all storage fluxes (W m-2); optional
#' - LE    Latent heat flux (W m-2)
#' - H     Sensible heat flux (W m-2)
#' - instantaneous    should the energy balance be calculated at the time step
#'                         of the observations (`true`), or over the entire time period
#'                         provided as input (`false`)
#' - missing_G_as_NA  if `true`, missing G are treated as `missing`s ,otherwise set to 0.
#' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0.
#'
#'
#' # Details
 The energy balance ratio (EBR) is calculated as:
#'
#'            ``EBR = sum(LE + H)/sum(Rn - G - S)``
#'
#'          the sum is taken for all time steps with complete observations (i.e. where
#'          all energy balance terms are available).
#'
#' # Value
 a named vector containing:
#'         - n: number of complete (all energy balance terms available) observations
#'         - intercept: intercept of the OLS regression
#'         - slope: slope of the OLS regression
#'         - r_squared: r^2 of the OLS regression
#'         - EBR: energy balance ratio
#'
#'         if `instantaneous = true`, only `EBR` is returned.
#'
#' #References
#' Wilson K., et al. 2002: Energy balance closure at FLUXNET sites.
#'             Agricultural and Forest Meteorology 113, 223-243.
#'
#' ```@example; output = false
#' ```
#' ## characterize energy balance closure for DE-Tha in June 2014
#' energy_closure(DE_Tha_Jun_2014,instantaneous=false)
#'
#' ## look at half-hourly closure
#' EBR_inst = energy_closure(DE_Tha_Jun_2014,instantaneous=true)
#' summary(EBR_inst)
#'
#' @importFrom stats complete_cases lm
"""
"""
function energy_closure(data,Rn="Rn",G=nothing,S=nothing,LE="LE",H="H",instantaneous=false,
                           missing_G_as_NA=false,missing_S_as_NA=false)

  check_input(data,list(Rn,LE,H,G,S))

  if(!isnothing(G))
    if (!missing_G_as_NA){G[ismissing(G)] = 0}
else
    cat("Ground heat flux G is not provided and set to 0.",fill=true)
    G = rep(0,nrow(data))
end

  if(!isnothing(S))
    if(!missing_S_as_NA){S[ismissing(S)] = 0 }
else
    cat("Energy storage fluxes S are not provided and set to 0.",fill=true)
    S = rep(0,nrow(data))
end

  if (!instantaneous)
    comp = complete_cases(Rn,G,S,LE,H)
    n    = sum(comp)

    EBR = sum(LE[comp] + H[comp]) / sum(Rn[comp] - G[comp] - S[comp])

    emod = lm(c(LE + H) ~ c(Rn - G - S))
    intercept = summary(emod)$coef[1,1]
    slope     = summary(emod)$coef[2,1]
    r_squared = summary(emod)$r_squared

    return(c("n"=n,"intercept"=round(intercept,3),"slope"=round(slope,3),"r^2"=round(r_squared,3),"EBR"=round(EBR,3)))

else

    EBR = (LE + H) /(Rn - G - S)

    return(EBR)
end
end

#' Isothermal Net Radiation
#'
#' Calculates the isothermal net radiation, i.e. the net radiation
#'              that the surface would receive if it had the same temperature than
#'              the air.
#'
#' - data       DataFrame or matrix containing all required variables
#' - Rn         Net radiation (W m-2)
#' - Tair       Air temperature (degC)
#' - Tsurf      Surface temperature (degC)
#' - emissivity Emissivity of the surface (-)
#' - constants  sigma - Stefan-Boltzmann constant (W m-2 K-4)
#'                   Kelvin - conversion degree Celsius to Kelvin
#'
#' # Details
 The isothermal net radiation (Rni) is given by:
#'
#'            ``Rni = Rn + \epsilon * \sigma * (Tsurf^4 - Tair^4)``
#'
#'          where ``\epsilon`` is the emissivity of the surface. Tsurf and Tair
#'          are in Kelvin.
#'
#' # Value
 - Rni -: isothermal net radiation (W m-2)
#'
#' #References
#' Jones, H. 2014: Plants and Microclimate. 3rd edition, Cambridge
#'             University Press.
#'
#' ```@example; output = false
#' ```
#' # calculate isothermal net radiation of a surface that is 2degC warmer than the air.
#' isothermal_Rn(Rn=400,Tair=25,Tsurf=27,emissivity=0.98)
#'
"""
"""
function isothermal_Rn(data,Rn="Rn",Tair="Tair",Tsurf="Tsurf",emissivity,
                          constants=BigleafConstants())

  check_input(data,list(Rn,Tair,Tsurf))

  Tair  = Tair + constants.Kelvin
  Tsurf = Tsurf + constants.Kelvin

  Rni = Rn + emissivity * constants.sigma * (Tsurf^4 - Tair^4)

  return(Rni)

end


text:
"""
Module Setup2D

Module to aid in setting up reference operators, meshes, geometric terms
Also provides functions to assemble global DG-SBP operators

"""

module SetupDG

# non-DG modules
using LinearAlgebra # for diagm
using SparseArrays # for sparse, droptol
using UnPack # for easy setting/getting in mutable structs

# matlab-like modules
using CommonUtils # for I matrix in geometricFactors
using Basis1D

# triangular routines
import Basis2DTri
import UniformTriMesh # for face vertices

# quadrilateral routines
import Basis2DQuad
import UniformQuadMesh # for face vertices

# hex routines
import Basis3DHex
import UniformHexMesh # for face vertices

# initialization of mesh/reference element data
export init_reference_interval, init_reference_tri
export init_reference_quad, init_reference_hex
export init_mesh
export MeshData, RefElemData

# do we realy need mutable here?  I think we can remove it
# TODO: add types to all of these, siwtch to immutable struct
mutable struct MeshData

    VX;VY;VZ # vertex coordinates
    K::Int # num elems
    EToV # mesh vertex array
    FToF # face connectivity

    x; y; z # physical points
    xf;yf;zf
    xq;yq;zq;wJq # phys quad points, Jacobian-scaled weights

    mapM; mapP; mapB # connectivity between face nodes

    # volume geofacs
    rxJ; sxJ; txJ
    ryJ; syJ; tyJ
    rzJ; szJ; tzJ; J

    # surface geofacs
    nxJ; nyJ; nzJ; sJ

    MeshData() = new() # empty initializer

end

mutable struct RefElemData

    Nfaces; fv # face vertex tuple list

    V1 # low order interp nodes and matrix

    # probably won't use nodes, but might as well keep them around
    r; s; t        # interpolation nodes
    rq; sq; tq; wq # volume quadrature
    rf; sf; tf; wf # surface quadrature
    rp; sp; tp     # plotting nodes

    VDM         # Vandermonde matrix
    Dr; Ds; Dt  # differentiation matrices
    Vq; Vf      # quadrature interpolation matrices
    M; Pq       # mass matrix, L2 projection matrix
    LIFT        # quadrature-based lift matrix
    Vp          # interp to equispaced nodes

    nrJ; nsJ; ntJ # reference normals

    RefElemData() = new() # empty initializer
end

function init_reference_interval(N)
    # initialize a new reference element data struct
    rd = RefElemData()

    # Construct matrices on reference elements
    r,_ = gauss_lobatto_quad(0,0,N)
    VDM = vandermonde_1D(N, r)
    Dr = grad_vandermonde_1D(N, r)/VDM
    @pack! rd = r,VDM,Dr

    V1 = vandermonde_1D(1,r)/vandermonde_1D(1,[-1;1])
    @pack! rd = V1

    rq,wq = gauss_quad(0,0,N+1)
    Vq = vandermonde_1D(N, rq)/VDM
    M = Vq'*diagm(wq)*Vq
    Pq = M\(Vq'*diagm(wq))
    @pack! rd = rq,wq,Vq,M,Pq

    rf = [-1;1]
    nrJ = [-1;1]
    Vf = vandermonde_1D(N,rf)/VDM
    LIFT = M\(transpose(Vf)) # lift matrix
    @pack! rd = rf,nrJ,Vf,LIFT

    # plotting nodes
    rp = LinRange(-1,1,50)
    Vp = vandermonde_1D(N,rp)/VDM
    @pack! rd = rp,Vp

    return rd

end

function init_reference_tri(N)

    # initialize a new reference element data struct
    rd = RefElemData()

    fv = UniformTriMesh.tri_face_vertices() # set faces for triangle
    Nfaces = length(fv)
    @pack! rd = fv, Nfaces

    # Construct matrices on reference elements
    r, s = Basis2DTri.nodes_2D(N)
    VDM = Basis2DTri.vandermonde_2D(N, r, s)
    Vr, Vs = Basis2DTri.grad_vandermonde_2D(N, r, s)
    Dr = Vr/VDM
    Ds = Vs/VDM
    @pack! rd = r,s,VDM,Dr,Ds

    # low order interpolation nodes
    r1,s1 = Basis2DTri.nodes_2D(1)
    V1 = Basis2DTri.vandermonde_2D(1,r,s)/Basis2DTri.vandermonde_2D(1,r1,s1)
    @pack! rd = V1

    #Nodes on faces, and face node coordinate
    r1D, w1D = gauss_quad(0,0,N)
    Nfp = length(r1D) # number of points per face
    e = ones(Nfp) # vector of all ones
    z = zeros(Nfp) # vector of all zeros
    rf = [r1D; -r1D; -e];
    sf = [-e; r1D; -r1D];
    wf = vec(repeat(w1D,3,1));
    nrJ = [z; e; -e]
    nsJ = [-e; e; z]
    @pack! rd = rf,sf,wf,nrJ,nsJ

    rq,sq,wq = Basis2DTri.quad_nodes_2D(2*N)
    Vq = Basis2DTri.vandermonde_2D(N,rq,sq)/VDM
    M = transpose(Vq)*diagm(wq)*Vq
    Pq = M\(transpose(Vq)*diagm(wq))
    @pack! rd = rq,sq,wq,Vq,M,Pq

    Vf = Basis2DTri.vandermonde_2D(N,rf,sf)/VDM # interpolates from nodes to face nodes
    LIFT = M\(transpose(Vf)*diagm(wf)) # lift matrix used in rhs evaluation
    @pack! rd = Vf,LIFT

    # plotting nodes
    rp, sp = Basis2DTri.equi_nodes_2D(15)
    Vp = Basis2DTri.vandermonde_2D(N,rp,sp)/VDM
    @pack! rd = rp,sp,Vp

    return rd
end

# default to full quadrature nodes
# if quad_nodes_1D=tuple of (r1D,w1D) is supplied, use those nodes
function init_reference_quad(N,quad_nodes_1D = gauss_quad(0,0,N))

    # initialize a new reference element data struct
    rd = RefElemData()

    fv = UniformQuadMesh.quad_face_vertices() # set faces for triangle
    Nfaces = length(fv)
    @pack! rd = fv, Nfaces

    # Construct matrices on reference elements
    r, s = Basis2DQuad.nodes_2D(N)
    VDM = Basis2DQuad.vandermonde_2D(N, r, s)
    Vr, Vs = Basis2DQuad.grad_vandermonde_2D(N, r, s)
    Dr = Vr/VDM
    Ds = Vs/VDM
    @pack! rd = r,s,VDM

    # low order interpolation nodes
    r1,s1 = Basis2DQuad.nodes_2D(1)
    V1 = Basis2DQuad.vandermonde_2D(1,r,s)/Basis2DQuad.vandermonde_2D(1,r1,s1)
    @pack! rd = V1

    #Nodes on faces, and face node coordinate
    # r1D,w1D = quad_nodes_1D(0,0,N)
    # r1D,w1D = gauss_lobatto_quad(0,0,N)
    # r1D,w1D = gauss_quad(0,0,N)
    r1D,w1D = quad_nodes_1D
    Nfp = length(r1D)
    e = ones(size(r1D))
    z = zeros(size(r1D))
    rf = [r1D; e; -r1D; -e]
    sf = [-e; r1D; e; -r1D]
    wf = vec(repeat(w1D,Nfaces,1));
    nrJ = [z; e; z; -e]
    nsJ = [-e; z; e; z]
    @pack! rd = rf,sf,wf,nrJ,nsJ

    # quadrature nodes - build from 1D nodes.
    # can also use "rq,sq,wq = Basis2DQuad.quad_nodes_2D(2*N)"
    rq,sq = (x->x[:]).(meshgrid(r1D))
    wr,ws = meshgrid(w1D)
    wq = wr[:] .* ws[:]
    Vq = Basis2DQuad.vandermonde_2D(N,rq,sq)/VDM
    M = transpose(Vq)*diagm(wq)*Vq
    Pq = M\(transpose(Vq)*diagm(wq))
    @pack! rd = rq,sq,wq,Vq,M,Pq

    Vf = Basis2DQuad.vandermonde_2D(N,rf,sf)/VDM # interpolates from nodes to face nodes
    LIFT = M\(transpose(Vf)*diagm(wf)) # lift matrix used in rhs evaluation

    # expose kronecker product sparsity
    Dr = droptol!(sparse(Dr), 1e-10)
    Ds = droptol!(sparse(Ds), 1e-10)
    Vf = droptol!(sparse(Vf),1e-10)
    LIFT = droptol!(sparse(LIFT),1e-10)
    @pack! rd = Dr,Ds,Vf,LIFT

    # plotting nodes
    rp, sp = Basis2DQuad.equi_nodes_2D(15)
    Vp = Basis2DQuad.vandermonde_2D(N,rp,sp)/VDM
    @pack! rd = rp,sp,Vp

    return rd
end

# dispatch to 2D or 3D version if tuple called
function init_mesh(VXYZ,EToV,rd::RefElemData)
    return init_mesh(VXYZ...,EToV,rd)
end

function init_mesh(VX,VY,EToV,rd::RefElemData)

    # initialize a new mesh data struct
    md = MeshData()

    @unpack fv = rd
    FToF = connect_mesh(EToV,fv)
    Nfaces,K = size(FToF)
    @pack! md = FToF,K,VX,VY,EToV

    #Construct global coordinates
    @unpack V1 = rd
    x = V1*VX[transpose(EToV)]
    y = V1*VY[transpose(EToV)]
    @pack! md = x,y

    #Compute connectivity maps: uP = exterior value used in DG numerical fluxes
    @unpack r,s,Vf = rd
    xf,yf = (x->Vf*x).((x,y))
    mapM,mapP,mapB = build_node_maps((xf,yf),FToF)
    Nfp = convert(Int,size(Vf,1)/Nfaces)
    mapM = reshape(mapM,Nfp*Nfaces,K)
    mapP = reshape(mapP,Nfp*Nfaces,K)
    @pack! md = xf,yf,mapM,mapP,mapB

    #Compute geometric factors and surface normals
    @unpack Dr,Ds = rd
    rxJ, sxJ, ryJ, syJ, J = geometric_factors(x, y, Dr, Ds)
    @pack! md = rxJ, sxJ, ryJ, syJ, J

    @unpack Vq,wq = rd
    xq,yq = (x->Vq*x).((x,y))
    wJq = diagm(wq)*(Vq*J)
    @pack! md = xq,yq,wJq

    #physical normals are computed via G*nhatJ, G = matrix of geometric terms
    @unpack nrJ,nsJ = rd
    nxJ = (Vf*rxJ).*nrJ + (Vf*sxJ).*nsJ;
    nyJ = (Vf*ryJ).*nrJ + (Vf*syJ).*nsJ;
    sJ = @. sqrt(nxJ^2 + nyJ^2)
    @pack! md = nxJ,nyJ,sJ

    return md
end

"========== 3D routines ============="

function init_reference_hex(N,quad_nodes_1D=gauss_quad(0,0,N))

    # initialize a new reference element data struct
    rd = RefElemData()

    fv = UniformHexMesh.hex_face_vertices() # set faces for triangle
    Nfaces = length(fv)
    @pack! rd = fv, Nfaces

    # Construct matrices on reference elements
    r,s,t = Basis3DHex.nodes_3D(N)
    VDM = Basis3DHex.vandermonde_3D(N,r,s,t)
    Vr,Vs,Vt = Basis3DHex.grad_vandermonde_3D(N,r,s,t)
    Dr,Ds,Dt = (A->A/VDM).(Basis3DHex.grad_vandermonde_3D(N,r,s,t))
    @pack! rd = r,s,t,VDM

    # low order interpolation nodes
    r1,s1,t1 = Basis3DHex.nodes_3D(1)
    V1 = Basis3DHex.vandermonde_3D(1,r,s,t)/Basis3DHex.vandermonde_3D(1,r1,s1,t1)
    @pack! rd = V1

    #Nodes on faces, and face node coordinate
    r1D,w1D = quad_nodes_1D
    rquad,squad = (x->x[:]).(meshgrid(r1D,r1D))
    wr,ws = (x->x[:]).(meshgrid(w1D,w1D))
    wquad = wr.*ws
    e = ones(size(rquad))
    zz = zeros(size(rquad))
    rf = [-e; e; rquad; rquad; rquad; rquad]
    sf = [rquad; rquad; -e; e; squad; squad]
    tf = [squad; squad; squad; squad; -e; e]
    wf = vec(repeat(wquad,Nfaces,1));
    nrJ = [-e; e; zz;zz; zz;zz]
    nsJ = [zz;zz; -e; e; zz;zz]
    ntJ = [zz;zz; zz;zz; -e; e]

    @pack! rd = rf,sf,tf,wf,nrJ,nsJ,ntJ

    # quadrature nodes - build from 1D nodes.
    rq,sq,tq = meshgrid(r1D,r1D,r1D)
    wr,ws,wt = meshgrid(w1D,w1D,w1D)
    wq = wr.*ws.*wt
    Vq = Basis3DHex.vandermonde_3D(N,rq,sq,tq)/VDM
    M = transpose(Vq)*diagm(wq)*Vq
    Pq = M\(transpose(Vq)*diagm(wq))
    @pack! rd = rq,sq,tq,wq,Vq,M,Pq

    Vf = Basis3DHex.vandermonde_3D(N,rf,sf,tf)/VDM
    LIFT = M\(transpose(Vf)*diagm(wf))

    # expose kronecker product sparsity
    Dr = droptol!(sparse(Dr),1e-12)
    Ds = droptol!(sparse(Ds),1e-12)
    Dt = droptol!(sparse(Dt),1e-12)
    Vf = droptol!(sparse(Vf),1e-12)
    LIFT = droptol!(sparse(LIFT),1e-12)
    @pack! rd = Dr,Ds,Dt,Vf,LIFT

    # plotting nodes
    rp,sp,tp = Basis3DHex.equi_nodes_3D(15)
    Vp = Basis3DHex.vandermonde_3D(N,rp,sp,tp)/VDM
    @pack! rd = rp,sp,tp,Vp

    return rd
end

function init_mesh(VX,VY,VZ,EToV,rd::RefElemData)

    # initialize a new mesh data struct
    md = MeshData()

    @unpack fv = rd
    FToF = connect_mesh(EToV,fv)
    Nfaces,K = size(FToF)
    @pack! md = FToF,K,VX,VY,VZ,EToV

    #Construct global coordinates
    @unpack V1 = rd
    x = V1*VX[transpose(EToV)]
    y = V1*VY[transpose(EToV)]
    z = V1*VZ[transpose(EToV)]
    @pack! md = x,y,z

    #Compute connectivity maps: uP = exterior value used in DG numerical fluxes
    @unpack r,s,t,Vf = rd
    xf,yf,zf = (x->Vf*x).((x,y,z))
    mapM,mapP,mapB = build_node_maps((xf,yf,zf),FToF)
    Nfp = convert(Int,size(Vf,1)/Nfaces)
    mapM = reshape(mapM,Nfp*Nfaces,K)
    mapP = reshape(mapP,Nfp*Nfaces,K)
    @pack! md = xf,yf,zf,mapM,mapP,mapB

    #Compute geometric factors and surface normals
    @unpack Dr,Ds,Dt = rd
    rxJ,sxJ,txJ,ryJ,syJ,tyJ,rzJ,szJ,tzJ,J = geometric_factors(x,y,z,Dr,Ds,Dt)
    @pack! md = rxJ,sxJ,txJ,ryJ,syJ,tyJ,rzJ,szJ,tzJ,J

    @unpack Vq,wq = rd
    xq,yq,zq = (x->Vq*x).((x,y,z))
    wJq = diagm(wq)*(Vq*J)
    @pack! md = xq,yq,zq,wJq

    #physical normals are computed via G*nhatJ, G = matrix of geometric terms
    @unpack nrJ,nsJ,ntJ = rd
    nxJ = nrJ.*(Vf*rxJ) + nsJ.*(Vf*sxJ) + ntJ.*(Vf*txJ)
    nyJ = nrJ.*(Vf*ryJ) + nsJ.*(Vf*syJ) + ntJ.*(Vf*tyJ)
    nzJ = nrJ.*(Vf*rzJ) + nsJ.*(Vf*szJ) + ntJ.*(Vf*tzJ)
    sJ = @. sqrt(nxJ.^2 + nyJ.^2 + nzJ.^2)
    @pack! md = nxJ,nyJ,nzJ,sJ

    return md
end

end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using Revise
using JDF
using CSV, DataFrames, Feather

@time a = CSV.read("c:/data/a.feather");

@time pmetadatas = psavejdf(a, "c:/data/a.jdf");

@time metadatas = savejdf(a, "c:/data/a.jdf");
#a = nothing
@time a2 = loadjdf("c:/data/a.jdf", metadatas);

all(names(a) .== names(a2))
all(skipmissing([all(a2[!,name] .== Array(a[!,name])) for name in names(a2)]))


text:
module SubModule1
    import AbstractModule.foo
    foo(x::Int) = println("Hello there Int")
end


text:
function sinpi(x::Real)
    if isinf(x)
        return throw(DomainError())
    elseif isnan(x)
        return nan(x)
    end

    rx = copysign(float(rem(x,2)),x)
    arx = abs(rx)

    if arx < oftype(rx,0.25)
        return sin(pi*rx)
    elseif arx <= oftype(rx,0.75)
        arx = oftype(rx,0.5) - arx
        return copysign(cos(pi*arx),rx)
    elseif arx < oftype(rx,1.25)
        rx = (one(rx) - arx)*sign(rx)
        return sin(pi*rx)
    elseif arx <= oftype(rx,1.75)
        arx = oftype(rx,1.5) - arx
        return -copysign(cos(pi*arx),rx)
    else
        rx = rx - copysign(oftype(rx,2.0),rx)
        return sin(pi*rx)
    end
end

function cospi(x::Real)
    if isinf(x)
        return throw(DomainError())
    elseif isnan(x)
        return nan(x)
    end

    rx = abs(float(rem(x,2)))

    if rx <= oftype(rx,0.25)
        return cos(pi*rx)
    elseif rx < oftype(rx,0.75)
        rx = oftype(rx,0.5) - rx
        return sin(pi*rx)
    elseif rx <= oftype(rx,1.25)
        rx = one(rx) - rx
        return -cos(pi*rx)
    elseif rx < oftype(rx,1.75)
        rx = rx - oftype(rx,1.5)
        return sin(pi*rx)
    else
        rx = oftype(rx,2.0) - rx
        return cos(pi*rx)
    end
end

sinpi(x::Integer) = zero(x)
cospi(x::Integer) = isodd(x) ? -one(x) : one(x)

function sinpi(z::Complex)
    zr, zi = reim(z)
    if !isfinite(zi) && zr == 0 return complex(zr, zi) end
    if isnan(zr) && !isfinite(zi) return complex(zr, zi) end
    if !isfinite(zr) && zi == 0 return complex(oftype(zr, NaN), zi) end
    if !isfinite(zr) && isfinite(zi) return complex(oftype(zr, NaN), oftype(zi, NaN)) end
    if !isfinite(zr) && !isfinite(zi) return complex(zr, oftype(zi, NaN)) end
    pizi = pi*zi
    complex(sinpi(zr)*cosh(pizi), cospi(zr)*sinh(pizi))
end

function cospi(z::Complex)
    zr, zi = reim(z)
    if !isfinite(zi) && zr == 0
        return complex(isnan(zi) ? zi : oftype(zi, Inf),
                       isnan(zi) ? zr : zr*-sign(zi))
    end
    if !isfinite(zr) && isinf(zi)
        return complex(oftype(zr, Inf), oftype(zi, NaN))
    end
    if isinf(zr)
        return complex(oftype(zr, NaN), zi==0 ? -copysign(zi, zr) : oftype(zi, NaN))
    end
    if isnan(zr) && zi==0 return complex(zr, abs(zi)) end
    pizi = pi*zi
    complex(cospi(zr)*cosh(pizi), -sinpi(zr)*sinh(pizi))
end
@vectorize_1arg Number sinpi
@vectorize_1arg Number cospi

sinc(x::Number) = x==0 ? one(x)  : oftype(x,sinpi(x)/(pi*x))
sinc(x::Integer) = x==0 ? one(x) : zero(x)
sinc{T<:Integer}(x::Complex{T}) = sinc(float(x))
@vectorize_1arg Number sinc
cosc(x::Number) = x==0 ? zero(x) : oftype(x,(cospi(x)-sinpi(x)/(pi*x))/x)
cosc(x::Integer) = cosc(float(x))
cosc{T<:Integer}(x::Complex{T}) = cosc(float(x))
@vectorize_1arg Number cosc

for (finv, f) in ((:sec, :cos), (:csc, :sin), (:cot, :tan),
                  (:sech, :cosh), (:csch, :sinh), (:coth, :tanh),
                  (:secd, :cosd), (:cscd, :sind), (:cotd, :tand))
    @eval begin
        ($finv){T<:Number}(z::T) = one(T) / (($f)(z))
        ($finv){T<:Number}(z::AbstractArray{T}) = one(T) ./ (($f)(z))
    end
end

for (fa, fainv) in ((:asec, :acos), (:acsc, :asin), (:acot, :atan),
                    (:asech, :acosh), (:acsch, :asinh), (:acoth, :atanh))
    @eval begin
        ($fa){T<:Number}(y::T) = ($fainv)(one(T) / y)
        ($fa){T<:Number}(y::AbstractArray{T}) = ($fainv)(one(T) ./ y)
    end
end

function sind(x::Real)
    if isinf(x)
        return throw(DomainError())
    elseif isnan(x)
        return nan(x)
    end

    rx = copysign(float(rem(x,360)),x)
    arx = abs(rx)

    if arx < oftype(rx,45.0)
        return sin(deg2rad(rx))
    elseif arx <= oftype(rx,135.0)
        arx = oftype(rx,90.0) - arx
        return copysign(cos(deg2rad(arx)),rx)
    elseif arx < oftype(rx,225.0)
        rx = (oftype(rx,180.0) - arx)*sign(rx)
        return sin(deg2rad(rx))
    elseif arx <= 315.0
        arx = oftype(rx,270.0) - arx
        return -copysign(cos(deg2rad(arx)),rx)
    else
        rx = rx - copysign(oftype(rx,360.0),rx)
        return sin(deg2rad(rx))
    end
end
@vectorize_1arg Real sind

function cosd(x::Real)
    if isinf(x)
        return throw(DomainError())
    elseif isnan(x)
        return nan(x)
    end

    rx = abs(float(rem(x,360)))

    if rx <= oftype(rx,45.0)
        return cos(deg2rad(rx))
    elseif rx < oftype(rx,135.0)
        rx = oftype(rx,90.0) - rx
        return sin(deg2rad(rx))
    elseif rx <= oftype(rx,225.0)
        rx = oftype(rx,180.0) - rx
        return -cos(deg2rad(rx))
    elseif rx < oftype(rx,315.0)
        rx = rx - oftype(rx,270.0)
        return sin(deg2rad(rx))
    else
        rx = oftype(rx,360.0) - rx
        return cos(deg2rad(rx))
    end
end
@vectorize_1arg Real cosd

tand(x::Real) = sind(x) / cosd(x)
@vectorize_1arg Real tand

for (fd, f) in ((:sind, :sin), (:cosd, :cos), (:tand, :tan))
    @eval begin
        ($fd)(z) = ($f)(deg2rad(z))
    end
end

for (fd, f) in ((:asind, :asin), (:acosd, :acos), (:atand, :atan),
                (:asecd, :asec), (:acscd, :acsc), (:acotd, :acot))
    @eval begin
        ($fd)(y) = rad2deg(($f)(y))
        @vectorize_1arg Real $fd
    end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
include("probabilities.jl")

type BloomFilter
    array::BitArray
    k::Int
    capacity::Int
    error_rate::Float64
    n_bits::Int
    mmap_location::String
end

### Hash functions (uses 2 hash method)
# Uses MurmurHash on 64-bit systems so sufficient randomness/speed
# Get the nth hash of a string using the formula hash_a + n * hash_b
# which uses 2 hash functions vs. k and has comparable properties
# See Kirsch and Mitzenmacher, 2008: http://www.eecs.harvard.edu/~kirsch/pubs/bbbf/rsa.pdf
function hash_n(key::String, k::Int, max::Int)
    a_hash = hash(key, uint(0))
    b_hash = hash(key, uint(170))
    hashes = Array(Uint, k)
    for i in 1:k
        hashes[i] = mod(a_hash + i * b_hash, max) + 1
    end
    return hashes
end

### Bloom filter constructors
# Constructor group #1: (optional IOStream / string to mmap-array), capacity, bits per element, k
# Note that this inserts NaN for the error rate, as this module does not include functionality
# for calculating error rates for any arbitrary bits per element and k (only a subset, see #2 below)
# USE CASES: Advanced use only. Not recommended as default constructor choice.
function BloomFilter(capacity::Int, bits_per_elem::Int, k_hashes::Int)
    n_bits = capacity * bits_per_elem
    BloomFilter(falses((1, n_bits)), k_hashes, capacity, NaN, n_bits, "")
end

function BloomFilter(mmap_handle::IOStream, capacity::Int, bits_per_elem::Int, k_hashes::Int)
    n_bits = capacity * bits_per_elem
    mb = mmap_bitarray((n_bits, 1), mmap_handle)
    BloomFilter(mb, k_hashes, capacity, NaN, n_bits, mmap_handle.name)
end

function BloomFilter(mmap_string::String, capacity::Int, bits_per_elem::Int, k_hashes::Int)
    if isfile(mmap_string)
        mmap_handle = open(mmap_string, "r+")
    else
        mmap_handle = open(mmap_string, "w+")
    end
    BloomFilter(mmap_handle, capacity, bits_per_elem, k_hashes)
end

# Constructor group #2: (optional IOStream / string to mmap-array), capacity, error rate, k hashes
# Looks up the optimal number of bits per element given an error rate and specified k in a pre-calculated
# probability table. Note that k must be <= 12 or an error will be thrown. Similarly, this method
# does not support Bloom filter construction where more than 4 bytes are required per element
# (though that can manually be accomplished with one of the above constrcutors).
# USE CASES: Recommended for most applications, but carries modest space-tradeoff for slightly faster operation with fewer hashes
function BloomFilter(capacity::Int, error_rate::Float64, k_hashes::Int)
    bits_per_elem, error_rate = get_k_error(error_rate, k_hashes)
    n_bits = capacity * bits_per_elem
    BloomFilter(falses((n_bits, 1)), k_hashes, capacity, error_rate, n_bits, "")
end

function BloomFilter(mmap_handle::IOStream, capacity::Int, error_rate::Float64, k_hashes::Int)
    bits_per_elem, error_rate = get_k_error(error_rate, k_hashes)
    n_bits = capacity * bits_per_elem
    mb = mmap_bitarray((n_bits, 1), mmap_handle)
    BloomFilter(mb, k_hashes, capacity, error_rate, n_bits, mmap_handle.name)
end

function BloomFilter(mmap_string::String, capacity::Int, error_rate::Float64, k_hashes::Int)
    if isfile(mmap_string)
        mmap_handle = open(mmap_string, "r+")
    else
        mmap_handle = open(mmap_string, "w+")
    end
    BloomFilter(mmap_handle, capacity, error_rate, k_hashes)
end

# Constructor group #3: (optional IOStream / string to mmap-array), capacity, error rate
# Computes an optimal k for a given error rate, where k is selected to minimize the overall
# space required for the Bloom filter. In practice, k may be larger than desired and require
# additional memory accesses.
# USE CASES: Recommended when extreme space efficiency is required, and modestly slower insertions
# and lookups are tolerable.
function BloomFilter(capacity::Int, error_rate::Float64)
    bits_per_elem = int(ceil(-1.0 * (log(error_rate) / (log(2) ^ 2))))
    k_hashes = int(round(log(2) * bits_per_elem))  # Note: ceil() would be strictly more conservative
    n_bits = capacity * bits_per_elem
    BloomFilter(falses((1, n_bits)), k_hashes, capacity, error_rate, n_bits, "")
end

function BloomFilter(mmap_handle::IOStream, capacity::Int, error_rate::Float64)
    bits_per_elem = int(ceil(-1.0 * (log(error_rate) / (log(2) ^ 2))))
    k_hashes = int(round(log(2) * bits_per_elem))  # Note: ceil() would be strictly more conservative
    n_bits = capacity * bits_per_elem
    mb = mmap_bitarray((n_bits, 1), mmap_handle)
    BloomFilter(mb, k_hashes, capacity, error_rate, n_bits, mmap_handle.name)
end

function BloomFilter(mmap_string::String, capacity::Int, error_rate::Float64)
    if isfile(mmap_string)
        mmap_handle = open(mmap_string, "r+")
    else
        mmap_handle = open(mmap_string, "w+")
    end
    BloomFilter(mmap_handle, capacity, error_rate)
end

### Bloom filter functions: insert!, add! (alias to insert), contains, and show
function add!(bf::BloomFilter, key::String)
    hashes = hash_n(key, bf.k, bf.n_bits)
    for h in hashes
        bf.array[h] = 1
    end
end

function contains(bf::BloomFilter, key::String)
    hashes = hash_n(key, bf.k, bf.n_bits)
    for h in hashes
        if bf.array[h] != 1
            return false
        end
    end
    return true
end

function show(io::IO, bf::BloomFilter)
    @printf "Bloom filter with capacity %d, " bf.capacity
    @printf "error rate of %.2f%%, and k of %d.\n" (bf.error_rate * 100) (bf.k)
    @printf "Total bits required: %d (%.1f / element).\n" bf.n_bits (bf.n_bits / bf.capacity)
    if bf.mmap_location != ""
        @printf "Bloom filter is backed by mmap at %s." bf.mmap_location
    else
        @printf "Bloom filter is in-memory."
    end
end


text:
####
#### Unregularized NNLS problem
####

struct NNLSProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W
end
function NNLSProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLS.NNLSWorkspace(C, d)
    NNLSProblem(C, d, m, n, nnls_work)
end

function solve!(work::NNLSProblem, C, d)
    # Solve NNLS problem
    NNLS.load!(work.nnls_work, C, d)
    NNLS.nnls!(work.nnls_work)
    return solution(work)
end
solve!(work::NNLSProblem) = solve!(work, work.C, work.d)

solution(work::NNLSProblem) = NNLS.solution(work.nnls_work)

chi2(work::NNLSProblem) = NNLS.residualnorm(work.nnls_work)^2

"""
    lsqnonneg(C::AbstractMatrix, d::AbstractVector)

Returns the nonnegative least-squares (NNLS) solution, X, of the equation:

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2
```

# Arguments
- `C::AbstractMatrix`: Left hand side matrix acting on `x`
- `d::AbstractVector`: Right hand side vector

# Outputs
- `X::AbstractVector`: NNLS solution
"""
lsqnonneg(C, d) = lsqnonneg!(lsqnonneg_work(C, d))
lsqnonneg_work(C, d) = NNLSProblem(C, d)
lsqnonneg!(work::NNLSProblem) = solve!(work)
lsqnonneg!(work::NNLSProblem{T}, C::AbstractMatrix{T}, d::AbstractVector{T}) where {T} = solve!(work, C, d)


text:
abstract type AbstractParameters end

#holds the main parameter object
mutable struct Parameters <: AbstractParameters
  σ²G::Float64
  ζ²G::Float64
  ζ²P::Float64
  SG::Float64
  SGP::Float64
end

#named constructor for testing purposes
function Parameters(;
  σ²G::Float64 = D_σ²G,
  ζ²G::Float64 = D_ζ²G,
  ζ²P::Float64 = D_ζ²P,
  SG::Float64 = D_SG,
  SGP::Float64 = D_SGP)::Parameters

  return Parameters(σ²G, ζ²G, ζ²P, SG, SGP)
end

function Base.print(Θ::Parameters)
  return  "σ²G = $(Θ.σ²G); ζ²G = $(Θ.ζ²G); ζ²P = $(Θ.ζ²P); SG = $(Θ.SG); SGP = $(Θ.SGP)"
end

Base.println(Θ::Parameters)::String = Base.print(Θ) * "\n"

copy2vector(Θ::Parameters)::Vector{Float64} = [Θ.σ²G, Θ.ζ²G, Θ.ζ²P, Θ.SG, Θ.SGP]

Base.names(Θ::Parameters) = ["sigma2G", "zeta2G", "zeta2P", "SG", "SGP"]

#holds the priors
struct Priors <: AbstractParameters
  θG::Float64
  δ²G::Float64
  αG::Float64
  βG::Float64
  αP::Float64
  βP::Float64
  θSG::Float64
  δ²SG::Float64
  θSGP::Float64
  δ²SGP::Float64
end

#named constructor for testing purposes
function Priors(;  θG::Float64 = D_θG,
  δ²G::Float64 = D_δ²G,
  αG::Float64 = D_αG,
  βG::Float64 = D_βG,
  αP::Float64 = D_αP,
  βP::Float64 = D_βP,
  θSG::Float64 = D_θSG,
  δ²SG::Float64 = D_δ²SG,
  θSGP::Float64 = D_θSGP,
  δ²SGP::Float64 = D_δ²SGP)::Priors

  return Priors(θG,  δ²G,  αG,  βG,  αP,  βP,  θSG,  δ²SG,  θSGP,  δ²SGP)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
    Pratt11 hopper policy for Hopper2D
"""

@with_kw mutable struct Pratt11Options
    live_plotting::Bool=false # Use the live plotting tool to debug
end

mutable struct Pratt11 <: Policy
	model::Flamingo
	phase::Symbol
	count::Int
	h_sim::Real
	u::AbstractVector
	w::AbstractVector
	contact::Vector{Bool}
	front::Symbol
	rear::Symbol
	q0::AbstractVector
	q1::AbstractVector
	qref::AbstractVector
	xdref::Real
	opts::Pratt11Options
end

function pratt_policy(model::Flamingo, h_sim;
		qref=[0.0, 0.849, -0.00, 0.1, 0.295, -0.3, 0.1, π/2, π/2],
		# qref = [0.0, 0.849, -0.00, -0.3, 0.1, 0.1, 0.295, π/2, π/2],
		# qref = [
		# 		0.0,
		# 		0.7957905840060919,
		# 		-0.0,
		# 		-0.5609986881410345,
		# 		0.20943951023931953,
		# 		-0.4487989505128276,
		# 		0.3830761319035729,
		# 		1.5707963267948966,
		# 		1.5707963267948966,
		# ],
		xdref=-0.10,
		mpc_opts = Pratt11Options(),
		)

	u  = zeros(model.dim.u)
	w  = zeros(model.dim.u+1)
	q0 = zeros(model.dim.q)
	q1 = zeros(model.dim.q)
	contact = [false for i=1:model.dim.c]
	Pratt11(model, :settle, 0, h_sim, u, w, contact, :foot_1, :foot_2,
		copy(q0), copy(q1), qref, xdref, mpc_opts)
end

function policy(p::Pratt11, x, traj, t)
	nc = p.model.dim.c
	iθ = Vector(3:9)
	il1 = [2,3,6]
	il2 = [4,5,7]
	m_flamingo = model.m_torso + model.m_thigh1 +
		model.m_thigh2 + model.m_calf1 + model.m_calf2
	# Initialization
	p.q0 = copy(traj.q[t])
	p.q1 = copy(x)
	qd = (p.q1 - p.q0)/traj.h

	# Detect contact
	p.contact = detect_contact(model, traj.γ[max(1,t-1)])
	# Compute feet positions
	xr = kinematics_3(model, p.q1, body=p.rear, mode=:com)[1]
	xf = kinematics_3(model, p.q1, body=p.front, mode=:com)[1]

	if p.front == :foot_1
		ilf = il1
		ilr = il2
	elseif p.front == :foot_2
		ilf = il2
		ilr = il1
	end

	# Update phase
	p.count += 1
	if p.phase == :settle && all(p.contact) && p.count*p.h_sim >= 0.278
		p.front, p.rear = foot_order(model, p.q1)
		p.phase = :translation
		p.count = 0
	# elseif p.phase == :translation && ((p.q1[1] - xr > 0.26) || (xf - p.q1[1] < 0.15)) && p.count*p.h_sim >= 10.926
	# 	p.front, p.rear = foot_order(model, p.q1)
	# 	p.phase = :swing
	# 	p.count = 0
	# elseif p.phase == :swing && ((xr - p.q1[1] > 0.20) || (p.q1[1] - xf > 0.10) || (p.q1[ilr[1]] > +0.0)) && p.count*p.h_sim >= 0.926
	# 	p.front, p.rear = foot_order(model, p.q1)
	# 	@show :translation
	# 	p.phase = :translation
	# 	p.count = 0
	end

	# @show t
	# @show p.count
	@show p.phase
	# @show p.front
	# @show p.contact

	if p.front == :foot_1
		ilf = il1
		ilr = il2
	elseif p.front == :foot_2
		ilf = il2
		ilr = il1
	end

	# # PD joint controller
	# kp = -50.0
	# kd =  0.04*kp
	# p.w[1] = kp * (p.q1[3] - p.qref[3]) + kd * qd[3]

	ang0 = config_to_angles(p.q0)
	ang1 = config_to_angles(p.q1)
	angref = config_to_angles(p.qref)
	angd = (ang1 - ang0) / traj.h
	if p.phase == :settle
		# # PD joint controller
		# kα = 100.0
		# kβ = 10.0
		# kp = -[kα, kα, kα, kα, kα, kβ, kβ]
		# kd =  0.04*kp
		# p.w[1:7] = kp .* (p.q1[iθ] - p.qref[iθ]) + kd .* (qd[iθ])

		# PD joint controller
		kα = 100.0
		kβ = 30.0
		kp = -[kα, kα, kα, kα, kβ, kβ]
		kd = 0.04*kp
		p.w[1:6] = kp .* (ang1 - angref) + kd .* angd
	end

	if p.phase == :translation
		kpfx = -200.0
		kdfx = 0.04*kpfx
		kpfz = -400.0
		kdfz = -200.0

		xref = p.q1[1] + p.xdref*p.h_sim
		fx = kpfx*(p.q1[1] - xref) + kdfx*(qd[1] - p.xdref)
		fz = kpfz*(p.q1[2] - p.qref[2]+0.02) + kdfz*qd[2] + model.g*m_flamingo
		f = [fx, fz]
		# α = 0.5
		α = 0.25 + 0.5*(p.q1[1] -0.10 - xr)/(xf-xr)
		α = clamp(α, 0.25, 0.75)
		# @show α
		ff = α*f
		fr = (1-α)*f
		τf = virtual_actuator_torque(p.model, p.q1, ff, body=p.front)
		τr = virtual_actuator_torque(p.model, p.q1, fr, body=p.rear)
		p.w[ilf] .= τf
		p.w[ilr] .= τr
		p.w[ilr[2]] += -0.5*qd[2 + ilr[2]]
	end

	# if p.phase == :swing
	# 	kpfx = -200.0
	# 	kdfx = 0.04*kpfx
	# 	kpfz = -400.0
	# 	kdfz = -200.0
	#
	# 	xref = p.q1[1] + p.xdref*p.h_sim
	# 	fx = kpfx*(p.q1[1] - xref) + kdfx*(qd[1] - p.xdref)
	# 	fz = kpfz*(p.q1[2] - p.qref[2]+0.02) + kdfz*qd[2] + model.g*m_flamingo
	# 	f = [fx, fz]
	# 	α = 1.0
	# 	ff = α*f
	# 	τf = virtual_actuator_torque(p.model, p.q1, ff, body=p.front)
	# 	p.w[ilf] .= τf
	#
	# 	p.w[ilf[1]] += -2*(p.q1[2 + ilf[1]] + 0.2)     - 0.2*qd[2 + ilf[1]]
	# 	p.w[ilf[2]] += -2*(p.q1[2 + ilf[2]] + 0.2)     - 0.2*qd[2 + ilf[2]]
	#
	# 	p.w[ilr[1]] = -3*(p.q1[2 + ilr[1]] + 1.0)     - 1.5*(qd[2 + ilr[1]] - 0.6)
	# 	p.w[ilr[2]] = -3*(p.q1[2 + ilr[2]] - 0.6)     - 0.2*qd[2 + ilr[2]]
	# 	p.w[ilr[3]] = -0.3*(p.q1[end] - π/2-0.3) - 0.2*qd[end]
	# end

	# # Ankle spring damper
	# p.w[ilr[3]] = -5*(p.q1[2 + ilr[3]] - π/2-p.q1[2 + ilr[2]]) - 0.05*qd[2 + ilr[3]]
	# p.w[ilf[3]] = -5*(p.q1[2 + ilf[3]] - π/2-p.q1[2 + ilf[2]]) - 0.05*qd[2 + ilf[3]]

	# B = B_func(model, p.q1)[:,3:end]
	# # @show size(p.w)
	# # @show size(B*B')
	# # @show size(B)
	# p.u = (B*B')\(B*p.w)
	# @show norm(p.w - B'*p.u)/norm(p.w)
	# @show scn.(p.w)
	# @show scn.(B'*p.u)
	p.u = p.w[1:6]
	# Rescale
	p.u .*= traj.h
    return p.u
end

function config_to_angles(q::AbstractVector{T}) where T
	ang = zeros(T,6)
	ang[1] = q[4] - q[3]
	ang[2] = q[5] - q[4]
	ang[3] = q[6] - q[3]
	ang[4] = q[7] - q[6]
	ang[5] = q[8] - q[5]
	ang[6] = q[9] - q[7]
	return ang
end

function detect_contact(model::Flamingo, γ::AbstractVector; contact_threshold=1.5e-2)
	nc = model.dim.c
	contact = [false for i=1:nc]
	for i = 1:nc
		if γ[i] .> contact_threshold
			contact[i] = true
		else
			contact[i] = false
		end
	end
	return contact
end

function foot_order(model::Flamingo,  q::AbstractVector)
	c1 = kinematics_3(model, q, body=:foot_1, mode=:com)
	c2 = kinematics_3(model, q, body=:foot_2, mode=:com)
	if c1[1] > c2[1]
		return :foot_1, :foot_2
	else
		return :foot_2, :foot_1
	end
	return nothing
end

function kinematic_map(model::Flamingo, q; body=:foot_1)
	#TODO need to compute wrt to the center of pressure not com
	x, z = q[1:2] - kinematics_3(model, q, body=body, mode=:com)
	return [x,z]
end
function jacobian_map(model::Flamingo, q; body=:foot_1)
	k(q) = kinematic_map(model, q, body=body)
	J = ForwardDiff.jacobian(k, q)[:,3:end]
	return J
end

function virtual_actuator_torque(model::Flamingo, q::AbstractVector,
		f::AbstractVector; body=:foot_1)
	J = jacobian_map(model, q, body=body)
	if body == :foot_1
		iJ = [2,3,6]
	elseif body == :foot_2
		iJ = [4,5,7]
	else
		@error "incorrect body specification"
	end
	τ = J[:,iJ]'*f
	return τ
end


text:
function renormalize{T}(tensorT1::Array{T,6}, tensorT2::Array{T,6}, dimM::Int; whichIsTensorT::Int = 1)
	if whichIsTensorT ==1
		tensorT = tensorT1
		tensorTtilde = tensorT2
	else
		tensorT = tensorT2
		tensorTtilde = tensorT1
	end
	tensorU = getTensorU(tensorT, dimM)
	tensorV = getTensorV(tensorT, dimM)
	newTensorT = getNewTensorT_2dQ(tensorT, tensorT, tensorU, tensorV)
	normalizationFactor = maximum(newTensorT)
	newTensorTtilde = getNewTensorT_2dQ(tensorT1, tensorT2, tensorU, tensorV)
	newTensorTtilde = newTensorTtilde./ normalizationFactor
	newTensorT = newTensorT./ normalizationFactor
	return newTensorT, newTensorTtilde, normalizationFactor
end

function getNewTensorT_2dQ{T}(tensorT1::Array{T,6}, tensorT2::Array{T,6}, tensorU::Array{T,3}, tensorV::Array{T,3})
	@tensor begin
		d1[x2,x,x1p,y1,y1p,z,i] := tensorT1[x1,x1p,y1,y1p,z,i] * tensorU[x1,x2,x]
		d2[x2,x,x2p,xp,y1,y1p,z,i] := d1[x2,x,x1p,y1,y1p,z,i] * tensorU[x1p,x2p,xp]
		d3[x,xp,y1,y1p,y2,y2p,z,zp] := d2[x2,x,x2p,xp,y1,y1p,z,i] * tensorT2[x2,x2p,y2,y2p,i,zp]
		d4[x,xp,y,y1p,y2p,z,zp] := d3[x,xp,y1,y1p,y2,y2p,z,zp] * tensorV[y1,y2,y]
		newTensorT[x,xp,y,yp,z,zp] := d4[x,xp,y,y1p,y2p,z,zp] * tensorV[y1p,y2p,yp]
	end
	return newTensorT
end

function getTensorU{T}(tensorT::Array{T,6}, dimM::Int)
	tenmatMMd = getTenmatMMd(tensorT)
	Ul, lambdaVector, Uld = svd(tenmatMMd.matrix)
	trunUl = truncMatrixU(Ul,dimM)
	tensorU = matU2tenU(trunUl)
	return tensorU
end

function getTensorV{T}(tensorT::Array{T,6}, dimM::Int)
	tensorV = getTensorU(permutedims(tensorT,[3,4,1,2,5,6]), dimM)
	return tensorV
end

function getTenmatMMd{T}(tensorT::Array{T,6})
	@tensor begin
		dummy1[x1,a1,i,j] := tensorT[x1,x1p,y1,y1p,z,i] *
			tensorT[a1,x1p,y1,y1p,z,j]
		dummy2[x2,a2,i,j] := tensorT[x2,x2p,y2,y2p,i,zp] *
			tensorT[a2,x2p,y2,y2p,j,zp]
		tensorMMd[x1,x2,a1,a2] := dummy1[x1,a1,i,j] *
			dummy2[x2,a2,i,j]
	end
	tenmatMMd = tensor2tenmat(tensorMMd,[1,2],[3,4])
	return tenmatMMd
end


text:
mtg = read_mtg("files/simple_plant.mtg", Dict);

@testset "DataFrame" begin
    df_mtg = DataFrame(mtg, [:scales, :Length])
    @test df_mtg[1, :scales] == [0, 1, 2, 3, 3]
    @test df_mtg[7, :Length] == 0.2
end

@testset "MetaGraph" begin
    meta_mtg = MetaGraph(mtg)
    @test meta_mtg[1][:scales] == [0, 1, 2, 3, 3]
    @test meta_mtg[7][:Length] == 0.2
end


text:
using Compose
import Cairo, Fontconfig

set_default_graphic_size(6.6inch, 3.3inch)

# See Picasso's "dog" sketch
dog = [[(183, 268), (186, 256), (189, 244)], [(290, 244), (300, 230), (339, 245)], [(350,290), (360, 300), (355, 210)],
[(370, 207), (380,196), (375, 193)], [(310, 220), (190, 220), (164, 205)], [(135, 194), (135, 265), (153, 275)],
[(168, 275), (170, 180), (150, 190)], [(122, 214), (142, 204), (85, 240)], [(100, 247), (125, 233), (140, 238)]]

ub = UnitBox(0,0, 500,500)
p = compose(context(), stroke("black"), fillopacity(0.1),
    (context(0,0, 0.5,1, units=ub,  rotation=Rotation(-π/8)), bezigon((180, 280), dog)),
    (context(0.5,0, 0.5,1, units=ub, rotation=Rotation(π/8)),  bezigon([(180, 280)], [dog]))
)

imgs = [SVG("bezigon.svg"), PNG("bezigon.png")]

draw.(imgs, [p])


text:
module AkazeFeatures

export
    AKAZEOptions,
    AKAZE,
    Create_Nonlinear_Scale_Space,
    Feature_Detection,
    Compute_Main_Orientation,
    Compute_Descriptors,

    DIFFUSIVITY_TYPE,
    PM_G1,
    PM_G2,
    WEICKERT,
    CHARBONNIER,

    load_image_as_grayscale,
    original_akaze_features,
    dump_keypoints_text,
    plot_features

include("akaze-config.jl")
include("akaze-descriptor.jl")
include("akaze.jl")
include("fed.jl")
include("nonlinear-diffusion.jl")
include("utils.jl")

end


text:
using FlightSims
const FS = FlightSims

using Plots
using Transducers

function get_traj_data(env, Δt, tf)
    x0 = State(env)()
    prob, df = sim(
                   x0,
                   Dynamics!(env);
                   tf=tf,
                   savestep=Δt,
                  )
    df
end

function plot_figures!(fig, df, pos_cmd_func)
    ts = df.time
    poss = df.sol |> Map(datum -> datum.multicopter.state.p) |> collect
    poss_cmd = ts |> Map(pos_cmd_func) |> collect
    plot!(fig,
          ts, hcat(poss_cmd...)';
          label="desired",
          color=:red,
         )
    plot!(fig,
          ts, hcat(poss...)';
          label="true",
          color=:blue,
         )
end

function main()
    pos_cmd_func = (t) -> [sin(t), cos(t), t]
    env = FS.BacksteppingPositionController_StaticAllocator_MulticopterEnv(pos_cmd_func)
    Δt = 0.01
    tf = 10.0
    df = get_traj_data(env, Δt, tf)
    # plot
    fig = plot(; legend=:topleft)
    plot_figures!(fig, df, pos_cmd_func)
    display(fig)
end


text:
export load_graphml

"""
    load_graphml(filename)

Read a graphML format file with file name `filename`.
"""
function load_graphml(filename)
    open(EzXML.StreamReader, filename) do reader
        serial = 0
        nodes = []
        edges = Edge[]

        # scan nodes
        for typ in reader
            if typ == EzXML.READER_ELEMENT
                elname = EzXML.name(reader)
                if elname == "node"
                    push!(nodes, reader["id"])
                elseif elname == "edge"
                    e = Edge(reader["source"], reader["target"])
                    push!(edges, e)
                end
            end
        end
    return digraph(nodes, edges)
    end
end


text:
function readstat_get_file_label(metadata::Ptr{Nothing})
    ptr = ccall((:readstat_get_file_label, libreadstat), Cstring, (Ptr{Nothing},), metadata)
    return ptr == C_NULL ? "" : unsafe_string(ptr)
end

function readstat_get_modified_time(metadata::Ptr{Nothing})
    return ccall((:readstat_get_modified_time, libreadstat), UInt, (Ptr{Nothing},), metadata)
end

function readstat_get_file_format_version(metadata::Ptr{Nothing})
    return ccall((:readstat_get_file_format_version, libreadstat), Cint, (Ptr{Nothing},), metadata)
end

function readstat_get_row_count(metadata::Ptr{Nothing})
    return ccall((:readstat_get_row_count, libreadstat), Cint, (Ptr{Nothing},), metadata)
end

function readstat_get_var_count(metadata::Ptr{Nothing})
    return ccall((:readstat_get_var_count, libreadstat), Cint, (Ptr{Nothing},), metadata)
end

function readstat_value_is_missing(value::ReadStatValue, variable::Ptr{Nothing})
    return Bool(ccall((:readstat_value_is_missing, libreadstat), Cint, (ReadStatValue,Ptr{Nothing}), value, variable))
end

function readstat_variable_get_index(variable::Ptr{Nothing})
    return ccall((:readstat_variable_get_index, libreadstat), Cint, (Ptr{Nothing},), variable)
end

function readstat_variable_get_name(variable::Ptr{Nothing})
    return unsafe_string(ccall((:readstat_variable_get_name, libreadstat), Cstring, (Ptr{Nothing},), variable))
end

function readstat_variable_get_type(variable::Ptr{Nothing})
    return ccall((:readstat_variable_get_type, libreadstat), Cint, (Ptr{Nothing},), variable)
end

function readstat_variable_get_storage_width(variable::Ptr{Nothing})
    return ccall((:readstat_variable_get_storage_width, libreadstat), Csize_t, (Ptr{Nothing},), variable)
end

function readstat_variable_get_measure(variable::Ptr{Nothing})
    return ccall((:readstat_variable_get_measure, libreadstat), Cint, (Ptr{Nothing},), variable)
end

function readstat_variable_get_alignment(variable::Ptr{Nothing})
    return ccall((:readstat_variable_get_alignment, libreadstat), Cint, (Ptr{Nothing},), variable)
end

function readstat_parser_free(parser::Ptr{Nothing})
    return ccall((:readstat_parser_free, libreadstat), Nothing, (Ptr{Nothing},), parser)
end

function readstat_value_type(val::Value)
    return ccall((:readstat_value_type, libreadstat), Cint, (Value,), val)
end

function readstat_parse(filename::String, type::Val{:dta}, parser::Ptr{Nothing}, ds::ReadStatDataFrame)
    return ccall((:readstat_parse_dta, libreadstat), Cint, (Ptr{Nothing}, Cstring, Any), parser, string(filename), ds)
end

function readstat_parse(filename::String, type::Val{:sav}, parser::Ptr{Nothing}, ds::ReadStatDataFrame)
    return ccall((:readstat_parse_sav, libreadstat), Cint, (Ptr{Nothing}, Cstring, Any), parser, string(filename), ds)
end

function readstat_parse(filename::String, type::Val{:por}, parser::Ptr{Nothing}, ds::ReadStatDataFrame)
    return ccall((:readstat_parse_por, libreadstat), Cint, (Ptr{Nothing}, Cstring, Any), parser, string(filename), ds)
end

function readstat_parse(filename::String, type::Val{:sas7bdat}, parser::Ptr{Nothing}, ds::ReadStatDataFrame)
    return ccall((:readstat_parse_sas7bdat, libreadstat), Cint, (Ptr{Nothing}, Cstring, Any), parser, string(filename), ds)
end

function readstat_parse(filename::String, type::Val{:xport}, parser::Ptr{Nothing}, ds::ReadStatDataFrame)
    return ccall((:readstat_parse_xport, libreadstat), Cint, (Ptr{Nothing}, Cstring, Any), parser, string(filename), ds)
end

function readstat_variable_get_missing_ranges_count(variable::Ptr{Nothing})
    return ccall((:readstat_variable_get_missing_ranges_count, libreadstat), Cint, (Ptr{Nothing},), variable)
end


text:
function VelocityPrior(t) #historical model
    #return 10.0 .+ 4.0 .* sin(t ./ 10)
    return MaxDriverSpeed .+ 1.0 .* sin(t ./ 10)
end

function VariancePrior(t) #historical model
    return 5.25
end

function Deviation(v) #deviation function
    return v >= MaxDriverSpeed ? 1.0 : -1.0
end

#I want to integrate the driver position, which is prior+deviation

function DriverVelocity(t, gp = nothing)
    v = VelocityPrior(t)
    if gp == nothing
        devfcn = Deviation(v)
    else
        devfcn = LearnedDeviation(gp, v)
    end
    return v + devfcn[1][1]
end

function DriverUncertainty(t, gp)
    v = VelocityPrior(t)
    μ, Σ = predict_y(gp, [v])
    return Σ[1]
end

function DriverPosition(ti, TimeSamples, gp = nothing, tol = 1e-3)
    if isa(TimeSamples, Int) || isa(TimeSamples, Float64)
        integral, err =
            quadgk(x -> DriverVelocity(x, gp), ti, TimeSamples, rtol = tol)
        return integral
    else
        nsamples, npaths = size(TimeSamples)
        PosArray = zeros(nsamples, npaths)
        for p = 1:npaths
            Threads.@threads for i = 1:nsamples
                PosArray[i, p], err = quadgk(
                    x -> DriverVelocity(x, gp),
                    ti,
                    TimeSamples[i, p],
                    rtol = tol,
                )
            end
        end
    end
    return PosArray
end

function DriverUncertainty(ti, TimeSamples, gp, tol = 1e-1)
    if isa(TimeSamples, Int) || isa(TimeSamples, Float64)
        integral, err =
            quadgk(x -> DriverUncertainty(x, gp), ti, TimeSamples, rtol = tol)
        return integral
    else
        nsamples, npaths = size(TimeSamples)
        UArray = zeros(nsamples, npaths)
        for p = 1:npaths
            Threads.@threads for i = 1:nsamples
                UArray[i, p], err = quadgk(
                    x -> DriverUncertainty(x, gp),
                    ti,
                    TimeSamples[i, p],
                    rtol = tol,
                )
            end
        end
    end
    return UArray
end

function DriverPositionVector(tv = [0.0, 1.0], gp = nothing)
    l = length(tv) - 1
    pos = zeros(l)
    for i = 1:l
        pos[i] = DriverPosition(tv[1], tv[i+1], gp)
    end
    return pos
end

max_stars_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_issues_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_forks_repo_name: gbarsih/Multi-Path-Safe-Rendezvous


text:
function mse(semantics::Vector{T}, data) where T
    res = 0
    @simd for i in 1:length(semantics)
        res += (data.target[i] - semantics[i])^2
    end
    res /= length(data.target)
end

function mse(tree, data)
    res = 0
    for i in 1:length(data.target)
        output = evaluate(tree, view(data.input, :, i))
        res += (data.target[i] - output)^2
    end
    res /= length(data.target)
end

rmse(x, data) = sqrt(mse(x, data))

function se(tree, data)
    res = 0
    for i in 1:length(data.target)
        output = evaluate(tree, view(data.input, :, i))
        res += (data.target[i] - output)^2
    end
    res
end

function se(semantics::Vector{T}, data) where T
    res = 0
    @simd for i in 1:length(semantics)
        res += (data.target[i] - semantics[i])^2
    end
    res
end

max_stars_repo_name: kogad/JuliaGSGP

max_issues_repo_name: kogad/JuliaGSGP

max_forks_repo_name: kogad/JuliaGSGP


text:
function pythagorean_triplets(input)
    triplets = []
    for a in 1:(ceil(input/2))
        for b in a:(ceil((input-a)/2))
            c = input - a - b
            triplet = sort([a, b, c])
            if (triplet[1]^2 + triplet[2]^2 == triplet[3]^2)
                push!(triplets, Tuple(Int.(triplet)))
            end
        end
    end
    return(unique(triplets))
end


text:
_featurepositions(a::Integer, b::Integer) = a:3:b
function _featurepositions(a::AbstractString, b::AbstractString)
    a = Int[ parse(Int,x) for x in split(a,',') ]
    b = Int[ parse(Int,x) for x in split(b,',') ]
    vcat( map((x,y)->x:3:y, a, b)... )
end

function featurepositions(featureTable::DataFrame, featureName)
    row = findfirst(isequal(featureName), featureTable[!,:Feature])
    @assert row!=0
    _featurepositions(featureTable[row,:Start], featureTable[row,:End])
end

function featuresegment(featureTable::DataFrame, featureName)
    row = findfirst(isequal(featureName), featureTable[!,:Feature])
    @assert row!=0
    featureTable[row,:Segment]
end

# returns the most specific (i.e. shortest) feature covering the given position
function featureat(featureTable::DataFrame, position::Integer)
    ind = findall(featureTable[!,:Start].<=position.<=featureTable[!,:End]) # matching features
    isempty(ind) && return ""
    i = argmin( featureTable[ind,:End] - featureTable[ind,:Start] )
    featureTable[ind[i],:Feature] # name of shortest feature
end

function codonposition(featureTable::DataFrame, featureName, position::Integer)
    i = findfirst(isequal(featureName), featureTable[!,:Feature])
    @assert i>0 "Feature $featureName not found in table"
    @assert featureTable[i,:Start]<=position<=featureTable[i,:End]
    p = position-featureTable[i,:Start]
    @assert mod(p,3)==0 "position does does not match the ORF"
    div(p,3)+1
end

function annotatevariants!(swarm::AnnotatedArray, featureTable::DataFrame, reference::AbstractString)
    pos = swarm[:position]
    feature = map(p->featureat(featureTable,p), pos)
    refAA = codon2aa(seq2codons(reference))[pos]
    fpos = map((f,p)->codonposition(featureTable,f,p), feature, pos)
    variantAA = swarm[:aminoacid]

    sz = map(max, size(pos), size(variantAA)) # replace with broadcast_shape?
    # variant in style: feature-refAA-featurePos-variantAA
    variant = broadcast!( (f,r,p,v)->"$f-$r$p$v", Vector{String}(undef,sz), feature, refAA, fpos, variantAA)

    annotate!(swarm,:feature,feature)
    annotate!(swarm,:variant,variant)
end

function variantdescriptions(pos, codons, featureTable::DataFrame, reference::AbstractString)
    pos    = reshape(pos,1,length(pos))
    codons = reshape(codons,length(codons),1)

    feature = map(p->featureat(featureTable,p), pos)
    refAA = codon2aa(seq2codons(reference))[pos]
    fpos = map((f,p)->codonposition(featureTable,f,p), feature, pos)
    variantAA = codon2aa(codons)

    sz = map(max, size(pos), size(variantAA)) # replace with broadcast_shape?
    # variant in style: feature-refAA-featurePos-variantAA
    variants = broadcast!( (f,r,p,v)->"$f-$r$p$v", Vector{String}(undef,sz), feature, refAA, fpos, variantAA)
    variants[:]
end


text:
to_mat(A::AbstractMatrix) = A
function to_mat(y::AbstractVector)
    ymat = zeros(eltype(y), length(y)-1, 2);
    ymat[:,1] .= @view y[1:end-1];
    ymat[:,2] .= @view y[2:end]
    return ymat
end
function to_mat(pairs::Vector{Pair{T1,T2}} where {T1<:Tuple{<:Real,<:Real}, T2<:Union{Real, Tuple{<:Real,<:Real}}})
    xmat, ymat = zeros(length(pairs), 2), zeros(length(pairs), 2)
    for i = 1:length(pairs)
        xmat[i,:] .= pairs[i].first
        ymat[i,:] .= pairs[i].second
    end
    return xmat, ymat
end

to_toeplitz(vals) = to_toeplitz(vals, vals)
function to_toeplitz(vals_left::AbstractVector{T1}, vals_top::AbstractVector{T2}) where {T1,T2}
    @assert vals_left[1] == vals_top[1]
    N_rows, N_cols = length(vals_left), length(vals_top)
    A = zeros(promote_type(T1,T2), N_rows, N_cols)
    vals_total = zeros(promote_type(T1,T2), N_rows+N_cols-1)
    vals_total[1:N_cols] .= @view(vals_top[end:-1:1])
    vals_total[N_rows+1:end] .= @view(vals_left[2:end])
    for j = 1:N_cols
        A[:,j] .= @view(vals_total[end-N_rows+2-j:end+1-j])
    end
    return A
end

to_hankel(vals) = to_hankel(vals, vals)
function to_hankel(vals_left::AbstractVector{T1}, vals_bottom::AbstractVector{T2}) where {T1,T2}
    @assert vals_left[end] == vals_bottom[1]
    N_rows, N_cols = length(vals_left), length(vals_bottom)
    A = zeros(promote_type(T1,T2), N_rows, N_cols)
    vals_total = zeros(promote_type(T1,T2), N_rows+N_cols-1)
    vals_total[1:N_rows] .= vals_left
    vals_total[N_rows+1:end] .= @view(vals_bottom[2:end])
    for j = 1:N_cols
        A[:,j] .= @view(vals_total[j:j+N_rows-1])
    end
    return A
end

function nan2zero!(x::Array{T}) where T
    for idx = 1:length(x)
        if isnan(x[idx])
            x[idx] = zero(T)
        end
    end
end

function freqz(h::Vector{T}; fs::Real = 1, N = 1000) where T
    ω = range(0, stop = π, length = N)
    H = zeros(Complex{T}, length(ω))
    for (idx,h_n) in enumerate(h)
        n = idx-1
        H .+= h_n * exp.(-n*1im*ω)
    end
    return H, ω * (fs/2π)
end


text:
@enum RTError::Int32 RT_None = 0 RT_Debug = 1 RT_Warning = 2 RT_Failure = 3 RT_Fatal = 4
@enum RTIndexType::Int32 RT_RTree = 0 RT_MVRTree = 1 RT_TPRTree = 2 RT_InvalidIndexType = -99
@enum RTStorageType::Int32 RT_Memory = 0 RT_Disk = 1 RT_Custom = 2 RT_InvalidStorageType = -99
@enum RTIndexVariant::Int32 RT_Linear = 0 RT_Quadratic = 1 RT_Star = 2 RT_InvalidIndexVariant = -99

IndexH = Ptr{Nothing}
IndexItemH = Ptr{Nothing}
IndexPropertyH = Ptr{Nothing}


text:
ifftc(x;dims=[1,2])=fftshift(ifft(ifftshift(x,dims),dims),dims)*prod(size(x)[dims])
fftc(x;dims=[1,2]) =fftshift(fft(ifftshift(x,dims),dims),dims)

#DCF
"""
Hoge, R.D., Kwan, R.K.S. and Bruce Pike, G. (1997), Density compensation functions for spiral MRI. Magn. Reson. Med., 38: 117-128. https://doi.org/10.1002/mrm.1910380117"""
function dcf_nufft(seq)
    _, k = get_kspace(seq)
	t_adc = get_sample_times(seq) |> t-> convert(Vector{Float64},t)
	Gx, Gy, Gz = get_grads(seq, t_adc)
	dcf = abs.(Gx .* k[:,1] .+ Gy .* k[:,2] .+ Gz .* k[:,3]) # dcf = γ|<k(t), g(t)>|
	dcf ./ maximum(dcf[:])
end

# using MRIReco: AcquisitionData
# function reconstruction(signal, recParams)
#     #Param unpack
#     Nx =  get(recParams, "Nx",  101)				|> x->floor(Int,x)
# 	Ny =  get(recParams, "Ny",  recParams["Nx"]) 	|> x->floor(Int,x)
# 	Nz =  get(recParams, "Nz",  1) 					|> x->floor(Int,x)
# 	name = get(recParams, "Name", "")
#     recon = get(recParams, "recon", "MIRT")

# 	println("")
# 	@info "Running reconstruction..."
# 	@time begin
# 	#K-data, only 2D for now
# 	signal_phase = signal .* get_sample_phase_compensation(seq) #ADC phase compensation
# 	kdata = reshape(signal_phase, Nx,Ny,Nz)
# 	image = zeros(2,2,1)
# 	kspace = zeros(2,2,1)
# 	if recon != "skip"
# 		cartesian = !(occursin("rad", name) || occursin("spiral", name))
# 		if recon == "MRIReco"
# 			#MRIReco.jl, WIP
# 			#Getting kspace coordinates
# 			times = Koma.get_sample_times(seq)
# 			_, traj = Koma.get_kspace(seq)
# 			TE = 0
# 			AQ = maximum(seq[Koma.is_ADC_on.(seq)].ADC.T)
# 			tr = Trajectory(name, traj', times, TE, AQ, Ny, Nx, Nz, cartesian, !cartesian)
# 			#RawData
# 			sequenceInfo = Dict([Symbol(key) => value for (key,value) = seq.DEF]...)
# 			kdata = [zeros(ComplexF64,size(nodes,2),nc) for echo=1:1, slice=1:nz, rep=1:1]
# 			raw = AcquisitionData(sequenceInfo,tr,kdata)
# 			#(1. dim k-spacenodes, 2. dim coils) the outer dims describe:
# 			#1. dim echoes, 2. dim slices, 3. dim repetitions
# 			println(raw)
# 			# sequenceInfo::Dict{Symbol,Any}
# 			# traj::Vector{Trajectory}
# 			# kdata::Array{Matrix{ComplexF64},3}
# 			# subsampleIndices::Vector{Vector{Int64}}
# 			# encodingSize::Vector{Int64}
# 			# fov::Vector{Float64}
# 		end
# 		if recon=="MIRT"
# 			#MIRT.jl
# 			#Getting kspace coordinates
# 			_, traj = get_kspace(seq)
# 			k2d = k[:,1:2]
# 			#Normalize
# 			kmax = maximum(abs.(k2d[:]))
# 			ω = π * k2d / kmax
# 			#Number of samples and Operator definition, TODO: consider vector FOV vectors
# 			if non_cartesian
# 				N = 2 * round(Int, kmax * seq.DEF["FOV"] )
# 				Nu = Nv = N
# 				dcf = dcf_nufft(seq)
# 			else
# 				Nu = Nx
# 				Nv = Ny
# 				dcf = 1
# 			end
# 			image = zeros(ComplexF64,Nu,Nv,Nz)
# 			#MIRT.jl NUFFT
# 			for i = 1:Nz
# 				n_shift=round.(Int, [Nu/2; Nv/2])
# 				p = nufft_init(ω, (Nu,Nv); n_shift, pi_error=false)
# 				#Reconstruction
# 				kdata_i = kdata[:,:,i]
# 				print(size(kdata_i))
# 				image[:,:,i] = p.A' * (dcf .* kdata_i[:])
# 			end
# 			println(size(image))
# 			kspace = fftc(image, dims=[1,2])
# 		end
# 	end
# 	end
# 	(image, kspace)
# end


text:
"""
calculates the index in the time series corresponding to the data. Assumes that the dates vector is sorted.
"""
function find_timestamp_index(
    dates::Union{Vector{Dates.DateTime}, StepRange{Dates.DateTime, Dates.Millisecond}},
    date::Dates.DateTime,
)
    if date == first(dates)
        index = 1
    elseif date == last(dates)
        index = length(dates)
    else
        dates_resolution = dates[2] - dates[1]
        index = 1 + ((date - first(dates)) ÷ dates_resolution)
    end
    # Uncomment for debugging. The method below is fool proof but slower
    # s_index = findlast(dates .<= date)
    # IS.@assert_op index == s_index
    if index < 1 || index > length(dates)
        error("Requested timestamp $date not in the provided dates $dates")
    end
    return index
end


text:
# Based on https://github.com/JuliaGPU/CUDA.jl/blob/master/examples/peakflops.jl

_kernel_fma_nfmas()::Int = 100_000
_kernel_fma_N()::Int = Int((_kernel_fma_nfmas() - 1) ÷ 3)

using CUDA: i32
"Dummy kernel doing `_kernel_fma_nfmas()` many FMAs (default: `100_000`)."
function kernel_fma(a, b, c, out)
    i = (blockIdx().x - 1i32) * blockDim().x + threadIdx().x
    @inbounds if i <= length(out)
        a_val = a[i]
        b_val = b[i]
        c_val = c[i]

        # We suggest LLVM to unroll the loop. For large values, however,
        # there might still be a loop (with a large loop body of ~1000 fmas).
        # Investigate with `@device_code_ptx @cuda kernel_fma(d_a, d_b, d_c, d_out)`.
        @unroll for j in 1:_kernel_fma_N()
            a_val = CUDA.fma(a_val, b_val, c_val)
            b_val = CUDA.fma(a_val, b_val, c_val)
            c_val = CUDA.fma(a_val, b_val, c_val)
        end

        out[i] = CUDA.fma(a_val, b_val, c_val)
    end

    return nothing
end

"""
    peakflops_gpu_fmas(; size::Integer=5_000_000, dtype=Float32, nbench=5, nkernel=5, device=CUDA.device(), verbose=true)
Tries to estimate the peak performance of a GPU in TFLOP/s by measuring the time
it takes to perform `_kernel_fma_nfmas() * size` many FMAs on CUDA cores.

**Keyword arguments:**
* `device` (default: `CUDA.device()`): CUDA device to be used.
* `dtype` (default: `Float32`): element type of the matrices.
* `size` (default: `5_000_000`): length of vectors.
* `nkernel` (default: `5`): number of kernel calls that make up one benchmarking sample.
* `nbench` (default: `5`): number of measurements to be performed the best of which is used for the TFLOP/s computation.
* `verbose` (default: `true`): toggle printing.
* `io` (default: `stdout`): set the stream where the results should be printed.
"""
function peakflops_gpu_fmas(;
    size::Integer=5_000_000,
    dtype=Float32,
    nbench=5,
    nkernel=5,
    device::CuDevice=CUDA.device(),
    verbose=true,
    io::IO=stdout,
)
    device!(device) do
        d_a = CUDA.rand(dtype, size)
        d_b = CUDA.rand(dtype, size)
        d_c = CUDA.rand(dtype, size)
        d_out = CUDA.zeros(dtype, size)

        kernel = @cuda launch = false kernel_fma(d_a, d_b, d_c, d_out)
        config = launch_configuration(kernel.fun)
        threads = min(size, config.threads)
        blocks = cld(size, threads)

        # warm-up
        CUDA.@elapsed kernel(d_a, d_b, d_c, d_out)

        t = Inf
        for _ in 1:nbench
            Δt = CUDA.@elapsed begin
                for _ in 1:nkernel
                    kernel(d_a, d_b, d_c, d_out; threads=threads, blocks=blocks)
                end
            end
            t = min(t, Δt)
        end
        t /= nkernel

        flopcount = 2 * _kernel_fma_nfmas() * size
        flops = (flopcount * 1e-12) / t

        if verbose
            printstyled(io,"Peakflops (TFLOP/s):\n"; bold=true)
            print(io," └ max: ")
            printstyled(io, round(flops; digits=2), "\n"; color=:green, bold=true)
        end
        return flops
    end
end


text:
module OpenStreetMapGraphs

import Dates
import HTTP
import JSON
import LightGraphs
import MetaGraphs

export overpassquery, overpass
export oneway
export maxspeed

include("load_save.jl")
include("oneway.jl")
include("maxspeed.jl")

include("parse_tag.jl")

end # module


text:
function main()
    for val in rand(10_000)
        sol = rand(10000)
        sort!(sol)
    end
end

main()

max_stars_repo_name: dmolina/DaemonMode

max_issues_repo_name: dmolina/DaemonMode

max_forks_repo_name: dmolina/DaemonMode


text:
include("./params.jl")

function pprocess_gamma(lambda,k)
    # Example usage:
    # lambda = 20.0; # Spike-rate (Hz).
    # k = 2.3;       # Shape parameter.
    # x, rs = pprocess_gamma(lambda, k)
    #
    # Input argument
    # y:    Spike-rate.
    # k:    The shape parameter of the ISI distribution.
    #       * To use small k << 1, use small sampling time dt.
    #       * Maximum k allowed is ~10.
    #
    # Output argument
    # x:    Spike time.
    # rs:   Intensity function.

    #Initialize parameters

    x = [0.0]; # push!(x, num)
    rs = zeros(L);

    P = 0; R = 0;

    rs[1] = lambda;

    eta = - log(rand());
    for i = 2:L
        tau = i*dt - x[end]

        p = 1.0 / gamma(k) * (k * lambda).^k * tau.^(k-1) .* exp(-k*lambda*tau); # ISI
        P = P + p*dt;

        r = p / (1.0 - P); # Intensity
        R = R + r*dt;

        if eta <= R
            push!(x,i*dt) # Spike time

            P = 0.0; R = 0.0;
            eta = - log(rand());
        end
        rs[i] = r;
    end

    return x, rs
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
######################################################################
# greedy2.jl
#
# considers bearing only where we have to rotate
# Only assume bearing only sensor at this time
#
# Defines a greedy, info-theoretic policy
# Currently only considers a small set of the possible actions
######################################################################
# eval_type = 1 (discrete greedy maximization)
# eval_type = 2 (determinant)
# eval_type = 3 (largest eigenvalue)
struct GreedyPolicy2 <: Policy
    eval_type::Int
end

# loop over all actions.
# The one with smallest expected entropy is best
# We know this depends on sensor
function action(m::SearchDomain, v::Vehicle, o, f::DF, p::GreedyPolicy2)
    best_mi = -Inf
    best_x = v.x
    best_y = v.y

    mu = centroid(f)
    mu_x = mu[1]
    mu_y = mu[2]
    Omega = inv(covariance(f))
    a1 = Omega[1,1]
    b1 = Omega[1,2]
    c1 = Omega[2,1]
    d1 = Omega[2,2]
    mi = 0.

    for xi = 1:f.n, yi = 1:f.n
		x = (xi - 0.5) * f.cell_size
        y = (yi - 0.5) * f.cell_size

        if p.eval_type == 1
            xp = (x,y,0.0)
            mi = mutual_information2(f, xp)
        else
            xr = mu_x - x
            yr = mu_y - y

            q = 180.0 / (pi * (xr*xr + yr*yr))
            q = q*q / v.sensor.noise_sigma;

            sum_a = a1 + q*(yr*yr)
            sum_b = b1 - q*(xr*yr)
            sum_c = c1 - q*(xr*yr)
            sum_d = d1 + q*(xr*xr)
            if p.eval_type == 2
                mi = sum_a*sum_d - sum_b*sum_c
            else
                mi = smallest_eig(sum_a, sum_b, sum_c, sum_d)
            end
        end

        #println("mi = ", mi)
        if mi > best_mi
            best_mi = mi
            best_x = x
            best_y = y
		end
	end
	return (best_x-v.x, best_y-v.y, 0.0)
end

function smallest_eig(a::Float64, b::Float64, c::Float64, d::Float64)
	B = -1.0 * (d+a)
	C = a*d - b*c
	a1 = (-B + sqrt(B*B - 4C)) / 2.0
	a2 = (-B - sqrt(B*B - 4C)) / 2.0

	if a1 < a2
		return a1
	end
	return a2
end


text:
module FeatherFiles

using FeatherLib, IteratorInterfaceExtensions, TableTraits, TableTraitsUtils,
    DataValues, Arrow, FileIO, TableShowUtils
import IterableTables

export load, save, File, @format_str

include("missing-conversion.jl")

struct FeatherFile
    filename::String
end

function Base.show(io::IO, source::FeatherFile)
    TableShowUtils.printtable(io, getiterator(source), "Feather file")
end

function Base.show(io::IO, ::MIME"text/html", source::FeatherFile)
    TableShowUtils.printHTMLtable(io, getiterator(source))
end

Base.showable(::MIME"text/html", source::FeatherFile) = true

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::FeatherFile)
    TableShowUtils.printdataresource(io, getiterator(source))
end

Base.showable(::MIME"application/vnd.dataresource+json", source::FeatherFile) = true

function fileio_load(f::FileIO.File{FileIO.format"Feather"})
    return FeatherFile(f.filename)
end

IteratorInterfaceExtensions.isiterable(x::FeatherFile) = true
TableTraits.isiterabletable(x::FeatherFile) = true
# TableTraits.supports_get_columns_view(x::FeatherFile) = true
TableTraits.supports_get_columns_copy_using_missing(x::FeatherFile) = true

function IteratorInterfaceExtensions.getiterator(file::FeatherFile)
    rs = featherread(file.filename)

    for i=1:length(rs.columns)
        col_eltype = eltype(rs.columns[i])
        if isa(col_eltype, Union) && col_eltype.b <: Missing
            T = DataValueArrowVector{col_eltype.a,typeof(rs.columns[i])}
            rs.columns[i] = T(rs.columns[i])
        end
    end

    it = create_tableiterator(rs.columns, rs.names)

    return it
end

# function TableTraits.get_columns_view(file::FeatherFile)
#     rs = featherread(file.filename)

#     for i=1:length(rs.columns)
#         col_eltype = eltype(rs.columns[i])
#         if isa(col_eltype, Union) && col_eltype.b <: Missing
#             T = DataValueArrowVector{col_eltype.a,typeof(rs.columns[i])}
#             rs.columns[i] = T(rs.columns[i])
#         end
#     end

#     T = eval(:(@NT($(Symbol.(rs.names)...)))){typeof.(rs.columns)...}

#     return T(rs.columns...)
# end

function TableTraits.get_columns_copy_using_missing(file::FeatherFile)
     rs = featherread(file.filename)
     return NamedTuple{(Symbol.(rs.names)...,)}(((convert(Vector{eltype(c)}, c) for c in rs.columns)...,))
end

function fileio_save(f::FileIO.File{FileIO.format"Feather"}, data)
    isiterabletable(data) || error("Can't write this data to a Feather file.")

    columns, colnames = create_columns_from_iterabletable(data)

    columns = Any[c for c in columns]

    for i=1:length(columns)
        if eltype(columns[i]) <: DataValue
            T = MissingDataValueVector{eltype(eltype(columns[i])),typeof(columns[i])}
            columns[i] = T(columns[i])
        end
    end

    featherwrite(f.filename, columns, colnames)
end

end # module