text:
## array.jl: Dense arrays

typealias Vector{T} Array{T,1}
typealias Matrix{T} Array{T,2}
typealias VecOrMat{T} Union(Vector{T}, Matrix{T})

typealias DenseVector{T} DenseArray{T,1}
typealias DenseMatrix{T} DenseArray{T,2}
typealias DenseVecOrMat{T} Union(DenseVector{T}, DenseMatrix{T})

typealias StridedArray{T,N,A<:DenseArray} Union(DenseArray{T,N}, SubArray{T,N,A})
typealias StridedVector{T,A<:DenseArray}  Union(DenseArray{T,1}, SubArray{T,1,A})
typealias StridedMatrix{T,A<:DenseArray}  Union(DenseArray{T,2}, SubArray{T,2,A})
typealias StridedVecOrMat{T} Union(StridedVector{T}, StridedMatrix{T})

## Basic functions ##

size(a::Array) = arraysize(a)
size(a::Array, d) = arraysize(a, d)
size(a::Matrix) = (arraysize(a,1), arraysize(a,2))
length(a::Array) = arraylen(a)
function sizeof{T}(a::Array{T})
    (isbits(T) ? sizeof(eltype(a)) : sizeof(Ptr)) * length(a)
end

strides{T}(a::Array{T,1}) = (1,)
strides{T}(a::Array{T,2}) = (1, size(a,1))
strides{T}(a::Array{T,3}) = (1, size(a,1), size(a,1)*size(a,2))

isassigned(a::Array, i::Int...) = isdefined(a, i...)

## copy ##

function unsafe_copy!{T}(dest::Ptr{T}, src::Ptr{T}, N)
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Uint),
          dest, src, N*sizeof(T))
    return dest
end

function unsafe_copy!{T}(dest::Array{T}, dsto, src::Array{T}, so, N)
    if isbits(T)
        unsafe_copy!(pointer(dest, dsto), pointer(src, so), N)
    else
        for i=0:N-1
            @inbounds arrayset(dest, src[i+so], i+dsto)
        end
    end
    return dest
end

function copy!{T}(dest::Array{T}, dsto::Integer, src::Array{T}, so::Integer, N::Integer)
    if so+N-1 > length(src) || dsto+N-1 > length(dest) || dsto < 1 || so < 1
        throw(BoundsError())
    end
    unsafe_copy!(dest, dsto, src, so, N)
end

copy!{T}(dest::Array{T}, src::Array{T}) = copy!(dest, 1, src, 1, length(src))

function reinterpret{T,S}(::Type{T}, a::Array{S,1})
    nel = int(div(length(a)*sizeof(S),sizeof(T)))
    return reinterpret(T, a, (nel,))
end

function reinterpret{T,S}(::Type{T}, a::Array{S})
    if sizeof(S) != sizeof(T)
        error("result shape not specified")
    end
    reinterpret(T, a, size(a))
end

function reinterpret{T,S,N}(::Type{T}, a::Array{S}, dims::NTuple{N,Int})
    if !isbits(T)
        error("cannot reinterpret to type ", T)
    end
    if !isbits(S)
        error("cannot reinterpret Array of type ", S)
    end
    nel = div(length(a)*sizeof(S),sizeof(T))
    if prod(dims) != nel
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(nel)"))
    end
    ccall(:jl_reshape_array, Array{T,N}, (Any, Any, Any), Array{T,N}, a, dims)
end
reinterpret(t::Type,x) = reinterpret(t,[x])[1]

# reshaping to same # of dimensions
function reshape{T,N}(a::Array{T,N}, dims::NTuple{N,Int})
    if prod(dims) != length(a)
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))"))
    end
    if dims == size(a)
        return a
    end
    ccall(:jl_reshape_array, Array{T,N}, (Any, Any, Any), Array{T,N}, a, dims)
end

# reshaping to different # of dimensions
function reshape{T,N}(a::Array{T}, dims::NTuple{N,Int})
    if prod(dims) != length(a)
        throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(length(a))"))
    end
    ccall(:jl_reshape_array, Array{T,N}, (Any, Any, Any), Array{T,N}, a, dims)
end

## Constructors ##

similar(a::Array, T, dims::Dims)      = Array(T, dims)
similar{T}(a::Array{T,1})             = Array(T, size(a,1))
similar{T}(a::Array{T,2})             = Array(T, size(a,1), size(a,2))
similar{T}(a::Array{T,1}, dims::Dims) = Array(T, dims)
similar{T}(a::Array{T,1}, m::Int)     = Array(T, m)
similar{T}(a::Array{T,1}, S)          = Array(S, size(a,1))
similar{T}(a::Array{T,2}, dims::Dims) = Array(T, dims)
similar{T}(a::Array{T,2}, m::Int)     = Array(T, m)
similar{T}(a::Array{T,2}, S)          = Array(S, size(a,1), size(a,2))

# T[x...] constructs Array{T,1}
function getindex(T::NonTupleType, vals...)
    a = Array(T,length(vals))
    for i = 1:length(vals)
        a[i] = vals[i]
    end
    return a
end

getindex(T::(Type...)) = Array(T,0)

# T[a:b] and T[a:s:b] also contruct typed ranges
function getindex{T<:Number}(::Type{T}, r::Range)
    copy!(Array(T,length(r)), r)
end

function getindex{T<:Number}(::Type{T}, r1::Range, rs::Range...)
    a = Array(T,length(r1)+sum(length,rs))
    o = 1
    copy!(a, o, r1)
    o += length(r1)
    for r in rs
        copy!(a, o, r)
        o += length(r)
    end
    return a
end

function fill!{T<:Union(Int8,Uint8)}(a::Array{T}, x::Integer)
    ccall(:memset, Ptr{Void}, (Ptr{Void}, Int32, Csize_t), a, x, length(a))
    return a
end
function fill!{T<:Union(Integer,FloatingPoint)}(a::Array{T}, x)
    # note: preserve -0.0 for floats
    if isbits(T) && T<:Integer && convert(T,x) == 0
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Int32, Csize_t), a,0,length(a)*sizeof(T))
    else
        for i = 1:length(a)
            @inbounds a[i] = x
        end
    end
    return a
end

fill(v, dims::Dims)       = fill!(Array(typeof(v), dims), v)
fill(v, dims::Integer...) = fill!(Array(typeof(v), dims...), v)

for (fname, felt) in ((:zeros,:zero), (:ones,:one))
    @eval begin
        ($fname){T}(::Type{T}, dims...)  = fill!(Array(T, dims...), ($felt)(T))
        ($fname)(dims...)                = fill!(Array(Float64, dims...), ($felt)(Float64))
        ($fname){T}(x::AbstractArray{T}) = ($fname)(T, size(x))
    end
end

function eye(T::Type, m::Integer, n::Integer)
    a = zeros(T,m,n)
    for i = 1:min(m,n)
        a[i,i] = one(T)
    end
    return a
end
eye(m::Integer, n::Integer) = eye(Float64, m, n)
eye(T::Type, n::Integer) = eye(T, n, n)
eye(n::Integer) = eye(Float64, n)
eye{T}(x::AbstractMatrix{T}) = eye(T, size(x, 1), size(x, 2))

function one{T}(x::AbstractMatrix{T})
    m,n = size(x)
    m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    eye(T, m)
end

linspace(start::Integer, stop::Integer, n::Integer) =
    linspace(float(start), float(stop), n)
function linspace(start::Real, stop::Real, n::Integer)
    (start, stop) = promote(start, stop)
    T = typeof(start)
    a = Array(T, int(n))
    if n == 1
        a[1] = start
        return a
    end
    n -= 1
    S = promote_type(T, Float64)
    for i=0:n
        a[i+1] = start*(convert(S, (n-i))/n) + stop*(convert(S, i)/n)
    end
    a
end
linspace(start::Real, stop::Real) = linspace(start, stop, 100)

logspace(start::Real, stop::Real, n::Integer) = 10.^linspace(start, stop, n)
logspace(start::Real, stop::Real) = logspace(start, stop, 50)

## Conversions ##

convert{T,n}(::Type{Array{T}}, x::Array{T,n}) = x
convert{T,n}(::Type{Array{T,n}}, x::Array{T,n}) = x
convert{T,n,S}(::Type{Array{T}}, x::Array{S,n}) = convert(Array{T,n}, x)
convert{T,n,S}(::Type{Array{T,n}}, x::Array{S,n}) = copy!(similar(x,T), x)

function collect(T::Type, itr)
    if applicable(length, itr)
        # when length() isn't defined this branch might pollute the
        # type of the other.
        a = Array(T,length(itr)::Integer)
        i = 0
        for x in itr
            a[i+=1] = x
        end
    else
        a = Array(T,0)
        for x in itr
            push!(a,x)
        end
    end
    return a
end

collect(itr) = collect(eltype(itr), itr)

## Indexing: getindex ##

getindex(a::Array) = arrayref(a,1)

getindex(A::Array, i0::Real) = arrayref(A,to_index(i0))
getindex(A::Array, i0::Real, i1::Real) = arrayref(A,to_index(i0),to_index(i1))
getindex(A::Array, i0::Real, i1::Real, i2::Real) =
    arrayref(A,to_index(i0),to_index(i1),to_index(i2))
getindex(A::Array, i0::Real, i1::Real, i2::Real, i3::Real) =
    arrayref(A,to_index(i0),to_index(i1),to_index(i2),to_index(i3))
getindex(A::Array, i0::Real, i1::Real, i2::Real, i3::Real,  i4::Real) =
    arrayref(A,to_index(i0),to_index(i1),to_index(i2),to_index(i3),to_index(i4))
getindex(A::Array, i0::Real, i1::Real, i2::Real, i3::Real,  i4::Real, i5::Real) =
    arrayref(A,to_index(i0),to_index(i1),to_index(i2),to_index(i3),to_index(i4),to_index(i5))

getindex(A::Array, i0::Real, i1::Real, i2::Real, i3::Real,  i4::Real, i5::Real, I::Real...) =
    arrayref(A,to_index(i0),to_index(i1),to_index(i2),to_index(i3),to_index(i4),to_index(i5),to_index(I)...)

# Fast copy using copy! for UnitRange
function getindex(A::Array, I::UnitRange{Int})
    lI = length(I)
    X = similar(A, lI)
    if lI > 0
        copy!(X, 1, A, first(I), lI)
    end
    return X
end

function getindex{T<:Real}(A::Array, I::AbstractVector{T})
    return [ A[i] for i in to_index(I) ]
end
function getindex{T<:Real}(A::Range, I::AbstractVector{T})
    return [ A[i] for i in to_index(I) ]
end
function getindex(A::Range, I::AbstractVector{Bool})
    checkbounds(A, I)
    return [ A[i] for i in to_index(I) ]
end

# logical indexing

function getindex_bool_1d(A::Array, I::AbstractArray{Bool})
    checkbounds(A, I)
    n = sum(I)
    out = similar(A, n)
    c = 1
    for i = 1:length(I)
        if I[i]
            out[c] = A[i]
            c += 1
        end
    end
    out
end

getindex(A::Vector, I::AbstractVector{Bool}) = getindex_bool_1d(A, I)
getindex(A::Vector, I::AbstractArray{Bool}) = getindex_bool_1d(A, I)
getindex(A::Array, I::AbstractVector{Bool}) = getindex_bool_1d(A, I)
getindex(A::Array, I::AbstractArray{Bool}) = getindex_bool_1d(A, I)

## Indexing: setindex! ##
setindex!{T}(A::Array{T}, x) = arrayset(A, convert(T,x), 1)

setindex!{T}(A::Array{T}, x, i0::Real) = arrayset(A, convert(T,x), to_index(i0))
setindex!{T}(A::Array{T}, x, i0::Real, i1::Real) =
    arrayset(A, convert(T,x), to_index(i0), to_index(i1))
setindex!{T}(A::Array{T}, x, i0::Real, i1::Real, i2::Real) =
    arrayset(A, convert(T,x), to_index(i0), to_index(i1), to_index(i2))
setindex!{T}(A::Array{T}, x, i0::Real, i1::Real, i2::Real, i3::Real) =
    arrayset(A, convert(T,x), to_index(i0), to_index(i1), to_index(i2), to_index(i3))
setindex!{T}(A::Array{T}, x, i0::Real, i1::Real, i2::Real, i3::Real, i4::Real) =
    arrayset(A, convert(T,x), to_index(i0), to_index(i1), to_index(i2), to_index(i3), to_index(i4))
setindex!{T}(A::Array{T}, x, i0::Real, i1::Real, i2::Real, i3::Real, i4::Real, i5::Real) =
    arrayset(A, convert(T,x), to_index(i0), to_index(i1), to_index(i2), to_index(i3), to_index(i4), to_index(i5))
setindex!{T}(A::Array{T}, x, i0::Real, i1::Real, i2::Real, i3::Real, i4::Real, i5::Real, I::Real...) =
    arrayset(A, convert(T,x), to_index(i0), to_index(i1), to_index(i2), to_index(i3), to_index(i4), to_index(i5), to_index(I)...)

function setindex!{T<:Real}(A::Array, x, I::AbstractVector{T})
    for i in I
        A[i] = x
    end
    return A
end

function setindex!{T}(A::Array{T}, X::Array{T}, I::UnitRange{Int})
    if length(X) != length(I)
        throw_setindex_mismatch(X, (I,))
    end
    copy!(A, first(I), X, 1, length(I))
    return A
end

function setindex!{T<:Real}(A::Array, X::AbstractArray, I::AbstractVector{T})
    if length(X) != length(I)
        throw_setindex_mismatch(X, (I,))
    end
    count = 1
    if is(X,A)
        X = copy(X)
    end
    for i in I
        A[i] = X[count]
        count += 1
    end
    return A
end

# logical indexing

function assign_bool_scalar_1d!(A::Array, x, I::AbstractArray{Bool})
    checkbounds(A, I)
    for i = 1:length(I)
        if I[i]
            A[i] = x
        end
    end
    A
end

function assign_bool_vector_1d!(A::Array, X::AbstractArray, I::AbstractArray{Bool})
    checkbounds(A, I)
    c = 1
    for i = 1:length(I)
        if I[i]
            A[i] = X[c]
            c += 1
        end
    end
    if length(X) != c-1
        throw(DimensionMismatch("assigned $(length(X)) elements to length $(c-1) destination"))
    end
    A
end

setindex!(A::Array, X::AbstractArray, I::AbstractVector{Bool}) = assign_bool_vector_1d!(A, X, I)
setindex!(A::Array, X::AbstractArray, I::AbstractArray{Bool}) = assign_bool_vector_1d!(A, X, I)
setindex!(A::Array, x, I::AbstractVector{Bool}) = assign_bool_scalar_1d!(A, x, I)
setindex!(A::Array, x, I::AbstractArray{Bool}) = assign_bool_scalar_1d!(A, x, I)

# efficiently grow an array

function _growat!(a::Vector, i::Integer, delta::Integer)
    n = length(a)
    if i < div(n,2)
        _growat_beg!(a, i, delta)
    else
        _growat_end!(a, i, delta)
    end
    return a
end

function _growat_beg!(a::Vector, i::Integer, delta::Integer)
    ccall(:jl_array_grow_beg, Void, (Any, Uint), a, delta)
    @inbounds for k = 1:i-1
        a[k] = a[k+delta]
    end
    return a
end

function _growat_end!(a::Vector, i::Integer, delta::Integer)
    ccall(:jl_array_grow_end, Void, (Any, Uint), a, delta)
    n = length(a)
    @inbounds for k = n:-1:(i+delta)
        a[k] = a[k-delta]
    end
    return a
end

# efficiently delete part of an array

function _deleteat!(a::Vector, i::Integer, delta::Integer)
    n = length(a)
    last = i+delta-1
    if i-1 < n-last
        _deleteat_beg!(a, i, delta)
    else
        _deleteat_end!(a, i, delta)
    end
    return a
end

function _deleteat_beg!(a::Vector, i::Integer, delta::Integer)
    @inbounds for k = i+delta-1:-1:1+delta
        a[k] = a[k-delta]
    end
    ccall(:jl_array_del_beg, Void, (Any, Uint), a, delta)
    return a
end

function _deleteat_end!(a::Vector, i::Integer, delta::Integer)
    n = length(a)
    @inbounds for k = i:n-delta
        a[k] = a[k+delta]
    end
    ccall(:jl_array_del_end, Void, (Any, Uint), a, delta)
    return a
end

## Dequeue functionality ##

const _grow_none_errmsg =
    "[] cannot grow. Instead, initialize the array with \"T[]\", where T is the desired element type."

function push!{T}(a::Array{T,1}, item)
    if is(T,None)
        error(_grow_none_errmsg)
    end
    # convert first so we don't grow the array if the assignment won't work
    item = convert(T, item)
    ccall(:jl_array_grow_end, Void, (Any, Uint), a, 1)
    a[end] = item
    return a
end

function push!(a::Array{Any,1}, item::ANY)
    ccall(:jl_array_grow_end, Void, (Any, Uint), a, 1)
    arrayset(a, item, length(a))
    return a
end

function append!{T}(a::Array{T,1}, items::AbstractVector)
    if is(T,None)
        error(_grow_none_errmsg)
    end
    n = length(items)
    ccall(:jl_array_grow_end, Void, (Any, Uint), a, n)
    copy!(a, length(a)-n+1, items, 1, n)
    return a
end

function prepend!{T}(a::Array{T,1}, items::AbstractVector)
    if is(T,None)
        error(_grow_none_errmsg)
    end
    n = length(items)
    ccall(:jl_array_grow_beg, Void, (Any, Uint), a, n)
    if a === items
        copy!(a, 1, items, n+1, n)
    else
        copy!(a, 1, items, 1, n)
    end
    return a
end

function resize!(a::Vector, nl::Integer)
    l = length(a)
    if nl > l
        ccall(:jl_array_grow_end, Void, (Any, Uint), a, nl-l)
    else
        if nl < 0
            throw(BoundsError())
        end
        ccall(:jl_array_del_end, Void, (Any, Uint), a, l-nl)
    end
    return a
end

function sizehint(a::Vector, sz::Integer)
    ccall(:jl_array_sizehint, Void, (Any, Uint), a, sz)
    a
end

function pop!(a::Vector)
    if isempty(a)
        error("array must be non-empty")
    end
    item = a[end]
    ccall(:jl_array_del_end, Void, (Any, Uint), a, 1)
    return item
end

function unshift!{T}(a::Array{T,1}, item)
    if is(T,None)
        error(_grow_none_errmsg)
    end
    item = convert(T, item)
    ccall(:jl_array_grow_beg, Void, (Any, Uint), a, 1)
    a[1] = item
    return a
end

function shift!(a::Vector)
    if isempty(a)
        error("array must be non-empty")
    end
    item = a[1]
    ccall(:jl_array_del_beg, Void, (Any, Uint), a, 1)
    return item
end

function insert!{T}(a::Array{T,1}, i::Integer, item)
    if i < 1
        throw(BoundsError())
    end
    item = convert(T, item)
    n = length(a)
    if i > n
        ccall(:jl_array_grow_end, Void, (Any, Uint), a, i-n)
    else
        _growat!(a, i, 1)
    end
    a[i] = item
    return a
end

function deleteat!(a::Vector, i::Integer)
    if !(1 <= i <= length(a))
        throw(BoundsError())
    end
    return _deleteat!(a, i, 1)
end

function deleteat!{T<:Integer}(a::Vector, r::UnitRange{T})
    n = length(a)
    f = first(r)
    l = last(r)
    if !(1 <= f && l <= n)
        throw(BoundsError())
    end
    return _deleteat!(a, f, length(r))
end

function deleteat!(a::Vector, inds)
    n = length(a)
    s = start(inds)
    done(inds, s) && return a
    (p, s) = next(inds, s)
    q = p+1
    while !done(inds, s)
        (i,s) = next(inds, s)
        if !(q <= i <= n)
            i < q && error("indices must be unique and sorted")
            throw(BoundsError())
        end
        while q < i
            @inbounds a[p] = a[q]
            p += 1; q += 1
        end
        q = i+1
    end
    while q <= n
        @inbounds a[p] = a[q]
        p += 1; q += 1
    end
    ccall(:jl_array_del_end, Void, (Any, Uint), a, n-p+1)
    return a
end

const _default_splice = []

function splice!(a::Vector, i::Integer, ins::AbstractArray=_default_splice)
    v = a[i]
    m = length(ins)
    if m == 0
        _deleteat!(a, i, 1)
    elseif m == 1
        a[i] = ins[1]
    else
        _growat!(a, i, m-1)
        for k = 1:m
            a[i+k-1] = ins[k]
        end
    end
    return v
end

function splice!{T<:Integer}(a::Vector, r::UnitRange{T}, ins::AbstractArray=_default_splice)
    v = a[r]
    m = length(ins)
    if m == 0
        deleteat!(a, r)
        return v
    end

    n = length(a)
    f = first(r)
    l = last(r)
    d = length(r)

    if m < d
        delta = d - m
        if f-1 < n-l
            _deleteat_beg!(a, f, delta)
        else
            _deleteat_end!(a, l-delta+1, delta)
        end
    elseif m > d
        delta = m - d
        if f-1 < n-l
            _growat_beg!(a, f, delta)
        else
            _growat_end!(a, l+1, delta)
        end
    end

    for k = 1:m
        a[f+k-1] = ins[k]
    end
    return v
end

function empty!(a::Vector)
    ccall(:jl_array_del_end, Void, (Any, Uint), a, length(a))
    return a
end

## Unary operators ##

function conj!{T<:Number}(A::AbstractArray{T})
    for i=1:length(A)
        A[i] = conj(A[i])
    end
    return A
end

for f in (:-, :~, :conj, :sign)
    @eval begin
        function ($f)(A::StridedArray)
            F = similar(A)
            for i=1:length(A)
                F[i] = ($f)(A[i])
            end
            return F
        end
    end
end

(-)(A::StridedArray{Bool}) = reshape([ -A[i] for i=1:length(A) ], size(A))

real(A::StridedArray) = reshape([ real(x) for x in A ], size(A))
imag(A::StridedArray) = reshape([ imag(x) for x in A ], size(A))
real{T<:Real}(x::StridedArray{T}) = x
imag{T<:Real}(x::StridedArray{T}) = zero(x)

function !(A::StridedArray{Bool})
    F = similar(A)
    for i=1:length(A)
        F[i] = !A[i]
    end
    return F
end

## Binary arithmetic operators ##

promote_array_type{Scalar, Arry}(::Type{Scalar}, ::Type{Arry}) = promote_type(Scalar, Arry)
promote_array_type{S<:Real, A<:FloatingPoint}(::Type{S}, ::Type{A}) = A
promote_array_type{S<:Union(Complex, Real), AT<:FloatingPoint}(::Type{S}, ::Type{Complex{AT}}) = Complex{AT}
promote_array_type{S<:Integer, A<:Integer}(::Type{S}, ::Type{A}) = A
promote_array_type{S<:Integer}(::Type{S}, ::Type{Bool}) = S

./{T<:Integer}(x::Integer, y::StridedArray{T}) =
    reshape([ x    ./ y[i] for i=1:length(y) ], size(y))
./{T<:Integer}(x::StridedArray{T}, y::Integer) =
    reshape([ x[i] ./ y    for i=1:length(x) ], size(x))

./{T<:Integer}(x::Integer, y::StridedArray{Complex{T}}) =
    reshape([ x    ./ y[i] for i=1:length(y) ], size(y))
./{T<:Integer}(x::StridedArray{Complex{T}}, y::Integer) =
    reshape([ x[i] ./ y    for i=1:length(x) ], size(x))
./{S<:Integer,T<:Integer}(x::Complex{S}, y::StridedArray{T}) =
    reshape([ x    ./ y[i] for i=1:length(y) ], size(y))
./{S<:Integer,T<:Integer}(x::StridedArray{S}, y::Complex{T}) =
    reshape([ x[i] ./ y    for i=1:length(x) ], size(x))

# ^ is difficult, since negative exponents give a different type

.^(x::Number, y::StridedArray) =
    reshape([ x ^ y[i] for i=1:length(y) ], size(y))

.^(x::StridedArray, y::Number      ) =
    reshape([ x[i] ^ y for i=1:length(x) ], size(x))

for f in (:+, :-, :div, :mod, :&, :|, :$)
    @eval begin
        function ($f){S,T}(A::StridedArray{S}, B::StridedArray{T})
            F = similar(A, promote_type(S,T), promote_shape(size(A),size(B)))
            for i=1:length(A)
                @inbounds F[i] = ($f)(A[i], B[i])
            end
            return F
        end
        # interaction with Ranges
        function ($f){S,T<:Real}(A::StridedArray{S}, B::Range{T})
            F = similar(A, promote_type(S,T), promote_shape(size(A),size(B)))
            i = 1
            for b in B
                @inbounds F[i] = ($f)(A[i], b)
                i += 1
            end
            return F
        end
        function ($f){S<:Real,T}(A::Range{S}, B::StridedArray{T})
            F = similar(B, promote_type(S,T), promote_shape(size(A),size(B)))
            i = 1
            for a in A
                @inbounds F[i] = ($f)(a, B[i])
                i += 1
            end
            return F
        end
    end
end
for f in (:.+, :.-, :.*, :./, :.\, :.%, :div, :mod, :rem, :&, :|, :$)
    @eval begin
        function ($f){T}(A::Number, B::StridedArray{T})
            F = similar(B, promote_array_type(typeof(A),T))
            for i=1:length(B)
                @inbounds F[i] = ($f)(A, B[i])
            end
            return F
        end
        function ($f){T}(A::StridedArray{T}, B::Number)
            F = similar(A, promote_array_type(typeof(B),T))
            for i=1:length(A)
                @inbounds F[i] = ($f)(A[i], B)
            end
            return F
        end
    end
end

# functions that should give an Int result for Bool arrays
for f in (:.+, :.-)
    @eval begin
        function ($f)(A::Bool, B::StridedArray{Bool})
            F = similar(B, Int, size(B))
            for i=1:length(B)
                @inbounds F[i] = ($f)(A, B[i])
            end
            return F
        end
        function ($f)(A::StridedArray{Bool}, B::Bool)
            F = similar(A, Int, size(A))
            for i=1:length(A)
                @inbounds F[i] = ($f)(A[i], B)
            end
            return F
        end
    end
end
for f in (:+, :-)
    @eval begin
        function ($f)(A::StridedArray{Bool}, B::StridedArray{Bool})
            F = similar(A, Int, promote_shape(size(A), size(B)))
            for i=1:length(A)
                @inbounds F[i] = ($f)(A[i], B[i])
            end
            return F
        end
    end
end

## promotion to complex ##

function complex{S<:Real,T<:Real}(A::Array{S}, B::Array{T})
    if size(A) != size(B); throw(DimensionMismatch("")); end
    F = similar(A, typeof(complex(zero(S),zero(T))))
    for i=1:length(A)
        @inbounds F[i] = complex(A[i], B[i])
    end
    return F
end

function complex{T<:Real}(A::Real, B::Array{T})
    F = similar(B, typeof(complex(A,zero(T))))
    for i=1:length(B)
        @inbounds F[i] = complex(A, B[i])
    end
    return F
end

function complex{T<:Real}(A::Array{T}, B::Real)
    F = similar(A, typeof(complex(zero(T),B)))
    for i=1:length(A)
        @inbounds F[i] = complex(A[i], B)
    end
    return F
end

# use memcmp for lexcmp on byte arrays
function lexcmp(a::Array{Uint8,1}, b::Array{Uint8,1})
    c = ccall(:memcmp, Int32, (Ptr{Uint8}, Ptr{Uint8}, Uint),
              a, b, min(length(a),length(b)))
    c < 0 ? -1 : c > 0 ? +1 : cmp(length(a),length(b))
end

## data movement ##

function slicedim(A::Array, d::Integer, i::Integer)
    d_in = size(A)
    leading = d_in[1:(d-1)]
    d_out = tuple(leading..., 1, d_in[(d+1):end]...)

    M = prod(leading)
    N = length(A)
    stride = M * d_in[d]

    B = similar(A, d_out)
    index_offset = 1 + (i-1)*M

    l = 1

    if M==1
        for j=0:stride:(N-stride)
            B[l] = A[j + index_offset]
            l += 1
        end
    else
        for j=0:stride:(N-stride)
            offs = j + index_offset
            for k=0:(M-1)
                B[l] = A[offs + k]
                l += 1
            end
        end
    end
    return B
end

function flipdim{T}(A::Array{T}, d::Integer)
    nd = ndims(A)
    sd = d > nd ? 1 : size(A, d)
    if sd == 1 || isempty(A)
        return copy(A)
    end

    B = similar(A)

    nnd = 0
    for i = 1:nd
        nnd += int(size(A,i)==1 || i==d)
    end
    if nnd==nd
        # flip along the only non-singleton dimension
        for i = 1:sd
            B[i] = A[sd+1-i]
        end
        return B
    end

    d_in = size(A)
    leading = d_in[1:(d-1)]
    M = prod(leading)
    N = length(A)
    stride = M * sd

    if M==1
        for j = 0:stride:(N-stride)
            for i = 1:sd
                ri = sd+1-i
                B[j + ri] = A[j + i]
            end
        end
    else
        if isbits(T) && M>200
            for i = 1:sd
                ri = sd+1-i
                for j=0:stride:(N-stride)
                    offs = j + 1 + (i-1)*M
                    boffs = j + 1 + (ri-1)*M
                    copy!(B, boffs, A, offs, M)
                end
            end
        else
            for i = 1:sd
                ri = sd+1-i
                for j=0:stride:(N-stride)
                    offs = j + 1 + (i-1)*M
                    boffs = j + 1 + (ri-1)*M
                    for k=0:(M-1)
                        B[boffs + k] = A[offs + k]
                    end
                end
            end
        end
    end
    return B
end

function rotl90(A::StridedMatrix)
    m,n = size(A)
    B = similar(A,(n,m))
    for i=1:m, j=1:n
        B[n-j+1,i] = A[i,j]
    end
    return B
end
function rotr90(A::StridedMatrix)
    m,n = size(A)
    B = similar(A,(n,m))
    for i=1:m, j=1:n
        B[j,m-i+1] = A[i,j]
    end
    return B
end
function rot180(A::StridedMatrix)
    m,n = size(A)
    B = similar(A)
    for i=1:m, j=1:n
        B[m-i+1,n-j+1] = A[i,j]
    end
    return B
end
function rotl90(A::AbstractMatrix, k::Integer)
    k = mod(k, 4)
    k == 1 ? rotl90(A) :
    k == 2 ? rot180(A) :
    k == 3 ? rotr90(A) : copy(A)
end
rotr90(A::AbstractMatrix, k::Integer) = rotl90(A,-k)
rot180(A::AbstractMatrix, k::Integer) = mod(k, 2) == 1 ? rot180(A) : copy(A)

# note: probably should be StridedVector or AbstractVector
function reverse(A::AbstractVector, s=1, n=length(A))
    B = similar(A)
    for i = 1:s-1
        B[i] = A[i]
    end
    for i = s:n
        B[i] = A[n+s-i]
    end
    for i = n+1:length(A)
        B[i] = A[i]
    end
    B
end

reverse(v::StridedVector) = (n=length(v); [ v[n-i+1] for i=1:n ])
reverse(v::StridedVector, s, n=length(v)) = reverse!(copy(v), s, n)
function reverse!(v::StridedVector, s=1, n=length(v))
    r = n
    for i=s:div(s+n-1,2)
        v[i], v[r] = v[r], v[i]
        r -= 1
    end
    v
end

function vcat{T}(arrays::Array{T,1}...)
    n = 0
    for a in arrays
        n += length(a)
    end
    arr = Array(T, n)
    ptr = pointer(arr)
    offset = 0
    if isbits(T)
        elsz = sizeof(T)
    else
        elsz = div(WORD_SIZE,8)
    end
    for a in arrays
        nba = length(a)*elsz
        ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Uint),
              ptr+offset, a, nba)
        offset += nba
    end
    return arr
end

## find ##

# returns the index of the next non-zero element, or 0 if all zeros
function findnext(A, start::Integer)
    for i = start:length(A)
        if A[i] != 0
            return i
        end
    end
    return 0
end
findfirst(A) = findnext(A,1)

# returns the index of the next matching element
function findnext(A, v, start::Integer)
    for i = start:length(A)
        if A[i] == v
            return i
        end
    end
    return 0
end
findfirst(A,v) = findnext(A,v,1)

# returns the index of the next element for which the function returns true
function findnext(testf::Function, A, start::Integer)
    for i = start:length(A)
        if testf(A[i])
            return i
        end
    end
    return 0
end
findfirst(testf::Function, A) = findnext(testf, A, 1)

function find(testf::Function, A::StridedArray)
    # use a dynamic-length array to store the indexes, then copy to a non-padded
    # array for the return
    tmpI = Array(Int, 0)
    for i = 1:length(A)
        if testf(A[i])
            push!(tmpI, i)
        end
    end
    I = similar(A, Int, length(tmpI))
    copy!(I, tmpI)
    I
end

function find(A::StridedArray)
    nnzA = countnz(A)
    I = similar(A, Int, nnzA)
    count = 1
    for i=1:length(A)
        if A[i] != 0
            I[count] = i
            count += 1
        end
    end
    return I
end

find(x::Number) = x == 0 ? Array(Int,0) : [1]
find(testf::Function, x) = find(testf(x))

findn(A::AbstractVector) = find(A)

function findn(A::StridedMatrix)
    nnzA = countnz(A)
    I = similar(A, Int, nnzA)
    J = similar(A, Int, nnzA)
    count = 1
    for j=1:size(A,2), i=1:size(A,1)
        if A[i,j] != 0
            I[count] = i
            J[count] = j
            count += 1
        end
    end
    return (I, J)
end

function findnz{T}(A::StridedMatrix{T})
    nnzA = countnz(A)
    I = zeros(Int, nnzA)
    J = zeros(Int, nnzA)
    NZs = zeros(T, nnzA)
    count = 1
    if nnzA > 0
        for j=1:size(A,2), i=1:size(A,1)
            Aij = A[i,j]
            if Aij != 0
                I[count] = i
                J[count] = j
                NZs[count] = Aij
                count += 1
            end
        end
    end
    return (I, J, NZs)
end

function findmax(a)
    if isempty(a)
        error("array must be non-empty")
    end
    m = a[1]
    mi = 1
    for i=2:length(a)
        ai = a[i]
        if ai > m || m!=m
            m = ai
            mi = i
        end
    end
    return (m, mi)
end

function findmin(a)
    if isempty(a)
        error("array must be non-empty")
    end
    m = a[1]
    mi = 1
    for i=2:length(a)
        ai = a[i]
        if ai < m || m!=m
            m = ai
            mi = i
        end
    end
    return (m, mi)
end

indmax(a) = findmax(a)[2]
indmin(a) = findmin(a)[2]

# similar to Matlab's ismember
# returns a vector containing the highest index in b for each value in a that is a member of b
function indexin{T}(a::AbstractArray{T}, b::AbstractArray{T})
    bdict = Dict(b, 1:length(b))
    [get(bdict, i, 0) for i in a]
end

# findin (the index of intersection)
function findin(a, b::UnitRange)
    ind = Array(Int, 0)
    f = first(b)
    l = last(b)
    for i = 1:length(a)
        if f <= a[i] <= l
            push!(ind, i)
        end
    end
    ind
end

function findin(a, b)
    ind = Array(Int, 0)
    bset = union!(Set(), b)
    for i = 1:length(a)
        if in(a[i], bset)
            push!(ind, i)
        end
    end
    ind
end

# Copying subregions
function indcopy(sz::Dims, I::Vector)
    n = length(I)
    s = sz[n]
    for i = n+1:length(sz)
        s *= sz[i]
    end
    dst = eltype(I)[findin(I[i], i < n ? (1:sz[i]) : (1:s)) for i = 1:n]
    src = eltype(I)[I[i][findin(I[i], i < n ? (1:sz[i]) : (1:s))] for i = 1:n]
    dst, src
end

function indcopy(sz::Dims, I::(RangeIndex...))
    n = length(I)
    s = sz[n]
    for i = n+1:length(sz)
        s *= sz[i]
    end
    dst::typeof(I) = ntuple(n, i-> findin(I[i], i < n ? (1:sz[i]) : (1:s)))::typeof(I)
    src::typeof(I) = ntuple(n, i-> I[i][findin(I[i], i < n ? (1:sz[i]) : (1:s))])::typeof(I)
    dst, src
end

## Filter ##

# given a function returning a boolean and an array, return matching elements
filter(f::Function, As::AbstractArray) = As[map(f, As)::AbstractArray{Bool}]

function filter!(f::Function, a::Vector)
    insrt = 1
    for curr = 1:length(a)
        if f(a[curr])
            a[insrt] = a[curr]
            insrt += 1
        end
    end
    splice!(a, insrt:length(a))
    return a
end

function filter(f::Function, a::Vector)
    r = Array(eltype(a), 0)
    for i = 1:length(a)
        if f(a[i])
            push!(r, a[i])
        end
    end
    return r
end

## Transpose ##
const transposebaselength=64
function transpose!(B::StridedMatrix,A::StridedMatrix)
    m, n = size(A)
    size(B) == (n,m) || throw(DimensionMismatch("transpose"))

    if m*n<=4*transposebaselength
        @inbounds begin
            for j = 1:n
                for i = 1:m
                    B[j,i] = A[i,j]
                end
            end
        end
    else
        transposeblock!(B,A,m,n,0,0)
    end
    return B
end
function transposeblock!(B::StridedMatrix,A::StridedMatrix,m::Int,n::Int,offseti::Int,offsetj::Int)
    if m*n<=transposebaselength
        @inbounds begin
            for j = offsetj+(1:n)
                for i = offseti+(1:m)
                    B[j,i] = A[i,j]
                end
            end
        end
    elseif m>n
        newm=m>>1
        transposeblock!(B,A,newm,n,offseti,offsetj)
        transposeblock!(B,A,m-newm,n,offseti+newm,offsetj)
    else
        newn=n>>1
        transposeblock!(B,A,m,newn,offseti,offsetj)
        transposeblock!(B,A,m,n-newn,offseti,offsetj+newn)
    end
    return B
end
function ctranspose!(B::StridedMatrix,A::StridedMatrix)
    m, n = size(A)
    size(B) == (n,m) || throw(DimensionMismatch("transpose"))

    if m*n<=4*transposebaselength
        @inbounds begin
            for j = 1:n
                for i = 1:m
                    B[j,i] = conj(A[i,j])
                end
            end
        end
    else
        ctransposeblock!(B,A,m,n,0,0)
    end
    return B
end
function ctransposeblock!(B::StridedMatrix,A::StridedMatrix,m::Int,n::Int,offseti::Int,offsetj::Int)
    if m*n<=transposebaselength
        @inbounds begin
            for j = offsetj+(1:n)
                for i = offseti+(1:m)
                    B[j,i] = conj(A[i,j])
                end
            end
        end
    elseif m>n
        newm=m>>1
        ctransposeblock!(B,A,newm,n,offseti,offsetj)
        ctransposeblock!(B,A,m-newm,n,offseti+newm,offsetj)
    else
        newn=n>>1
        ctransposeblock!(B,A,m,newn,offseti,offsetj)
        ctransposeblock!(B,A,m,n-newn,offseti,offsetj+newn)
    end
    return B
end

function transpose(A::StridedMatrix)
    B = similar(A, size(A, 2), size(A, 1))
    transpose!(B, A)
end
function ctranspose(A::StridedMatrix)
    B = similar(A, size(A, 2), size(A, 1))
    ctranspose!(B, A)
end
ctranspose{T<:Real}(A::StridedVecOrMat{T}) = transpose(A)

transpose(x::StridedVector) = [ x[j] for i=1, j=1:size(x,1) ]
ctranspose{T}(x::StridedVector{T}) = T[ conj(x[j]) for i=1, j=1:size(x,1) ]

# set-like operators for vectors
# These are moderately efficient, preserve order, and remove dupes.

function intersect(v1, vs...)
    ret = Array(eltype(v1),0)
    for v_elem in v1
        inall = true
        for i = 1:length(vs)
            if !in(v_elem, vs[i])
                inall=false; break
            end
        end
        if inall
            push!(ret, v_elem)
        end
    end
    ret
end

function union(vs...)
    ret = Array(promote_eltype(vs...),0)
    seen = Set()
    for v in vs
        for v_elem in v
            if !in(v_elem, seen)
                push!(ret, v_elem)
                push!(seen, v_elem)
            end
        end
    end
    ret
end
# setdiff only accepts two args
function setdiff(a, b)
    args_type = promote_type(eltype(a), eltype(b))
    bset = Set(b)
    ret = Array(args_type,0)
    seen = Set()
    for a_elem in a
        if !in(a_elem, seen) && !in(a_elem, bset)
            push!(ret, a_elem)
            push!(seen, a_elem)
        end
    end
    ret
end
# symdiff is associative, so a relatively clean
# way to implement this is by using setdiff and union, and
# recursing. Has the advantage of keeping order, too, but
# not as fast as other methods that make a single pass and
# store counts with a Dict.
symdiff(a) = a
symdiff(a, b) = union(setdiff(a,b), setdiff(b,a))
symdiff(a, b, rest...) = symdiff(a, symdiff(b, rest...))

_cumsum_type{T<:Number}(v::AbstractArray{T}) = typeof(+zero(T))
_cumsum_type(v) = typeof(v[1]+v[1])

for (f, fp, op) = ((:cumsum, :cumsum_pairwise, :+),
                   (:cumprod, :cumprod_pairwise, :*) )
    # in-place cumsum of c = s+v(i1:n), using pairwise summation as for sum
    @eval function ($fp)(v::AbstractVector, c::AbstractVector, s, i1, n)
        if n < 128
            @inbounds c[i1] = ($op)(s, v[i1])
            for i = i1+1:i1+n-1
                @inbounds c[i] = $(op)(c[i-1], v[i])
            end
        else
            n2 = div(n,2)
            ($fp)(v, c, s, i1, n2)
            ($fp)(v, c, c[(i1+n2)-1], i1+n2, n-n2)
        end
    end

    @eval function ($f)(v::AbstractVector)
        n = length(v)
        c = $(op===:+ ? (:(similar(v,_cumsum_type(v)))) :
                        (:(similar(v))))
        if n == 0; return c; end
        ($fp)(v, c, $(op==:+ ? :(zero(v[1])) : :(one(v[1]))), 1, n)
        return c
    end

    @eval function ($f)(A::StridedArray, axis::Integer)
        dimsA = size(A)
        ndimsA = ndims(A)
        axis_size = dimsA[axis]
        axis_stride = 1
        for i = 1:(axis-1)
            axis_stride *= size(A,i)
        end

        B = $(op===:+ ? (:(similar(A,_cumsum_type(A)))) :
                        (:(similar(A))))

        if axis_size < 1
            return B
        end

        for i = 1:length(A)
            if div(i-1, axis_stride) % axis_size == 0
               B[i] = A[i]
            else
               B[i] = ($op)(B[i-axis_stride], A[i])
            end
        end

        return B
    end

    @eval ($f)(A::AbstractArray) = ($f)(A, 1)
end

for (f, op) = ((:cummin, :min), (:cummax, :max))
    @eval function ($f)(v::AbstractVector)
        n = length(v)
        cur_val = v[1]
        res = similar(v, n)
        res[1] = cur_val
        for i in 2:n
            cur_val = ($op)(v[i], cur_val)
            res[i] = cur_val
        end
        return res
    end

    @eval function ($f)(A::StridedArray, axis::Integer)
        dimsA = size(A)
        ndimsA = ndims(A)
        axis_size = dimsA[axis]
        axis_stride = 1
        for i = 1:(axis-1)
            axis_stride *= size(A,i)
        end

        if axis_size < 1
            return A
        end

        B = similar(A)

        for i = 1:length(A)
            if div(i-1, axis_stride) % axis_size == 0
               B[i] = A[i]
            else
               B[i] = ($op)(A[i], B[i-axis_stride])
            end
        end

        return B
    end

    @eval ($f)(A::AbstractArray) = ($f)(A, 1)
end

max_stars_repo_name: aplund/julia

max_issues_repo_name: aplund/julia

max_forks_repo_name: aplund/julia


text:
using ..Atmos
using ..Atmos: MoistureModel, thermo_state, turbulence_tensors
using ..Mesh.Topologies
using ..Mesh.Grids
using ..Thermodynamics
using LinearAlgebra

"""
    atmos_les_default_init(bl, currtime)

Initialize the 'AtmosLESDefault' diagnostics group.
"""
function atmos_les_default_init(dgngrp::DiagnosticsGroup, currtime)
    atmos_collect_onetime(Settings.mpicomm, Settings.dg, Settings.Q)

    return nothing
end

# Simple horizontal averages
function vars_atmos_les_default_simple(m::AtmosModel, FT)
    @vars begin
        u::FT
        v::FT
        w::FT
        avg_rho::FT             # ρ
        rho::FT                 # ρρ
        temp::FT
        pres::FT
        thd::FT                 # θ_dry
        et::FT                  # e_tot
        ei::FT                  # e_int
        ht::FT
        hi::FT
        w_ht_sgs::FT

        moisture::vars_atmos_les_default_simple(m.moisture, FT)
    end
end
vars_atmos_les_default_simple(::MoistureModel, FT) = @vars()
function vars_atmos_les_default_simple(m::EquilMoist, FT)
    @vars begin
        qt::FT                  # q_tot
        ql::FT                  # q_liq
        qv::FT                  # q_vap
        thv::FT                 # θ_vir
        thl::FT                 # θ_liq
        w_qt_sgs::FT
    end
end
num_atmos_les_default_simple_vars(m, FT) =
    varsize(vars_atmos_les_default_simple(m, FT))
atmos_les_default_simple_vars(m, array) =
    Vars{vars_atmos_les_default_simple(m, eltype(array))}(array)

function atmos_les_default_simple_sums!(
    atmos::AtmosModel,
    state_conservative,
    state_gradient_flux,
    state_auxiliary,
    thermo,
    currtime,
    MH,
    sums,
)
    sums.u += MH * state_conservative.ρu[1]
    sums.v += MH * state_conservative.ρu[2]
    sums.w += MH * state_conservative.ρu[3]
    sums.avg_rho += MH * state_conservative.ρ
    sums.rho += MH * state_conservative.ρ * state_conservative.ρ
    sums.temp += MH * thermo.temp * state_conservative.ρ
    sums.pres += MH * thermo.pres * state_conservative.ρ
    sums.thd += MH * thermo.θ_dry * state_conservative.ρ
    sums.et += MH * state_conservative.ρe
    sums.ei += MH * thermo.e_int * state_conservative.ρ
    sums.ht += MH * thermo.h_tot * state_conservative.ρ
    sums.hi += MH * thermo.h_int * state_conservative.ρ

    ν, D_t, _ = turbulence_tensors(
        atmos,
        state_conservative,
        state_gradient_flux,
        state_auxiliary,
        currtime,
    )
    d_h_tot = -D_t .* state_gradient_flux.∇h_tot
    sums.w_ht_sgs += MH * d_h_tot[end] * state_conservative.ρ

    atmos_les_default_simple_sums!(
        atmos.moisture,
        state_conservative,
        state_gradient_flux,
        thermo,
        MH,
        D_t,
        sums,
    )

    return nothing
end
function atmos_les_default_simple_sums!(
    ::MoistureModel,
    state_conservative,
    state_gradient_flux,
    thermo,
    MH,
    D_t,
    sums,
)
    return nothing
end
function atmos_les_default_simple_sums!(
    moist::EquilMoist,
    state_conservative,
    state_gradient_flux,
    thermo,
    MH,
    D_t,
    sums,
)
    sums.moisture.qt += MH * state_conservative.moisture.ρq_tot
    sums.moisture.ql += MH * thermo.moisture.q_liq * state_conservative.ρ
    sums.moisture.qv += MH * thermo.moisture.q_vap * state_conservative.ρ
    sums.moisture.thv += MH * thermo.moisture.θ_vir * state_conservative.ρ
    sums.moisture.thl += MH * thermo.moisture.θ_liq_ice * state_conservative.ρ
    d_q_tot = (-D_t) .* state_gradient_flux.moisture.∇q_tot
    sums.moisture.w_qt_sgs += MH * d_q_tot[end] * state_conservative.ρ

    return nothing
end

# Variances and covariances
function vars_atmos_les_default_ho(m::AtmosModel, FT)
    @vars begin
        var_u::FT               # u′u′
        var_v::FT               # v′v′
        var_w::FT               # w′w′
        w3::FT                  # w′w′w′
        tke::FT
        var_ei::FT              # e_int′e_int′

        cov_w_u::FT             # w′u′
        cov_w_v::FT             # w′v′
        cov_w_rho::FT           # w′ρ′
        cov_w_thd::FT           # w′θ_dry′
        cov_w_ei::FT            # w′e_int′

        moisture::vars_atmos_les_default_ho(m.moisture, FT)
    end
end
vars_atmos_les_default_ho(::MoistureModel, FT) = @vars()
function vars_atmos_les_default_ho(m::EquilMoist, FT)
    @vars begin
        var_qt::FT              # q_tot′q_tot′
        var_thl::FT             # θ_liq_ice′θ_liq_ice′

        cov_w_qt::FT            # w′q_tot′
        cov_w_ql::FT            # w′q_liq′
        cov_w_qv::FT            # w′q_vap′
        cov_w_thv::FT           # w′θ_v′
        cov_w_thl::FT           # w′θ_liq_ice′
        cov_qt_thl::FT          # q_tot′θ_liq_ice′
        cov_qt_ei::FT           # q_tot′e_int′
    end
end
num_atmos_les_default_ho_vars(m, FT) = varsize(vars_atmos_les_default_ho(m, FT))
atmos_les_default_ho_vars(m, array) =
    Vars{vars_atmos_les_default_ho(m, eltype(array))}(array)

function atmos_les_default_ho_sums!(
    atmos::AtmosModel,
    state_conservative,
    thermo,
    MH,
    ha,
    sums,
)
    u = state_conservative.ρu[1] / state_conservative.ρ
    u′ = u - ha.u
    v = state_conservative.ρu[2] / state_conservative.ρ
    v′ = v - ha.v
    w = state_conservative.ρu[3] / state_conservative.ρ
    w′ = w - ha.w
    e_int′ = thermo.e_int - ha.ei
    θ_dry′ = thermo.θ_dry - ha.thd

    sums.var_u += MH * u′^2 * state_conservative.ρ
    sums.var_v += MH * v′^2 * state_conservative.ρ
    sums.var_w += MH * w′^2 * state_conservative.ρ
    sums.w3 += MH * w′^3 * state_conservative.ρ
    sums.tke +=
        0.5 * (
            MH * u′^2 * state_conservative.ρ +
            MH * v′^2 * state_conservative.ρ +
            MH * w′^2 * state_conservative.ρ
        )
    sums.var_ei += MH * e_int′^2 * state_conservative.ρ

    sums.cov_w_u += MH * w′ * u′ * state_conservative.ρ
    sums.cov_w_v += MH * w′ * v′ * state_conservative.ρ
    sums.cov_w_rho +=
        MH * w′ * (state_conservative.ρ - ha.avg_rho) * state_conservative.ρ
    sums.cov_w_thd += MH * w′ * θ_dry′ * state_conservative.ρ
    sums.cov_w_ei += MH * w′ * e_int′ * state_conservative.ρ

    atmos_les_default_ho_sums!(
        atmos.moisture,
        state_conservative,
        thermo,
        MH,
        ha,
        w′,
        e_int′,
        sums,
    )

    return nothing
end
function atmos_les_default_ho_sums!(
    ::MoistureModel,
    state_conservative,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    return nothing
end
function atmos_les_default_ho_sums!(
    moist::EquilMoist,
    state_conservative,
    thermo,
    MH,
    ha,
    w′,
    e_int′,
    sums,
)
    q_tot = state_conservative.moisture.ρq_tot / state_conservative.ρ
    q_tot′ = q_tot - ha.moisture.qt
    q_liq′ = thermo.moisture.q_liq - ha.moisture.ql
    q_vap′ = thermo.moisture.q_vap - ha.moisture.qv
    θ_vir′ = thermo.moisture.θ_vir - ha.moisture.thv
    θ_liq_ice′ = thermo.moisture.θ_liq_ice - ha.moisture.thl

    sums.moisture.var_qt += MH * q_tot′^2 * state_conservative.ρ
    sums.moisture.var_thl += MH * θ_liq_ice′^2 * state_conservative.ρ

    sums.moisture.cov_w_qt += MH * w′ * q_tot′ * state_conservative.ρ
    sums.moisture.cov_w_ql += MH * w′ * q_liq′ * state_conservative.ρ
    sums.moisture.cov_w_qv += MH * w′ * q_vap′ * state_conservative.ρ
    sums.moisture.cov_w_thv += MH * w′ * θ_vir′ * state_conservative.ρ
    sums.moisture.cov_w_thl += MH * w′ * θ_liq_ice′ * state_conservative.ρ
    sums.moisture.cov_qt_thl += MH * q_tot′ * θ_liq_ice′ * state_conservative.ρ
    sums.moisture.cov_qt_ei += MH * q_tot′ * e_int′ * state_conservative.ρ

    return nothing
end

"""
    atmos_les_default_collect(bl, currtime)

Collect the various 'AtmosLESDefault' diagnostic variables for the
current timestep and write them into a file.
"""
function atmos_les_default_collect(dgngrp::DiagnosticsGroup, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    bl = dg.balance_law
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get needed arrays onto the CPU
    if array_device(Q) isa CPU
        host_state_conservative = Q.realdata
        host_state_auxiliary = dg.state_auxiliary.realdata
        host_vgeo = grid.vgeo
        host_state_gradient_flux = dg.state_gradient_flux.realdata
    else
        host_state_conservative = Array(Q.realdata)
        host_state_auxiliary = Array(dg.state_auxiliary.realdata)
        host_vgeo = Array(grid.vgeo)
        host_state_gradient_flux = Array(dg.state_gradient_flux.realdata)
    end
    FT = eltype(host_state_conservative)

    zvals = AtmosCollected.zvals
    repdvsr = AtmosCollected.repdvsr

    # Visit each node of the state variables array and:
    # - generate and store the thermo variables,
    # - accumulate the simple horizontal sums, and
    # - determine the cloud fraction, top and base
    #
    thermo_array =
        [zeros(FT, num_thermo(bl, FT)) for _ in 1:npoints, _ in 1:nrealelem]
    simple_sums = [
        zeros(FT, num_atmos_les_default_simple_vars(bl, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    ql_gt_0_z = [zeros(FT, (Nq * Nq * nhorzelem)) for _ in 1:(Nqk * nvertelem)]
    ql_gt_0_full = zeros(FT, (Nq * Nq * nhorzelem))
    # In honor of PyCLES!
    cld_top = FT(-100000)
    cld_base = FT(100000)
    @visitQ nhorzelem nvertelem Nqk Nq begin
        evk = Nqk * (ev - 1) + k

        state_conservative =
            extract_state_conservative(dg, host_state_conservative, ijk, e)
        state_gradient_flux =
            extract_state_gradient_flux(dg, host_state_gradient_flux, ijk, e)
        state_auxiliary =
            extract_state_auxiliary(dg, host_state_auxiliary, ijk, e)
        MH = host_vgeo[ijk, grid.MHid, e]

        thermo = thermo_vars(bl, thermo_array[ijk, e])
        compute_thermo!(bl, state_conservative, state_auxiliary, thermo)

        simple = atmos_les_default_simple_vars(bl, simple_sums[evk])
        atmos_les_default_simple_sums!(
            bl,
            state_conservative,
            state_gradient_flux,
            state_auxiliary,
            thermo,
            currtime,
            MH,
            simple,
        )

        # FIXME properly
        if isa(bl.moisture, EquilMoist)
            if !iszero(thermo.moisture.q_liq)
                idx = (Nq * Nq * (eh - 1)) + (Nq * (j - 1)) + i
                ql_gt_0_z[evk][idx] = one(FT)
                ql_gt_0_full[idx] = one(FT)

                z = zvals[evk]
                cld_top = max(cld_top, z)
                cld_base = min(cld_base, z)
            end
        end
    end

    # reduce horizontal sums and cloud data across ranks and compute averages
    simple_avgs = [
        zeros(FT, num_atmos_les_default_simple_vars(bl, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    cld_frac = zeros(FT, 1, Nqk * nvertelem)
    for evk in 1:(Nqk * nvertelem)
        MPI.Allreduce!(simple_sums[evk], simple_avgs[evk], +, mpicomm)
        simple_avgs[evk] .= simple_avgs[evk] ./ repdvsr[evk]

        # FIXME properly
        if isa(bl.moisture, EquilMoist)
            tot_ql_gt_0_z = MPI.Reduce(sum(ql_gt_0_z[evk]), +, 0, mpicomm)
            tot_horz_z = MPI.Reduce(length(ql_gt_0_z[evk]), +, 0, mpicomm)
            if mpirank == 0
                cld_frac[evk] = tot_ql_gt_0_z / tot_horz_z
            end
        end
    end
    # FIXME properly
    if isa(bl.moisture, EquilMoist)
        cld_top = MPI.Reduce(cld_top, MPI.MAX, 0, mpicomm)
        if cld_top == FT(-100000)
            cld_top = NaN
        end
        cld_base = MPI.Reduce(cld_base, MPI.MIN, 0, mpicomm)
        if cld_base == FT(100000)
            cld_base = NaN
        end
        tot_ql_gt_0_full = MPI.Reduce(sum(ql_gt_0_full), +, 0, mpicomm)
        tot_horz_full = MPI.Reduce(length(ql_gt_0_full), +, 0, mpicomm)
        cld_cover = zero(FT)
        if mpirank == 0
            cld_cover = tot_ql_gt_0_full / tot_horz_full
        end
    end

    # complete density averaging
    simple_varnames = map(
        s -> startswith(s, "moisture.") ? s[10:end] : s,
        flattenednames(vars_atmos_les_default_simple(bl, FT)),
    )
    for vari in 1:length(simple_varnames)
        for evk in 1:(Nqk * nvertelem)
            simple_ha = atmos_les_default_simple_vars(bl, simple_avgs[evk])
            avg_rho = simple_ha.avg_rho
            if simple_varnames[vari] != "avg_rho"
                simple_avgs[evk][vari] /= avg_rho
            end
        end
    end

    # compute the variances and covariances
    ho_sums = [
        zeros(FT, num_atmos_les_default_ho_vars(bl, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    @visitQ nhorzelem nvertelem Nqk Nq begin
        evk = Nqk * (ev - 1) + k

        state_conservative =
            extract_state_conservative(dg, host_state_conservative, ijk, e)
        thermo = thermo_vars(bl, thermo_array[ijk, e])
        MH = host_vgeo[ijk, grid.MHid, e]

        simple_ha = atmos_les_default_simple_vars(bl, simple_avgs[evk])
        ho = atmos_les_default_ho_vars(bl, ho_sums[evk])
        atmos_les_default_ho_sums!(
            bl,
            state_conservative,
            thermo,
            MH,
            simple_ha,
            ho,
        )
    end

    # reduce across ranks and compute averages
    ho_avgs = [
        zeros(FT, num_atmos_les_default_ho_vars(bl, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    for evk in 1:(Nqk * nvertelem)
        MPI.Reduce!(ho_sums[evk], ho_avgs[evk], +, 0, mpicomm)
        if mpirank == 0
            ho_avgs[evk] .= ho_avgs[evk] ./ repdvsr[evk]
        end
    end

    # complete density averaging and prepare output
    if mpirank == 0
        varvals = OrderedDict()
        for vari in 1:length(simple_varnames)
            davg = zeros(FT, 1, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                davg[evk] = simple_avgs[evk][vari]
            end
            varvals[simple_varnames[vari]] = (("time", "z"), davg, Dict())
        end

        ho_varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_les_default_ho(bl, FT)),
        )
        for vari in 1:length(ho_varnames)
            davg = zeros(FT, 1, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                simple_ha = atmos_les_default_simple_vars(bl, simple_avgs[evk])
                avg_rho = simple_ha.avg_rho
                davg[evk] = ho_avgs[evk][vari] / avg_rho
            end
            varvals[ho_varnames[vari]] = (("time", "z"), davg, Dict())
        end

        if isa(bl.moisture, EquilMoist)
            varvals["cld_frac"] = (("time", "z"), cld_frac, Dict())
            varvals["cld_top"] = (("time",), cld_top, Dict())
            varvals["cld_base"] = (("time",), cld_base, Dict())
            varvals["cld_cover"] = (("time",), cld_cover, Dict())
        end

        # write output
        dprefix = @sprintf(
            "%s_%s_%s_num%04d",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
            dgngrp.num
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        write_data(
            dgngrp.writer,
            dfilename,
            OrderedDict("z" => (zvals, Dict())),
            varvals,
            currtime,
        )
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_les_default_fini(dgngrp::DiagnosticsGroup, currtime) end

ext: jl

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# This file was generated, do not modify it. # hide
index!(G; parallel_block=512)
nothing # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
include("lobpcg_hyper_impl.jl")

function lobpcg_hyper(A, X0; maxiter=100, prec=nothing,
                      tol=20size(A, 2)*eps(real(eltype(A))),
                      largest=false, n_conv_check=nothing, kwargs...)
    prec === nothing && (prec = I)

    @assert !largest "Only seeking the smallest eigenpairs is implemented."
    result = LOBPCG(A, X0, I, prec, tol, maxiter; n_conv_check=n_conv_check, kwargs...)

    n_conv_check === nothing && (n_conv_check = size(X0, 2))
    converged = maximum(result.residual_norms[1:n_conv_check]) < tol
    iterations = size(result.residual_history, 2) - 1

    merge(result, (iterations=iterations, converged=converged))
end


text:
"""
    T2mapOptions(; kwargs...)
    T2mapOptions(image::Array{T,4}; kwargs...) where {T}

Options structure for [`T2mapSEcorr`](@ref).
This struct collects keyword arguments passed to `T2mapSEcorr`, performs checks on parameter types and values, and assigns default values to unspecified parameters.

# Arguments

    $TYPEDFIELDS

!!! note
    The 5D array that is saved when `SaveNNLSBasis` is set to `true` has dimensions `MatrixSize x nTE x nT2`, and therefore is typically extremely large.
    If the flip angle is fixed via `SetFlipAngle`, however, this is not an issue as only the unique `nTE x nT2` 2D basis matrix is saved.

See also:
* [`T2mapSEcorr`](@ref)
"""
@with_kw_noshow struct T2mapOptions{T<:Real}
    "Perform T2-mapping using legacy algorithms."
    legacy::Bool = false

    "Size of first 3 dimensions of input 4D image. This argument has no default, but is inferred automatically as `size(image)[1:3]` when calling `T2mapSEcorr(image; kwargs...)`."
    MatrixSize::NTuple{3,Int}
    @assert all(MatrixSize .>= 1)

    "Number of echoes in input signal. This argument is has no default, but is inferred automatically as `size(image, 4)` when calling `T2mapSEcorr(image; kwargs...)`."
    nTE::Int
    @assert nTE >= 4

    "Interecho spacing (Units: seconds). This argument has no default."
    TE::T # seconds
    @assert TE > 0.0

    "Number of T2 times to estimate in the multi-exponential analysis. This argument has no default."
    nT2::Int
    @assert nT2 >= 2

    "Tuple of min and max T2 values (Units: seconds). This argument has no default."
    T2Range::NTuple{2,T} # seconds
    @assert 0.0 < T2Range[1] < T2Range[2]

    # "Variable interecho spacing"
    # vTEparam::Union{Tuple{T,T,Int}, Nothing} = nothing
    # @assert vTEparam === nothing || begin
    #     error("Variable interecho spacing is not yet supported")
    #     TE1, TE2, nTE1 = vTEparam
    #     0.0 < TE1 < TE2 && nTE1 < nTE && round(Int, TE2/TE1) ≈ TE2/TE1
    # end

    "Assumed value of T1 (Units: seconds)."
    T1::T = 1.0 # seconds
    @assert T1 > 0.0

    "First echo intensity cutoff for empty voxels."
    Threshold::T = !legacy ? 0.0 : 200.0
    @assert Threshold >= 0.0

    "Minimum refocusing angle for flip angle optimization (Units: degrees)."
    MinRefAngle::T = 50.0 # degrees
    @assert 0.0 <= MinRefAngle <= 180.0

    "During flip angle optimization, goodness of fit is checked for up to `nRefAngles` angles in the range `[MinRefAngle, 180]`. The optimal angle is then determined through interpolation from these samples."
    nRefAngles::Int = !legacy ? 32 : 8
    @assert nRefAngles >= 2

    "Initial number of angles to check during flip angle optimization before refinement near likely optima. Setting `nRefAnglesMin` equal to `nRefAngles` forces all angles to be checked."
    nRefAnglesMin::Int = !legacy ? min(5, nRefAngles) : nRefAngles
    @assert 2 <= nRefAnglesMin <= nRefAngles

    "Regularization routine to use. One of \"none\", \"chi2\", \"gcv\", or \"lcurve\", representing no regularization, `Chi2Factor`-based Tikhonov regularization, Generalized Cross-Validation method, or L-Curve based regularization, respectively."
    Reg::String
    @assert Reg ∈ ("none", "chi2", "gcv", "lcurve")

    "Constraint on ``\\chi^2`` used for regularization when `Reg == \"chi2\"`."
    Chi2Factor::Union{T,Nothing} = nothing
    @assert (Reg == "chi2" && Chi2Factor !== nothing && Chi2Factor > 1.0) || Reg != "chi2"

    "Refocusing pulse control angle for stimulated echo correction (Units: degrees)."
    SetRefConAngle::Union{T,Nothing} = 180.0 # degrees
    @assert SetRefConAngle === nothing || 0.0 <= SetRefConAngle <= 180.0

    "Instead of optimizing flip angle, use `SetFlipAngle` for all voxels (Units: degrees)."
    SetFlipAngle::Union{T,Nothing} = nothing
    @assert SetFlipAngle === nothing || 0.0 <= SetFlipAngle <= 180.0

    "Boolean flag to include a 3D array of the ``\\ell^2``-norms of the residuals from the NNLS fits in the output maps dictionary."
    SaveResidualNorm::Bool = false

    "Boolean flag to include a 4D array of the time domain decay curves resulting from the NNLS fits in the output maps dictionary."
    SaveDecayCurve::Bool = false

    "Boolean flag to include 3D arrays of the regularization parameters ``\\mu`` and resulting ``\\chi^2``-factors in the output maps dictionary."
    SaveRegParam::Bool = false

    "Boolean flag to include a 5D (or 2D if `SetFlipAngle` is used) array of NNLS basis matrices in the output maps dictionary."
    SaveNNLSBasis::Bool = SetFlipAngle !== nothing

    "Suppress printing to the console."
    Silent::Bool = false
end
T2mapOptions(args...; kwargs...) = T2mapOptions{Float64}(args...; kwargs...)
T2mapOptions(image::Array{T,4}; kwargs...) where {T} = T2mapOptions{T}(; kwargs..., MatrixSize = size(image)[1:3], nTE = size(image)[4])

t2_times(o::T2mapOptions{T}) where {T} = logrange(o.T2Range..., o.nT2)
flip_angles(o::T2mapOptions{T}) where {T} = o.SetFlipAngle === nothing ? collect(range(o.MinRefAngle, T(180); length = o.nRefAngles)) : T[o.SetFlipAngle]
refcon_angles(o::T2mapOptions{T}) where {T} = o.SetRefConAngle === nothing ? collect(range(o.MinRefAngle, T(180); length = o.nRefAngles)) : T[o.SetRefConAngle]

function _show_string(o::T2mapOptions)
    io = IOBuffer()
    print(io, "T2-distribution analysis settings:")
    fields = fieldsof(typeof(o), Vector)
    fields = fields[sortperm(uppercase.(string.(fields)))] # sort alphabetically, ignoring case
    padlen = 1 + maximum(f -> length(string(f)), fields)
    for f in fields
        (f == :vTEparam) && continue # skip
        print(io, "\n$(lpad(rpad(f, padlen), padlen + 4)): $(getfield(o,f))")
    end
    return String(take!(io))
end
Base.show(io::IO, ::MIME"text/plain", o::T2mapOptions) = print(io, _show_string(o))

"""
    T2partOptions(; kwargs...)
    T2partOptions(t2dist::Array{T,4}; kwargs...) where {T}

Options structure for [`T2partSEcorr`](@ref).
This struct collects keyword arguments passed to `T2partSEcorr`, performs checks on parameter types and values, and assigns default values to unspecified parameters.

# Arguments

    $TYPEDFIELDS

See also:
* [`T2partSEcorr`](@ref)
"""
@with_kw_noshow struct T2partOptions{T<:Real}
    "Calculate T2-parts using legacy algorithms."
    legacy::Bool = false

    "Size of first 3 dimensions of input 4D T2 distribution. This argument is has no default, but is inferred automatically as `size(t2dist)[1:3]` when calling `T2partSEcorr(t2dist; kwargs...)`."
    MatrixSize::NTuple{3,Int}
    @assert all(MatrixSize .>= 1)

    "Number of T2 times to use. This argument has no default."
    nT2::Int
    @assert nT2 >= 2

    "Tuple of min and max T2 values (Units: seconds). This argument has no default."
    T2Range::NTuple{2,T} # seconds
    @assert 0.0 < T2Range[1] < T2Range[2]

    "Tuple of min and max T2 values of the short peak window (Units: seconds). This argument has no default."
    SPWin::NTuple{2,T} # seconds
    @assert SPWin[1] < SPWin[2]

    "Tuple of min and max T2 values of the middle peak window (Units: seconds). This argument has no default."
    MPWin::NTuple{2,T} # seconds
    @assert MPWin[1] < MPWin[2]

    "Apply sigmoidal weighting to the upper limit of the short peak window in order to smooth the hard small peak window cutoff time. `Sigmoid` is the delta-T2 parameter, which is the distance in seconds on either side of the `SPWin` upper limit where the sigmoid curve reaches 10% and 90% (Units: seconds)."
    Sigmoid::Union{T,Nothing} = nothing
    @assert Sigmoid === nothing || Sigmoid > 0

    "Suppress printing to the console."
    Silent::Bool = false
end
T2partOptions(args...; kwargs...) = T2partOptions{Float64}(args...; kwargs...)
T2partOptions(t2dist::Array{T,4}; kwargs...) where {T} = T2partOptions{T}(; kwargs..., MatrixSize = size(t2dist)[1:3], nT2 = size(t2dist)[4])

function _show_string(o::T2partOptions)
    io = IOBuffer()
    print(io, "T2-parts analysis settings:")
    fields = fieldsof(typeof(o), Vector)
    fields = fields[sortperm(uppercase.(string.(fields)))] # sort alphabetically
    padlen = 1 + maximum(f -> length(string(f)), fields)
    for f in fields
        print(io, "\n$(lpad(rpad(f, padlen), padlen + 4)): $(getfield(o,f))")
    end
    return String(take!(io))
end
Base.show(io::IO, ::MIME"text/plain", o::T2partOptions) = print(io, _show_string(o))


text:
module TestDimacs10
using Test
import GraphIdx.Io.Dimacs10:
  parse_dimacs10_idx, parse_dimacs10_header, parse_dimacs10_edges
import GraphIdx

const test_str = """% some comments
%  ... to be ignored
% integer node numbers start at one
% number of nodes, number of edges
4 3
2 4
2
4
"""

@testset "DIMACS10 reader                " begin
    idx = parse_dimacs10_idx(IOBuffer(test_str))
    n, m, head, tail = parse_dimacs10_edges(IOBuffer(test_str), Int)
    @test parse_dimacs10_header(IOBuffer(test_str)) == (4, 3)
    @test GraphIdx.num_nodes(idx) == 4
    @test GraphIdx.num_edges(idx) == 3
    @test n == 4
    @test m == 3
    @test head == [1, 1, 3]
    @test tail == [2, 4, 4]
    @test idx == GraphIdx.BiAdjacentIndex([0, 2, 3, 4, 6] .+ 1,
                                          [1, 3, 0, 3, 0, 2] .+ 1)
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
"""
    Stacked(bs)
    Stacked(bs, ranges)
    stack(bs::Bijector{0}...) # where `0` means 0-dim `Bijector`

A `Bijector` which stacks bijectors together which can then be applied to a vector
where `bs[i]::Bijector` is applied to `x[ranges[i]]::UnitRange{Int}`.

# Arguments
- `bs` can be either a `Tuple` or an `AbstractArray` of 0- and/or 1-dimensional bijectors
  - If `bs` is a `Tuple`, implementations are type-stable using generated functions
  - If `bs` is an `AbstractArray`, implementations are _not_ type-stable and use iterative methods
- `ranges` needs to be an iterable consisting of `UnitRange{Int}`
  - `length(bs) == length(ranges)` needs to be true.

# Examples
```
b1 = Logit(0.0, 1.0)
b2 = Identity{0}()
b = stack(b1, b2)
b([0.0, 1.0]) == [b1(0.0), 1.0]  # => true
```
"""
struct Stacked{Bs, Rs} <: Bijector{1}
    bs::Bs
    ranges::Rs
end
Stacked(bs::Tuple) = Stacked(bs, ntuple(i -> i:i, length(bs)))
Stacked(bs::AbstractArray) = Stacked(bs, [i:i for i in 1:length(bs)])

# define nested numerical parameters
# TODO: replace with `Functors.@functor Stacked (bs,)` when
# https://github.com/FluxML/Functors.jl/pull/7 is merged
function Functors.functor(::Type{<:Stacked}, x)
    function reconstruct_stacked(xs)
        return Stacked(xs.bs, x.ranges)
    end
    return (bs = x.bs,), reconstruct_stacked
end

function Base.:(==)(b1::Stacked, b2::Stacked)
    bs1, bs2 = b1.bs, b2.bs
    if !(bs1 isa Tuple && bs2 isa Tuple || bs1 isa Vector && bs2 isa Vector)
        return false
    end
    return all(bs1 .== bs2) && all(b1.ranges .== b2.ranges)
end

isclosedform(b::Stacked) = all(isclosedform, b.bs)

stack(bs::Bijector{0}...) = Stacked(bs)

# For some reason `inv.(sb.bs)` was unstable... This works though.
inv(sb::Stacked) = Stacked(map(inv, sb.bs), sb.ranges)
# map is not type stable for many stacked bijectors as a large tuple
# hence the generated function
@generated function inv(sb::Stacked{A}) where {A <: Tuple}
    exprs = []
    for i = 1:length(A.parameters)
        push!(exprs, :(inv(sb.bs[$i])))
    end
    :(Stacked(($(exprs...), ), sb.ranges))
end

@generated function _transform(x, rs::NTuple{N, UnitRange{Int}}, bs::Bijector...) where N
    exprs = []
    for i = 1:N
        push!(exprs, :(bs[$i](x[rs[$i]])))
    end
    return :(vcat($(exprs...)))
end
function _transform(x, rs::NTuple{1, UnitRange{Int}}, b::Bijector)
    @assert rs[1] == 1:length(x)
    return b(x)
end
function (sb::Stacked{<:Tuple,<:Tuple})(x::AbstractVector{<:Real})
    y = _transform(x, sb.ranges, sb.bs...)
    @assert size(y) == size(x) "x is size $(size(x)) but y is $(size(y))"
    return y
end
# The Stacked{<:AbstractArray} version is not TrackedArray friendly
function (sb::Stacked)(x::AbstractVector{<:Real})
    N = length(sb.bs)
    N == 1 && return sb.bs[1](x[sb.ranges[1]])

    y = mapvcat(1:N) do i
        sb.bs[i](x[sb.ranges[i]])
    end
    @assert size(y) == size(x) "x is size $(size(x)) but y is $(size(y))"
    return y
end

(sb::Stacked)(x::AbstractMatrix{<:Real}) = eachcolmaphcat(sb, x)

function logabsdetjac(
    b::Stacked,
    x::AbstractVector{<:Real}
)
    N = length(b.bs)
    init = sum(logabsdetjac(b.bs[1], x[b.ranges[1]]))

    return if N > 1
        init + sum(2:N) do i
            sum(logabsdetjac(b.bs[i], x[b.ranges[i]]))
        end
    else
        init
    end
end

function logabsdetjac(
    b::Stacked{<:Tuple{Vararg{<:Any, N}}, <:Tuple{Vararg{<:Any, N}}},
    x::AbstractVector{<:Real}
) where {N}
    init = sum(logabsdetjac(b.bs[1], x[b.ranges[1]]))

    return if N == 1
        init
    else
        init + sum(2:N) do i
            sum(logabsdetjac(b.bs[i], x[b.ranges[i]]))
        end
    end
end

function logabsdetjac(b::Stacked, x::AbstractMatrix{<:Real})
    return map(eachcol(x)) do c
        logabsdetjac(b, c)
    end
end

# Generates something similar to:
#
# quote
#     (y_1, _logjac) = forward(b.bs[1], x[b.ranges[1]])
#     logjac = sum(_logjac)
#     (y_2, _logjac) = forward(b.bs[2], x[b.ranges[2]])
#     logjac += sum(_logjac)
#     return (rv = vcat(y_1, y_2), logabsdetjac = logjac)
# end
@generated function forward(b::Stacked{<:Tuple{Vararg{<:Any, N}}, <:Tuple{Vararg{<:Any, N}}}, x::AbstractVector) where {N}
    expr = Expr(:block)
    y_names = []

    push!(expr.args, :((y_1, _logjac) = forward(b.bs[1], x[b.ranges[1]])))
    # TODO: drop the `sum` when we have dimensionality
    push!(expr.args, :(logjac = sum(_logjac)))
    push!(y_names, :y_1)
    for i = 2:N
        y_name = Symbol("y_$i")
        push!(expr.args, :(($y_name, _logjac) = forward(b.bs[$i], x[b.ranges[$i]])))

        # TODO: drop the `sum` when we have dimensionality
        push!(expr.args, :(logjac += sum(_logjac)))

        push!(y_names, y_name)
    end

    push!(expr.args, :(return (rv = vcat($(y_names...)), logabsdetjac = logjac)))
    return expr
end

function forward(sb::Stacked, x::AbstractVector)
    N = length(sb.bs)
    yinit, linit = forward(sb.bs[1], x[sb.ranges[1]])
    logjac = sum(linit)
    ys = mapvcat(drop(sb.bs, 1), drop(sb.ranges, 1)) do b, r
        y, l = forward(b, x[r])
        logjac += sum(l)
        y
    end
    return (rv = vcat(yinit, ys), logabsdetjac = logjac)
end


text:
function prepare_demand_data(out_dir, files_ID, time_instances, week_day_list, month_w)
    ############
    #Read in the demand file
    for instance in time_instances
        for day in week_day_list
            file = open( out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * string(day) * "_" * instance * ".txt")
            demands = Dict()
            s = 0
            cnt = 0
            for line in eachline(file)
                if contains(line, "Origin")
                    s = parse(Int, split(line)[2])

                elseif contains(line, " : ")
                    pairs_ = split(line, ";")

                    for pair in pairs_
                        if contains(pair, ":")
                            pair_vals = split(pair, ":")
                            t = parse(Int, pair_vals[1])
                            demand = parse(Float64, pair_vals[2])
                            demands[(s,t)] = demand
                        end
                    end
                end
            end

            close(file)

            outfile = open(out_dir * "/OD_demands/demands_" * month_w * "_" * string(day) * "_" * instance * ".json", "w")

            JSON.print(outfile, demands)

            close(outfile)
        end

    end
end

using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np
@pyimport json
@pyimport os
@pyimport pickle
@pyimport pandas as pd
@pyimport collections

@pyimport parameters_julia
@pyimport utils_julia

out_dir = parameters_julia.out_dir;
files_ID = parameters_julia.files_ID;
month_w = parameters_julia.month_w;
year = parameters_julia.year;
time_instances = parameters_julia.instances_ID;
deg_grid = parameters_julia.deg_grid;
c_grid = parameters_julia.c_grid;
lamb_grid = parameters_julia.lamb_grid;
week_day_list = parameters_julia.week_day_list;

prepare_demand_data(out_dir, files_ID, time_instances, week_day_list, month_w)


text:
using Test
using MPI

MPI.Init()

info = MPI.Info()
@test typeof(info) == MPI.Info
@test info.val == MPI.INFO_NULL.val

testinfo(;kwargs...) = MPI.Info(kwargs...)

info = testinfo(foo="fast", bar=true, baz=[10, -2])
@test info.val != MPI.INFO_NULL.val

@test length(info) == 3
@test info[:foo] == "fast"

@test sort!(collect(keys(info))) == [:bar, :baz, :foo]
@test sort!(collect(info)) == [:bar=>"true", :baz=>"10,-2", :foo=>"fast"]

delete!(info, :bar)
@test_throws KeyError info[:bar]

MPI.free(info)
@test info.val == MPI.INFO_NULL.val

MPI.Finalize()
@test MPI.Finalized()


text:
n = 10 # rand(1:10)
A = matrixdepot("invhilb", n)
@test issymmetric(A)
@test isposdef(A)
println("'invhilb' passed test...")


text:
using Documenter, CharibdeOptim

makedocs(
    modules = [CharibdeOptim],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CharibdeOptim.jl",
    authors  = "Yashvardhan Sharma",
    pages = [
        "Home" => "index.md"
         ]
)

deploydocs(
    repo = "github.com/yashcodes/CharibdeOptim.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
# ---
# title: "`mlaplacian`"
# id: op_mlaplacian
# cover: assets/mlaplacian.png
# ---

# Laplacian operator is defined as the difference between external gradient and internal
# gradient -- `mgradient(img; mode=:external) - mgradient(img; mode=:internal)`.

using ImageMorphology
using TestImages
using ImageBase
using ImageShow

img = restrict(testimage_dip3e("fig0940")) # rice
img01 = @. Gray(img > 0.5)
mosaic(img, img01; nrow=1)

#-

mosaic(mlaplacian(img), mlaplacian(img01); nrow=1)

# Note that laplacian can produce negative values:

A = falses(7, 7)
A[3:5, 3:5] .= true
A[4, 4] = false
Int.(mlaplacian(A))

# ## See also

# ImageBase.jl provides the finite-difference version of laplace operator `flaplacian`.
# ImageFiltering.jl provides `Laplacian` and `LaplacianOfGaussian` kernels.

fout = abs.(ImageBase.FiniteDiff.flaplacian(img))
mout = abs.(mlaplacian(img))
mosaic(fout, mout; nrow=1)

# For a comprehensive and more accurate documentation, please check the [`mlaplacian`](@ref)
# reference page.

# save cover image #src
using FileIO #src
mkpath("assets") #src
img = Gray.(restrict(testimage("blobs"))) #src
save("assets/mlaplacian.png", abs.(mlaplacian(img; r=2))) #src


text:
module MuttsVersionedTrees

using MuttsInterface
using AbstractTrees  # For printing

const insert!, delete! = Base.insert!, Base.delete!  # For export
export VTree, insert!, delete!

const mark_immutable! = MuttsInterface.mark_immutable!  # For export
export mark_immutable!  # From MuttsInterface

"""
    VTree{T}(value; left = nothing, right = nothing)

A versioned, mutable until shared binary branching tree.

```jldoctest
julia> begin
           head = VTree{Int}()
           # Insert elements
           head = insert!(head, 5);
           head = insert!(head, 1);
           head = insert!(head, 10);
           head = insert!(head, 2);
           # Print
           println(head)
           # Freeze `head` -- ready to share to other threads.
           mark_immutable!(head)
           println(head)
           # Keep branching from head
           head = insert!(head, 7);
           head = delete!(head, 1);
           println(head)
       end
VTreeNode{Int64}(5) ✔
├─ VTreeNode{Int64}(1) ✔
│  ├─ ∅
│  └─ VTreeNode{Int64}(2) ✔
└─ VTreeNode{Int64}(10) ✔

VTreeNode{Int64}(5) 🔒
├─ VTreeNode{Int64}(1) 🔒
│  ├─ ∅
│  └─ VTreeNode{Int64}(2) 🔒
└─ VTreeNode{Int64}(10) 🔒

VTreeNode{Int64}(5) ✔
├─ VTreeNode{Int64}(2) 🔒
└─ VTreeNode{Int64}(10) ✔
   ├─ VTreeNode{Int64}(7) ✔
   └─ ∅

```
"""
VTree

struct EmptyVTree{T} end
Base.convert(::Type{EmptyVTree{T}}, ::Nothing) where T = EmptyVTree{T}()
Base.show(io::IO, ::EmptyVTree) = print(io, "∅")

@mutt struct VTreeNode{T}
    value :: T
    left  :: Union{EmptyVTree, VTreeNode{T}}
    right :: Union{EmptyVTree, VTreeNode{T}}
    VTreeNode{T}(v, l,r) where T = new{T}(v, l,r)
    VTreeNode{T}(v; left = nothing, right = nothing) where T = new{T}(v, left, right)
end

const VTree{T} = Union{EmptyVTree, VTreeNode{T}}
Base.convert(::Type{VTree{T}}, ::Nothing) where T = EmptyVTree{T}()

# ----- Helper Constructor --------
VTree{T}() where T = EmptyVTree{T}()
VTree{T}(x, args...; kwargs...) where T = VTreeNode{T}(x, args...; kwargs...)

AbstractTrees.children(v::EmptyVTree{T}) where T = VTree{T}[]
function AbstractTrees.children(v::VTreeNode)
    if (v.left isa EmptyVTree && v.right isa EmptyVTree)
        ()
    else
        (v.left, v.right)
    end
end
function AbstractTrees.printnode(io::IO, v::VTreeNode{T}) where T
    mutable_emoji = is_mutts_mutable(v) ? "✔" : "🔒"
    print(io,"VTreeNode{$T}($(v.value)) $mutable_emoji")
end

Base.show(io::IO, v::VTree) = AbstractTrees.print_tree(io, v)

# ----- Copy needed for MuttsInterface.branch!() ---------------
MuttsInterface.make_mutable_copy(v::EmptyVTree{T}) where T = EmptyVTree{T}()
MuttsInterface.make_mutable_copy(v::VTreeNode{T}) where T = VTreeNode{T}(v.value, v.left, v.right)

# --- Insert! -----------------
"""
    head = insert!(head::VTree, value)

Return a new VTree with `value` inserted into the sorted tree `head`, mutating it if
possible, or returning a new branched copy if needed.

```jldoctest
julia> head1 = mark_immutable!(head)  # Store backup of head
VTreeNode{Int64}(5) 🔒
├─ VTreeNode{Int64}(1) 🔒
└─ VTreeNode{Int64}(10) 🔒
   ├─ VTreeNode{Int64}(7) 🔒
   └─ ∅

julia> head = insert!(head, 2)
VTreeNode{Int64}(5) ✔
├─ VTreeNode{Int64}(1) ✔
│  ├─ ∅
│  └─ VTreeNode{Int64}(2) ✔
└─ VTreeNode{Int64}(10) 🔒
   ├─ VTreeNode{Int64}(7) 🔒
   └─ ∅

julia> head1
VTreeNode{Int64}(5) 🔒
├─ VTreeNode{Int64}(1) 🔒
└─ VTreeNode{Int64}(10) 🔒
   ├─ VTreeNode{Int64}(7) 🔒
   └─ ∅
```
"""
Base.insert!(::EmptyVTree{T}, x) where T = VTreeNode{T}(x)
function Base.insert!(head::VTreeNode{T}, x) where T
    head = mutable_version(head)

    if x <= head.value
        if head.left isa EmptyVTree
            head.left = VTreeNode{T}(x)
        else
            head.left = insert!(head.left, x)
        end
    else
        if head.right isa EmptyVTree
            head.right = VTreeNode{T}(x)
        else
            head.right = insert!(head.right, x)
        end
    end
    return head
end

# --- delete! -----------------
"""
    head = delete!(head::VTree, x)

Return a new VTree with `value` deleted from the sorted tree `head`, mutating it if
possible, or returning a new branched copy if needed.

```jldoctest
julia> head = head2;  # Starting with a different example tree

julia> head1 = mark_immutable!(head)  # Store backup of head
VTreeNode{Int64}(5) 🔒
├─ VTreeNode{Int64}(1) 🔒
│  ├─ ∅
│  └─ VTreeNode{Int64}(2) 🔒
└─ VTreeNode{Int64}(10) 🔒
   ├─ VTreeNode{Int64}(7) 🔒
   │  ├─ ∅
   │  └─ VTreeNode{Int64}(8) 🔒
   └─ ∅

julia> head = delete!(head, 5)
VTreeNode{Int64}(7) ✔
├─ VTreeNode{Int64}(1) 🔒
│  ├─ ∅
│  └─ VTreeNode{Int64}(2) 🔒
└─ VTreeNode{Int64}(10) ✔
   ├─ VTreeNode{Int64}(8) 🔒
   └─ ∅

julia> head1
VTreeNode{Int64}(5) 🔒
├─ VTreeNode{Int64}(1) 🔒
│  ├─ ∅
│  └─ VTreeNode{Int64}(2) 🔒
└─ VTreeNode{Int64}(10) 🔒
   ├─ VTreeNode{Int64}(7) 🔒
   │  ├─ ∅
   │  └─ VTreeNode{Int64}(8) 🔒
   └─ ∅
```
"""
Base.delete!(e::EmptyVTree, _) = e
function Base.delete!(head::VTreeNode, x)
    head = mutable_version(head)

    if x == head.value
        if head.left isa EmptyVTree
            # Promote right child to parent
            head = head.right
        elseif head.right isa EmptyVTree
            head = head.left
        else
            #
            head.right, popped_val = _pop_leftmost!(head.right)
            head.value = popped_val
        end
    elseif x < head.value
        head.left = delete!(head.left, x)
    else
        head.right = delete!(head.right, x)
    end
    return head
end

# Used in delete!, above
_pop_leftmost!(head::EmptyVTree) = head
function _pop_leftmost!(head::VTreeNode)
    head = mutable_version(head)

    if head.left isa EmptyVTree
        popped_val = head.value
        head = delete!(head, popped_val)
        return (head, popped_val)
    else
        (head.left, popped_val) = _pop_leftmost!(head.left)
        return (head, popped_val)
    end
end

end  # module


text:
# vectors-tests.jl -
#
# Tests for vectorized operations.
#

if ! isdefined(:OptimPackNextGen)
    include("../src/OptimPackNextGen.jl")
end

module VectorsTests

using Base.Test
using OptimPackNextGen.Algebra

@testset "vfill" begin
    for T in (Float16, Float32, Float64)
        dims = (3,4,5)
        a = randn(T, dims)
        @test (0,0) == extrema(vfill!(a, 0) .- zeros(T,dims))
        @test (0,0) == extrema(vfill!(a, 1) .- ones(T,dims))
        @test (0,0) == extrema(vfill!(a,pi) .- fill!(similar(a), pi))
    end
end

@testset "vscale" begin
    for T in (Float16, Float32, Float64)
        dims = (3,4,5)
        a = randn(T, dims)
        b = vcreate(a)
        @test (0,0) == extrema(vscale!(b, 0,a) .- 0*a)
        @test (0,0) == extrema(vscale!(b, 1,a) .- a)
        @test (0,0) == extrema(vscale!(b,-1,a) .+ a)
        @test (0,0) == extrema(vscale!(b,pi,a) .- T(pi)*a)
    end
end

@testset "vupdate" begin
    for T in (Float16, Float32, Float64)
        dims = (3,4,5)
        a = randn(T, dims)
        d = randn(T, dims)
        b = vcreate(a)
        @test (0,0) == extrema(vupdate!(vcopy!(b,a), 0,d) .- a)
        @test (0,0) == extrema(vupdate!(vcopy!(b,a), 1,d) .- (a .+ d))
        @test (0,0) == extrema(vupdate!(vcopy!(b,a),-1,d) .- (a .- d))
        @test (0,0) == extrema(vupdate!(vcopy!(b,a),pi,d) .- (a .+ T(pi)*d))
    end
end

@testset "vproduct" begin
    for T in (Float16, Float32, Float64)
        dims = (3,4,5)
        a = randn(T, dims)
        b = randn(T, dims)
        c = vcreate(a)
        @test (0,0) == extrema(vproduct!(c,a,b) .- (a .* b))
    end
end

@testset "vcombine" begin
    for T in (Float16, Float32, Float64)
        dims = (3,4,5)
        a = randn(T, dims)
        b = randn(T, dims)
        d = vcreate(a)
        α = randn(T)
        β = randn(T)
        @test (0,0) == extrema(vcombine!(d, 0,a, 0,b) .- zeros(T,dims))
        @test (0,0) == extrema(vcombine!(d, 0,a, 1,b) .- b)
        @test (0,0) == extrema(vcombine!(d, 0,a,-1,b) .+ b)
        @test (0,0) == extrema(vcombine!(d, 0,a,pi,b) .- T(pi)*b)
        @test (0,0) == extrema(vcombine!(d, 1,a, 0,b) .- a)
        @test (0,0) == extrema(vcombine!(d,-1,a, 0,b) .+ a)
        @test (0,0) == extrema(vcombine!(d,pi,a, 0,b) .- T(pi)*a)
        @test (0,0) == extrema(vcombine!(d, 1,a, 1,b) .- (a .+ b))
        @test (0,0) == extrema(vcombine!(d,-1,a, 1,b) .- (b .- a))
        @test (0,0) == extrema(vcombine!(d, 1,a,-1,b) .- (a .- b))
        @test (0,0) == extrema(vcombine!(d,-1,a,-1,b) .+ (a .+ b))
        @test (0,0) == extrema(vcombine!(d, α,a, β,b) .- (T(α)*a + T(β)*b))
    end
end

end # module


text:
###################################################
## geodesy.jl
## fork of some code in geodesy.jl
###################################################

### Ellipsoid
# Specify datum for translation between LLA and other coordinate systems
struct Ellipsoid
    a::Float64        # Semi-major axis
    b::Float64        # Semi-minor axis
    e²::Float64       # Eccentricity squared
    e′²::Float64      # Second eccentricity squared
end

function Ellipsoid(; a::AbstractString="", b::AbstractString="", f_inv::AbstractString="")
    if isempty(a) || isempty(b) == isempty(f_inv)
        throw(ArgumentError("Specify parameter 'a' and either 'b' or 'f_inv'"))
    end
    if isempty(b)
        _ellipsoid_af(parse(BigFloat,a), parse(BigFloat,f_inv))
    else
        _ellipsoid_ab(parse(BigFloat,a), parse(BigFloat,b))
    end
end

function _ellipsoid_ab(a::BigFloat, b::BigFloat)
    e² = (a^2 - b^2) / a^2
    e′² = (a^2 - b^2) / b^2

    Ellipsoid(a, b, e², e′²)
end
function _ellipsoid_af(a::BigFloat, f_inv::BigFloat)
    b = a * (1 - inv(f_inv))

    _ellipsoid_ab(a, b)
end

### World Geodetic Coordinate System of 1984 (WGS 84)
# Standardized coordinate system for Earth
# Global ellipsoidal reference surface
const WGS84  = Ellipsoid(a = "6378137.0", f_inv = "298.257223563")

### Point in Latitude-Longitude-Altitude (LLA) coordinates
# Used to store node data in OpenStreetMap XML files
struct LLA
    lat::Float64
    lon::Float64
    alt::Float64
end
LLA(lat, lon) = LLA(lat, lon, 0.0)

### Point in Earth-Centered-Earth-Fixed (ECEF) coordinates
# Global cartesian coordinate system rotating with the Earth
struct ECEF
    x::Float64
    y::Float64
    z::Float64
end

### Point in East-North-Up (ENU) coordinates
# Local cartesian coordinate system
# Linearized about a reference point
struct ENU
    east::Float64
    north::Float64
    up::Float64
end
ENU(x, y) = ENU(x, y, 0.0)

# Given a reference point for linarization
function ENU(lla::LLA, lla_ref::LLA, datum::Ellipsoid = WGS84)
    ecef = ECEF(lla, datum)
    return ENU(ecef, lla_ref, datum)
end

###############################
### LLA to ECEF coordinates ###
###############################

function ECEF(lla::LLA, datum::Ellipsoid = WGS84)
    ϕdeg, λdeg, h = lla.lat, lla.lon, lla.alt
    d = datum

    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)
    sinλ, cosλ = sind(λdeg), cosd(λdeg)

    N = d.a / sqrt(1 - d.e² * sinϕ^2)  # Radius of curvature (meters)

    x = (N + h) * cosϕ * cosλ
    y = (N + h) * cosϕ * sinλ
    z = (N * (1 - d.e²) + h) * sinϕ

    return ECEF(x, y, z)
end

function ENU(ecef::ECEF, lla_ref::LLA, datum::Ellipsoid = WGS84)
    ϕdeg, λdeg = lla_ref.lat, lla_ref.lon

    ecef_ref = ECEF(lla_ref, datum)
    ∂x = ecef.x - ecef_ref.x
    ∂y = ecef.y - ecef_ref.y
    ∂z = ecef.z - ecef_ref.z

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [     -sinλ       cosλ  0.0
    #      -cosλ*sinϕ -sinλ*sinϕ cosϕ
    #       cosλ*cosϕ  sinλ*cosϕ sinϕ]
    #
    # east, north, up = R * [∂x, ∂y, ∂z]
    east  = ∂x * -sinλ      + ∂y * cosλ       + ∂z * 0.0
    north = ∂x * -cosλ*sinϕ + ∂y * -sinλ*sinϕ + ∂z * cosϕ
    up    = ∂x * cosλ*cosϕ  + ∂y * sinλ*cosϕ  + ∂z * sinϕ

    return ENU(east, north, up)
end

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
module XORShift

import Base: rand, rand!

type XORShiftStar1024 <: AbstractRNG
    x00::UInt64
    x01::UInt64
    x02::UInt64
    x03::UInt64
    x04::UInt64
    x05::UInt64
    x06::UInt64
    x07::UInt64
    x08::UInt64
    x09::UInt64
    x10::UInt64
    x11::UInt64
    x12::UInt64
    x13::UInt64
    x14::UInt64
    x15::UInt64
    p::UInt64
end
XORShiftStar() = XORShiftStar1024(rand(UInt64, 16)..., 0)

const globalXORShiftStar1024 = XORShiftStar()

fromUInt64(::Type{UInt64}, u::UInt64) = u
fromUInt64(::Type{Float64}, u::UInt64) = 5.421010862427522e-20 * u

function rand(rng::XORShiftStar1024, ::Type{UInt64})
    @inbounds begin
        p  = rng.p
        s0 = getfield(rng, reinterpret(Int, p)+1)
        p = (p + 1) % 16
        s1 = getfield(rng, reinterpret(Int, p)+1)
        s1 $= s1 << 31
        s1 $= s1 >> 11
        s0 $= s0 >> 30
        s0 $= s1
        unsafe_store!(convert(Ptr{UInt64}, pointer_from_objref(rng))+8,
                      s0, reinterpret(Int, p)+1)
        rng.p = p
    end
    return s0 * 1181783497276652981
end

rand{T<:Number}(rng::XORShiftStar1024, ::Type{T}) = fromUInt64(T, rand(rng, UInt64))

@eval function rand!{T}(rng::XORShiftStar1024, a::AbstractArray{T})
    p = rng.p
    n = length(a)
    q = min((16 - p) % 16, n)
    r = (n - q) % 16
    @inbounds for i = 1:q
        a[i] = rand(rng, T)
    end
    $(Expr(:block, [
        let x = symbol(@sprintf("x%02d", k % 16))
            :($x = rng.$x)
        end
    for k = 0:15 ]...))
    r = q + (n-q) >> 4 << 4
    @inbounds for i = q+1:16:r
        $(Expr(:block, [
            let x0 = symbol(@sprintf("x%02d", k % 16)),
                x1 = symbol(@sprintf("x%02d", (k + 1) % 16))
                quote
                    s0 = $x0
                    s1 = $x1
                    s1 $= s1 << 31
                    s1 $= s1 >> 11
                    s0 $= s0 >> 30
                    s0 $= s1
                    $x1 = s0
                    a[i + $k] = fromUInt64(T, s0 * 1181783497276652981)
                end
            end
        for k = 0:15 ]...))
    end
    $(Expr(:block, [
        let x = symbol(@sprintf("x%02d", k % 16))
            :(rng.$x = $x)
        end
    for k = 0:15 ]...))
    rng.p = 0
    @inbounds for i = r+1:n
        a[i] = rand(rng, T)
    end
    return a
end

end # module

max_stars_repo_name: JuliaLang/XORShiftRNG

max_issues_repo_name: JuliaLang/XORShiftRNG

max_forks_repo_name: JuliaLang/XORShiftRNG


text:
module State

const Iterate           = Union{Number, Vector, Nothing}
const FloatVoid         = Union{Number, Nothing}
const MatrixType        = Any #Union{Number, AbstractArray, Nothing}

include("GenericStatemod.jl")
include("LSAtTmod.jl")
include("NLPAtXmod.jl")

export AbstractState, GenericState, update!
export LSAtT, copy, update!
export NLPAtX, update! #, convert_nlp, convert_ls

end # module


text:
################################################################################
######  Structure holding information computed at the end of a sampling   ######
################################################################################
mutable struct SampledObservables{T} <: EvaluatedAlgorithm
    ObsNames::Vector{Symbol}
    ObsAve::Vector{T}
    ObsVals::Vector{Vector{T}}
end

function SampledObservables(oprob, net::NeuralNetwork)
    names = deepcopy(oprob.Names)
    type  = out_type(net)
    for el=oprob.ObservablesTransposed
        type = promote_type(type, eltype(el))
    end
    ObsAve  = zeros(type, length(oprob.ObservablesTransposed))
    ObsVals = [Vector{type}() for i=1:length(oprob.ObservablesTransposed)]

    return SampledObservables(names, ObsAve, ObsVals)
end

EvaluatedNetwork(alg::ObservablesProblem, net) =
    SampledObservables(alg, net)

function add!(acc::SampledObservables, o::SampledObservables)
    for i=1:length(acc.ObsNames)
        @assert acc.ObsNames[i] == o.ObsNames[i]
        acc.ObsAve[i] = acc.ObsAve[i] + o.ObsAve[i]
    end
    acc.ObsAve[i] ./= length(acc.ObsNames)
    acc
end

################################################################################
######      Cache holding the information generated along a sampling      ######
################################################################################
mutable struct MCMCObsEvaluationCache{T,T2,S} <: EvaluationSamplingCache
    ObsNames::Vector{Symbol}
    ObsAve::T
    ObsVals::T2
    Zave::Real

    # caches
    σ::S
end

function MCMCObsEvaluationCache(net::NeuralNetwork, obs_prob::ObservablesProblem)
    n_obs    = length(obs_prob.ObservablesTransposed)
    obs_ave  = Vector{ComplexF64}(undef, n_obs)
    obs_vals = [Vector{ComplexF64}() for i=1:n_obs]

    σ        = state(obs_prob, net)

    cache=MCMCObsEvaluationCache(obs_prob.Names, obs_ave, obs_vals, 0.0, σ)
    zero!(cache)
    cache
end
SamplingCache(alg::ObservablesProblem, prob::ObservablesProblem, net) = MCMCObsEvaluationCache(net, prob)
SamplingCache(prob::ObservablesProblem, net) = MCMCObsEvaluationCache(net, prob)

function zero!(comp_vals::MCMCObsEvaluationCache)
    comp_vals.ObsAve   .= 0.0
    comp_vals.Zave      = 0.0
    resize!.(comp_vals.ObsVals, 0)
end

function add!(acc::MCMCObsEvaluationCache, o::MCMCObsEvaluationCache)
   acc.ObsAve   .+= o.ObsAve
   acc.Zave      += o.Zave
   append!.(acc.ObsVals, o.ObsVals)
end

function evaluation_post_sampling(vals::MCMCObsEvaluationCache)
    obs_ave  = vals.ObsAve./vals.Zave
    obs_vals = [cumsum(ob_t)./(1:length(ob_t)) for
                                    ob_t=vals.ObsVals]

    SampledObservables(vals.ObsNames, obs_ave, obs_vals)
end

function evaluation_post_sampling!(out::SampledObservables,
                                   vals::MCMCObsEvaluationCache,
                                   sampler_steps = vals.Zave)
    out.ObsAve .= vals.ObsAve./vals.Zave

    for (obv_t, ob_t) = zip(out.ObsVals, vals.ObsVals)
        resize!(obv_t, length(ob_t))
        cumsum!(obv_t, ob_t)
        obv_t ./= 1:length(ob_t)
    end
    return out
end


text:
export ResidueField

################################################################################
#
#  Residue field construction for arbitrary prime ideals
#
################################################################################

# Assume that m is a surjective morphism pi: O -> A, where A is a simple algebra
# over the prime field Z/pZ.
# This functions returns - a in O, such that pi(a) is a primitive element
#                        - f in Z[X], such that f is the minimal polynomial of
#                          pi(a)
#                        - a vector of fmpz_mat B, such that
#                          pi(basis(O)[i]) = sum_j B[i][1, j] * pi(a)^j
function compute_residue_field_data(m)
  O = domain(m)
  basisO = basis(O, copy = false)
  B = codomain(m)
  primB, minprimB, getcoordpowerbasis = _as_field(B)
  f = degree(minprimB)
  prim_elem = m\(primB)
  min_poly_prim_elem = fmpz_poly(fmpz[lift(coeff(minprimB, i)) for i in 0:degree(minprimB)])
  min_poly_prim_elem.parent = FmpzPolyRing(FlintZZ, :$, false)
  basis_in_prim = Vector{fmpz_mat}(undef, degree(O))
  for i in 1:degree(O)
    basis_in_prim[i] = zero_matrix(FlintZZ, 1, f)
    t = getcoordpowerbasis(m(basisO[i]))
    for j in 1:f
      basis_in_prim[i][1, j] = lift(t[1, j])
    end
  end
  return prim_elem, min_poly_prim_elem, basis_in_prim
end

# Compute the residue field data and store it in the prime P given the map m
function compute_residue_field_data!(P, m)
  P.prim_elem, P.min_poly_prim_elem, P.basis_in_prim = compute_residue_field_data(m)
  return nothing
end

# Compute the residue field data given the prime P
function compute_residue_field_data!(P)
  p = minimum(P)
  if fits(Int, p)
    smallp = Int(p)
    A, m = AlgAss(order(P), P, smallp)
    compute_residue_field_data!(P, m)
  else
    AA, mm = AlgAss(order(P), P, p)
    compute_residue_field_data!(P, mm)
  end
  return nothing
end

# Get the residue field data. This is the function one should use, since the
# data is often cached.
function get_residue_field_data(P)
  if isdefined(P, :prim_elem)
    return P.prim_elem, P.min_poly_prim_elem, P.basis_in_prim
  else
    compute_residue_field_data!(P)
    get_residue_field_data(P)
  end
end

################################################################################
#
#  Residue field construction for nonindex divisors
#
################################################################################

# It is assumed that p is not an index divisor
function _residue_field_nonindex_divisor_helper(f::fmpq_poly, g::fmpq_poly, p, degree_one::Type{Val{S}} = Val{false}) where S
  R = GF(p, cached = false)

  Zy, y = PolynomialRing(FlintZZ, "y", cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)

  gmodp = Rx(g)
  fmodp = Rx(f)

  h = gcd(gmodp,fmodp)

	if degree_one === Val{true}
    return R, h
	else
    if isa(p, Int)
    	F3 = FqNmodFiniteField(h, :$, false)
      return F3, h
    elseif isa(p, fmpz)
    	F4 = FqFiniteField(h, :$, false)
      return F4, h
  	end
	end
end

function _residue_field_nonindex_divisor(O, P, small::Type{Val{T}} = Val{false}, degree_one::Type{Val{S}} = Val{false}) where {S, T}
  # This code assumes that P comes from prime_decomposition
  @assert has_2_elem(P) && isprime_known(P) && isprime(P)

  gtwo = P.gen_two

  f = nf(O).pol
  g = parent(f)(elem_in_nf(gtwo))

  if small === Val{true}
    @assert fits(Int, minimum(P, copy = false))
    F, h = _residue_field_nonindex_divisor_helper(f, g, Int(minimum(P)), degree_one)
    mF = Mor(O, F, h)
    mF.P = P
    return F, mF
  elseif small === Val{false}
    F, h = _residue_field_nonindex_divisor_helper(f, g, minimum(P), degree_one)
    mF = Mor(O, F, h)
    mF.P = P
    return F, mF
  end
end

################################################################################
#
#  Residue field construction for index divisors
#
################################################################################

function _residue_field_generic(O, P, small::Type{Val{T}} = Val{false}, degree_one::Type{Val{S}} = Val{false}) where {S, T}
  if small == Val{true}
    @assert fits(Int, minimum(P, copy = false))
    if degree_one === Val{true}
			f1 = NfOrdToGFMor(O, P)
      return codomain(f1), f1
		else
    	f = NfOrdToFqNmodMor(O, P)
    	return codomain(f), f
    end
  elseif small === Val{false}
    if degree_one === Val{true}
    	f3 = NfOrdToGFFmpzMor(O, P)
      return codomain(f3), f3
		else
			f2 = NfOrdToFqMor(O, P)
    	return codomain(f2), f2
		end
  end
end

################################################################################
#
#  High level functions
#
################################################################################
@doc Markdown.doc"""
    ResidueField(O::NfOrd, P::NfOrdIdl, check::Bool = true) -> Field, Map
Returns the residue field of the prime ideal $P$ together with th
projection map. If ```check``` is true, the ideal is checked for
being prime.
"""
function ResidueField(O::NfOrd, P::NfOrdIdl, check::Bool = true)
  if check
    !isprime(P) && error("Ideal must be prime")
  end
  if !ismaximal_known(O) || !ismaximal(O)
    return _residue_field_generic(O, P)
  end
  if !isindex_divisor(O, minimum(P)) && has_2_elem(P)
    return _residue_field_nonindex_divisor(O, P)
  else
    return _residue_field_generic(O, P)
  end
end

function ResidueFieldSmall(O::NfOrd, P::NfOrdIdl)
  p = minimum(P)
  !fits(Int, p) && error("Minimum of prime ideal must be small (< 64 bits)")
  if !ismaximal_known(O) || !ismaximal(O) || !isdefining_polynomial_nice(nf(O))
    return _residue_field_generic(O, P, Val{true})
  end
  if !isindex_divisor(O, minimum(P))
    return _residue_field_nonindex_divisor(O, P, Val{true})
  else
    return _residue_field_generic(O, P, Val{true})
  end
end

function ResidueFieldDegree1(O::NfOrd, P::NfOrdIdl)
  @assert degree(P) == 1
  if !ismaximal_known(O) || !ismaximal(O)
    return _residue_field_generic(O, P, Val{false}, Val{true})
  end
  if !isindex_divisor(O, minimum(P)) && has_2_elem(P)
    return _residue_field_nonindex_divisor(O, P, Val{false}, Val{true})
  else
    return _residue_field_generic(O, P, Val{false}, Val{true})
  end
end

function ResidueFieldSmallDegree1(O::NfOrd, P::NfOrdIdl)
  p = minimum(P)
  !fits(Int, p) && error("Minimum of prime ideal must be small (< 64 bits)")
  @assert degree(P) == 1
  if !ismaximal_known(O) || !ismaximal(O) || !isdefining_polynomial_nice(nf(O))
    return _residue_field_generic(O, P, Val{true}, Val{true})
  end
  if !isindex_divisor(O, minimum(P))
    return _residue_field_nonindex_divisor(O, P, Val{true}, Val{true})
  else
    return _residue_field_generic(O, P, Val{true}, Val{true})
  end
end