text:
"""
    BetaPrime(α,β)

The *Beta prime distribution* has probability density function

```math
f(x; \\alpha, \\beta) = \\frac{1}{B(\\alpha, \\beta)}
x^{\\alpha - 1} (1 + x)^{- (\\alpha + \\beta)}, \\quad x > 0
```

The Beta prime distribution is related to the [`Beta`](@ref) distribution via the
relation ship that if ``X \\sim \\operatorname{Beta}(\\alpha, \\beta)`` then ``\\frac{X}{1 - X}
\\sim \\operatorname{BetaPrime}(\\alpha, \\beta)``

```julia
BetaPrime()        # equivalent to BetaPrime(1, 1)
BetaPrime(a)       # equivalent to BetaPrime(a, a)
BetaPrime(a, b)    # Beta prime distribution with shape parameters a and b

params(d)          # Get the parameters, i.e. (a, b)
```

External links

* [Beta prime distribution on Wikipedia](http://en.wikipedia.org/wiki/Beta_prime_distribution)

"""
struct BetaPrime{T<:Real} <: ContinuousUnivariateDistribution
    α::T
    β::T
    BetaPrime{T}(α::T, β::T) where {T} = new{T}(α, β)
end

function BetaPrime(α::T, β::T; check_args=true) where {T<:Real}
    check_args && @check_args(BetaPrime, α > zero(α) && β > zero(β))
    return BetaPrime{T}(α, β)
end

BetaPrime(α::Real, β::Real) = BetaPrime(promote(α, β)...)
BetaPrime(α::Integer, β::Integer) = BetaPrime(float(α), float(β))
BetaPrime(α::Real) = BetaPrime(α, α)
BetaPrime() = BetaPrime(1.0, 1.0, check_args=false)

@distr_support BetaPrime 0.0 Inf

#### Conversions
function convert(::Type{BetaPrime{T}}, α::Real, β::Real) where T<:Real
    BetaPrime(T(α), T(β))
end
function convert(::Type{BetaPrime{T}}, d::BetaPrime{S}) where {T <: Real, S <: Real}
    BetaPrime(T(d.α), T(d.β), check_args=false)
end

#### Parameters

params(d::BetaPrime) = (d.α, d.β)
@inline partype(d::BetaPrime{T}) where {T<:Real} = T

#### Statistics

function mean(d::BetaPrime{T}) where T<:Real
    ((α, β) = params(d); β > 1 ? α / (β - 1) : T(NaN))
end

function mode(d::BetaPrime{T}) where T<:Real
    ((α, β) = params(d); α > 1 ? (α - 1) / (β + 1) : zero(T))
end

function var(d::BetaPrime{T}) where T<:Real
    (α, β) = params(d)
    β > 2 ? α * (α + β - 1) / ((β - 2) * (β - 1)^2) : T(NaN)
end

function skewness(d::BetaPrime{T}) where T<:Real
    (α, β) = params(d)
    if β > 3
        s = α + β - 1
        2(α + s) / (β - 3) * sqrt((β - 2) / (α * s))
    else
        return T(NaN)
    end
end

#### Evaluation

function logpdf(d::BetaPrime, x::Real)
    α, β = params(d)
    _x = max(0, x)
    z = xlogy(α - 1, _x) - (α + β) * log1p(_x) - logbeta(α, β)
    return x < 0 ? oftype(z, -Inf) : z
end

function zval(::BetaPrime, x::Real)
    y = max(x, 0)
    z = y / (1 + y)
    # map `Inf` to `Inf` (otherwise it returns `NaN`)
    return isinf(x) && x > 0 ? oftype(z, Inf) : z
end
xval(::BetaPrime, z::Real) = z / (1 - z)

for f in (:cdf, :ccdf, :logcdf, :logccdf)
    @eval $f(d::BetaPrime, x::Real) = $f(Beta(d.α, d.β; check_args=false), zval(d, x))
end

for f in (:quantile, :cquantile, :invlogcdf, :invlogccdf)
    @eval $f(d::BetaPrime, p::Real) = xval(d, $f(Beta(d.α, d.β; check_args=false), p))
end

#### Sampling

function rand(rng::AbstractRNG, d::BetaPrime)
    (α, β) = params(d)
    rand(rng, Gamma(α)) / rand(rng, Gamma(β))
end


text:
## Type aliases for convenience ##

typealias AbstractVector{T} AbstractArray{T,1}
typealias AbstractMatrix{T} AbstractArray{T,2}
typealias AbstractVecOrMat{T} Union(AbstractVector{T}, AbstractMatrix{T})

## Basic functions ##

size{T,n}(t::AbstractArray{T,n}, d) = (d>n ? 1 : size(t)[d])
eltype(x) = Any
eltype{T,n}(::AbstractArray{T,n}) = T
eltype{T,n}(::Type{AbstractArray{T,n}}) = T
eltype{T<:AbstractArray}(::Type{T}) = eltype(super(T))
iseltype(x,T) = eltype(x) <: T
isinteger(x::AbstractArray) = all(isinteger,x)
isinteger{T<:Integer,n}(x::AbstractArray{T,n}) = true
isreal(x::AbstractArray) = all(isreal,x)
isreal{T<:Real,n}(x::AbstractArray{T,n}) = true
ndims{T,n}(::AbstractArray{T,n}) = n
ndims{T,n}(::Type{AbstractArray{T,n}}) = n
ndims{T<:AbstractArray}(::Type{T}) = ndims(super(T))
nfilled(t::AbstractArray) = length(t)
length(t::AbstractArray) = prod(size(t))::Int
endof(a::AbstractArray) = length(a)
first(a::AbstractArray) = a[1]
first(a) = next(a,start(a))[1]
last(a) = a[end]

function stride(a::AbstractArray, i::Integer)
    if i > ndims(a)
        return length(a)
    end
    s = 1
    for n=1:(i-1)
        s *= size(a, n)
    end
    return s
end

strides(a::AbstractArray) = ntuple(ndims(a), i->stride(a,i))::Dims

function isassigned(a::AbstractArray, i::Int...)
    # TODO
    try
        a[i...]
        true
    catch
        false
    end
end

# used to compute "end" for last index
function trailingsize(A, n)
    s = 1
    for i=n:ndims(A)
        s *= size(A,i)
    end
    return s
end

## Bounds checking ##
checkbounds(sz::Int, i::Int) = 1 <= i <= sz || throw(BoundsError())
checkbounds(sz::Int, i::Real) = checkbounds(sz, to_index(i))
checkbounds(sz::Int, I::AbstractVector{Bool}) = length(I) == sz || throw(BoundsError())
checkbounds(sz::Int, r::Range{Int}) = isempty(r) || (minimum(r) >= 1 && maximum(r) <= sz) || throw(BoundsError())
checkbounds{T<:Real}(sz::Int, r::Range{T}) = checkbounds(sz, to_index(r))

function checkbounds{T <: Real}(sz::Int, I::AbstractArray{T})
    for i in I
        checkbounds(sz, i)
    end
end

checkbounds(A::AbstractArray, I::AbstractArray{Bool}) = size(A) == size(I) || throw(BoundsError())

checkbounds(A::AbstractArray, I) = checkbounds(length(A), I)

function checkbounds(A::AbstractMatrix, I::Union(Real,AbstractArray), J::Union(Real,AbstractArray))
    checkbounds(size(A,1), I)
    checkbounds(size(A,2), J)
end

function checkbounds(A::AbstractArray, I::Union(Real,AbstractArray), J::Union(Real,AbstractArray))
    checkbounds(size(A,1), I)
    checkbounds(trailingsize(A,2), J)
end

function checkbounds(A::AbstractArray, I::Union(Real,AbstractArray)...)
    n = length(I)
    if n > 0
        for dim = 1:(n-1)
            checkbounds(size(A,dim), I[dim])
        end
        checkbounds(trailingsize(A,n), I[n])
    end
end

## Bounds-checking without errors ##
in_bounds(l::Int, i::Integer) = 1 <= i <= l
function in_bounds(sz::Dims, I::Int...)
    n = length(I)
    for dim = 1:(n-1)
        1 <= I[dim] <= sz[dim] || return false
    end
    s = sz[n]
    for i = n+1:length(sz)
        s *= sz[i]
    end
    1 <= I[n] <= s
end

## Constructors ##

# default arguments to similar()
similar{T}(a::AbstractArray{T})               = similar(a, T, size(a))
similar   (a::AbstractArray, T)               = similar(a, T, size(a))
similar{T}(a::AbstractArray{T}, dims::Dims)   = similar(a, T, dims)
similar{T}(a::AbstractArray{T}, dims::Int...) = similar(a, T, dims)
similar   (a::AbstractArray, T, dims::Int...) = similar(a, T, dims)

function reshape(a::AbstractArray, dims::Dims)
    if prod(dims) != length(a)
        error("dimensions must be consistent with array size")
    end
    copy!(similar(a, dims), a)
end
reshape(a::AbstractArray, dims::Int...) = reshape(a, dims)

vec(a::AbstractArray) = reshape(a,length(a))
vec(a::AbstractVector) = a

function squeeze(A::AbstractArray, dims)
    d = ()
    for i in 1:ndims(A)
        if in(i,dims)
            if size(A,i) != 1
                error("squeezed dims must all be size 1")
            end
        else
            d = tuple(d..., size(A,i))
        end
    end
    reshape(A, d)
end

function copy!(dest::AbstractArray, src)
    i = 1
    for x in src
        dest[i] = x
        i += 1
    end
    return dest
end

# copy with minimal requirements on src
# if src is not an AbstractArray, moving to the offset might be O(n)
function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer=1)
    st = start(src)
    for j = 1:(soffs-1)
        _, st = next(src, st)
    end
    i = doffs
    while !done(src,st)
        val, st = next(src, st)
        dest[i] = val
        i += 1
    end
    return dest
end

# NOTE: this is to avoid ambiguity with the deprecation of
#   copy!(dest::AbstractArray, src, doffs::Integer)
# Remove this when that deprecation is removed.
function copy!(dest::AbstractArray, doffs::Integer, src::Integer)
    dest[doffs] = src
    return dest
end

# this method must be separate from the above since src might not have a length
function copy!(dest::AbstractArray, doffs::Integer, src, soffs::Integer, n::Integer)
    n == 0 && return dest
    st = start(src)
    for j = 1:(soffs-1)
        _, st = next(src, st)
    end
    for i = doffs:(doffs+n-1)
        done(src,st) && throw(BoundsError())
        val, st = next(src, st)
        dest[i] = val
    end
    return dest
end

# if src is an AbstractArray and a source offset is passed, use indexing
function copy!(dest::AbstractArray, doffs::Integer, src::AbstractArray, soffs::Integer, n::Integer=length(src))
    for i = 0:(n-1)
        dest[doffs+i] = src[soffs+i]
    end
    return dest
end

copy(a::AbstractArray) = copy!(similar(a), a)
copy(a::AbstractArray{None}) = a # cannot be assigned into so is immutable

function copy!{R,S}(B::AbstractMatrix{R}, ir_dest::Range{Int}, jr_dest::Range{Int}, A::AbstractMatrix{S}, ir_src::Range{Int}, jr_src::Range{Int})
    if length(ir_dest) != length(ir_src) || length(jr_dest) != length(jr_src)
        error("source and destination must have same size")
    end
    checkbounds(B, ir_dest, jr_dest)
    checkbounds(A, ir_src, jr_src)
    jdest = first(jr_dest)
    for jsrc in jr_src
        idest = first(ir_dest)
        for isrc in ir_src
            B[idest,jdest] = A[isrc,jsrc]
            idest += step(ir_dest)
        end
        jdest += step(jr_dest)
    end
    return B
end

function copy_transpose!{R,S}(B::AbstractMatrix{R}, ir_dest::Range{Int}, jr_dest::Range{Int}, A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int})
    if length(ir_dest) != length(jr_src) || length(jr_dest) != length(ir_src)
        error("source and destination must have same size")
    end
    checkbounds(B, ir_dest, jr_dest)
    checkbounds(A, ir_src, jr_src)
    idest = first(ir_dest)
    for jsrc in jr_src
        jdest = first(jr_dest)
        for isrc in ir_src
            B[idest,jdest] = A[isrc,jsrc]
            jdest += step(jr_dest)
        end
        idest += step(ir_dest)
    end
    return B
end

zero{T}(x::AbstractArray{T}) = fill!(similar(x), zero(T))

## iteration support for arrays as ranges ##

start(a::AbstractArray) = 1
next(a::AbstractArray,i) = (a[i],i+1)
done(a::AbstractArray,i) = (i > length(a))
isempty(a::AbstractArray) = (length(a) == 0)

## Conversions ##

for (f,t) in ((:char,   Char),
              (:int,    Int),
              (:int8,   Int8),
              (:int16,  Int16),
              (:int32,  Int32),
              (:int64,  Int64),
              (:int128, Int128),
              (:uint,   Uint),
              (:uint8,  Uint8),
              (:uint16, Uint16),
              (:uint32, Uint32),
              (:uint64, Uint64),
              (:uint128,Uint128))
    @eval begin
        ($f)(x::AbstractArray{$t}) = x
        ($f)(x::AbstractArray{$t}) = x

        function ($f)(x::AbstractArray)
            y = similar(x,$t)
            i = 1
            for e in x
                y[i] = ($f)(e)
                i += 1
            end
            y
        end
    end
end

for (f,t) in ((:integer, Integer),
              (:unsigned, Unsigned))
    @eval begin
        ($f){T<:$t}(x::AbstractArray{T}) = x
        ($f){T<:$t}(x::AbstractArray{T}) = x

        function ($f)(x::AbstractArray)
            y = similar(x,typeof(($f)(one(eltype(x)))))
            i = 1
            for e in x
                y[i] = ($f)(e)
                i += 1
            end
            y
        end
    end
end

bool(x::AbstractArray{Bool}) = x
bool(x::AbstractArray) = copy!(similar(x,Bool), x)

convert{T,S,N}(::Type{AbstractArray{T,N}}, A::Array{S,N}) = convert(Array{T,N}, A)

for (f,T) in ((:float16,    Float16),
              (:float32,    Float32),
              (:float64,    Float64),
              (:complex64,  Complex64),
              (:complex128, Complex128))
    @eval ($f){S,N}(x::AbstractArray{S,N}) = convert(AbstractArray{$T,N}, x)
end

float{T<:FloatingPoint}(x::AbstractArray{T}) = x
complex{T<:Complex}(x::AbstractArray{T}) = x

float{T,N}(x::AbstractArray{T,N}) = convert(AbstractArray{typeof(float(one(T))),N},x)
complex{T,N}(x::AbstractArray{T,N}) = convert(AbstractArray{typeof(complex(one(T))),N}, x)

full(x::AbstractArray) = x

## range conversions ##

for fn in _numeric_conversion_func_names
    @eval begin
        $fn(r::StepRange) = $fn(r.start):$fn(r.step):$fn(last(r))
        $fn(r::UnitRange) = $fn(r.start):$fn(last(r))
    end
end

for fn in (:float,:float16,:float32,:float64)
    @eval begin
        $fn(r::FloatRange) = FloatRange($fn(r.start), $fn(r.step), r.len, $fn(r.divisor))
    end
end

## Unary operators ##

conj{T<:Real}(x::AbstractArray{T}) = x
conj!{T<:Real}(x::AbstractArray{T}) = x

real{T<:Real}(x::AbstractArray{T}) = x
imag{T<:Real}(x::AbstractArray{T}) = zero(x)

+{T<:Number}(x::AbstractArray{T}) = x
*{T<:Number}(x::AbstractArray{T}) = x

## Binary arithmetic operators ##

*(A::Number, B::AbstractArray) = A .* B
*(A::AbstractArray, B::Number) = A .* B

/(A::AbstractArray, B::Number) = A ./ B

\(A::Number, B::AbstractArray) = B ./ A

## Indexing: getindex ##

getindex(t::AbstractArray, i::Real) = error("indexing not defined for ", typeof(t))

# linear indexing with a single multi-dimensional index
function getindex(A::AbstractArray, I::AbstractArray)
    x = similar(A, size(I))
    for i=1:length(I)
        x[i] = A[I[i]]
    end
    return x
end

# index A[:,:,...,i,:,:,...] where "i" is in dimension "d"
# TODO: more optimized special cases
slicedim(A::AbstractArray, d::Integer, i) =
    A[[ n==d ? i : (1:size(A,n)) for n in 1:ndims(A) ]...]

function flipdim(A::AbstractVector, d::Integer)
    d > 0 || error("dimension to flip must be positive")
    d == 1 || return copy(A)
    reverse(A)
end

function flipdim(A::AbstractArray, 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
    alli = [ 1:size(B,n) for n in 1:nd ]
    for i = 1:sd
        B[[ n==d ? sd+1-i : alli[n] for n in 1:nd ]...] = slicedim(A, d, i)
    end
    return B
end

flipud(A::AbstractArray) = flipdim(A, 1)
fliplr(A::AbstractArray) = flipdim(A, 2)

circshift(a, shiftamt::Real) = circshift(a, [integer(shiftamt)])
function circshift(a, shiftamts)
    n = ndims(a)
    I = cell(n)
    for i=1:n
        s = size(a,i)
        d = i<=length(shiftamts) ? shiftamts[i] : 0
        I[i] = d==0 ? (1:s) : mod([-d:s-1-d], s).+1
    end
    a[I...]::typeof(a)
end

## Indexing: setindex! ##

# 1-d indexing is assumed defined on subtypes
setindex!(t::AbstractArray, x, i::Real) =
    error("setindex! not defined for ",typeof(t))
setindex!(t::AbstractArray, x) = throw(MethodError(setindex!, (t, x)))

## Indexing: handle more indices than dimensions if "extra" indices are 1

# Don't require vector/matrix subclasses to implement more than 1/2 indices,
# respectively, by handling the extra dimensions in AbstractMatrix.

function getindex(A::AbstractVector, i1,i2,i3...)
    if i2*prod(i3) != 1
        throw(BoundsError())
    end
    A[i1]
end
function getindex(A::AbstractMatrix, i1,i2,i3,i4...)
    if i3*prod(i4) != 1
        throw(BoundsError())
    end
    A[i1,i2]
end

function setindex!(A::AbstractVector, x, i1,i2,i3...)
    if i2*prod(i3) != 1
        throw(BoundsError())
    end
    A[i1] = x
end
function setindex!(A::AbstractMatrix, x, i1,i2,i3,i4...)
    if i3*prod(i4) != 1
        throw(BoundsError())
    end
    A[i1,i2] = x
end

## get (getindex with a default value) ##

typealias RangeVecIntList{A<:AbstractVector{Int}} Union((Union(Range, AbstractVector{Int})...), AbstractVector{UnitRange{Int}}, AbstractVector{Range{Int}}, AbstractVector{A})

get(A::AbstractArray, i::Integer, default) = in_bounds(length(A), i) ? A[i] : default
get(A::AbstractArray, I::(), default) = similar(A, typeof(default), 0)
get(A::AbstractArray, I::Dims, default) = in_bounds(size(A), I...) ? A[I...] : default

function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::Union(Range, AbstractVector{Int}), default::T)
    ind = findin(I, 1:length(A))
    X[ind] = A[I[ind]]
    X[1:first(ind)-1] = default
    X[last(ind)+1:length(X)] = default
    X
end

get(A::AbstractArray, I::Range, default) = get!(similar(A, typeof(default), length(I)), A, I, default)

function get!{T}(X::AbstractArray{T}, A::AbstractArray, I::RangeVecIntList, default::T)
    fill!(X, default)
    dst, src = indcopy(size(A), I)
    X[dst...] = A[src...]
    X
end

get(A::AbstractArray, I::RangeVecIntList, default) = get!(similar(A, typeof(default), map(length, I)...), A, I, default)

## Concatenation ##

promote_eltype() = None
promote_eltype(v1, vs...) = promote_type(eltype(v1), promote_eltype(vs...))

#TODO: ERROR CHECK
cat(catdim::Integer) = Array(None, 0)

vcat() = Array(None, 0)
hcat() = Array(None, 0)

## cat: special cases
hcat{T}(X::T...)         = T[ X[j] for i=1, j=1:length(X) ]
hcat{T<:Number}(X::T...) = T[ X[j] for i=1, j=1:length(X) ]
vcat{T}(X::T...)         = T[ X[i] for i=1:length(X) ]
vcat{T<:Number}(X::T...) = T[ X[i] for i=1:length(X) ]

function vcat(X::Number...)
    T = None
    for x in X
        T = promote_type(T,typeof(x))
    end
    hvcat_fill(Array(T,length(X)), X)
end

function hcat(X::Number...)
    T = None
    for x in X
        T = promote_type(T,typeof(x))
    end
    hvcat_fill(Array(T,1,length(X)), X)
end

function hcat{T}(V::AbstractVector{T}...)
    height = length(V[1])
    for j = 2:length(V)
        if length(V[j]) != height; error("vector must have same lengths"); end
    end
    [ V[j][i]::T for i=1:length(V[1]), j=1:length(V) ]
end

function vcat{T}(V::AbstractVector{T}...)
    n = 0
    for Vk in V
        n += length(Vk)
    end
    a = similar(full(V[1]), n)
    pos = 1
    for k=1:length(V)
        Vk = V[k]
        for i=1:length(Vk)
            a[pos] = Vk[i]
            pos += 1
        end
    end
    a
end

function hcat{T}(A::AbstractVecOrMat{T}...)
    nargs = length(A)
    nrows = size(A[1], 1)
    ncols = 0
    dense = true
    for j = 1:nargs
        Aj = A[j]
        dense &= isa(Aj,Array)
        nd = ndims(Aj)
        ncols += (nd==2 ? size(Aj,2) : 1)
        if size(Aj, 1) != nrows; error("number of rows must match"); end
    end
    B = similar(full(A[1]), nrows, ncols)
    pos = 1
    if dense
        for k=1:nargs
            Ak = A[k]
            n = length(Ak)
            copy!(B, pos, Ak, 1, n)
            pos += n
        end
    else
        for k=1:nargs
            Ak = A[k]
            p1 = pos+(isa(Ak,AbstractMatrix) ? size(Ak, 2) : 1)-1
            B[:, pos:p1] = Ak
            pos = p1+1
        end
    end
    return B
end

function vcat{T}(A::AbstractMatrix{T}...)
    nargs = length(A)
    nrows = sum(a->size(a, 1), A)::Int
    ncols = size(A[1], 2)
    for j = 2:nargs
        if size(A[j], 2) != ncols; error("number of columns must match"); end
    end
    B = similar(full(A[1]), nrows, ncols)
    pos = 1
    for k=1:nargs
        Ak = A[k]
        p1 = pos+size(Ak,1)-1
        B[pos:p1, :] = Ak
        pos = p1+1
    end
    return B
end

## cat: general case

function cat(catdim::Integer, X...)
    nargs = length(X)
    dimsX = map((a->isa(a,AbstractArray) ? size(a) : (1,)), X)
    ndimsX = map((a->isa(a,AbstractArray) ? ndims(a) : 1), X)
    d_max = maximum(ndimsX)

    if catdim > d_max + 1
        for i=1:nargs
            if dimsX[1] != dimsX[i]
                error("all inputs must have same dimensions when concatenating along a higher dimension");
            end
        end
    elseif nargs >= 2
        for d=1:d_max
            if d == catdim; continue; end
            len = d <= ndimsX[1] ? dimsX[1][d] : 1
            for i = 2:nargs
                if len != (d <= ndimsX[i] ? dimsX[i][d] : 1)
                    error("mismatch in dimension ", d)
                end
            end
        end
    end

    cat_ranges = [ catdim <= ndimsX[i] ? dimsX[i][catdim] : 1 for i=1:nargs ]

    function compute_dims(d)
        if d == catdim
            if catdim <= d_max
                return sum(cat_ranges)
            else
                return nargs
            end
        else
            if d <= ndimsX[1]
                return dimsX[1][d]
            else
                return 1
            end
        end
    end

    ndimsC = max(catdim, d_max)
    dimsC = ntuple(ndimsC, compute_dims)::(Int...)
    typeC = promote_type(map(x->isa(x,AbstractArray) ? eltype(x) : typeof(x), X)...)
    C = similar(isa(X[1],AbstractArray) ? full(X[1]) : [X[1]], typeC, dimsC)

    range = 1
    for k=1:nargs
        nextrange = range+cat_ranges[k]
        cat_one = [ i != catdim ? (1:dimsC[i]) : (range:nextrange-1) for i=1:ndimsC ]
        C[cat_one...] = X[k]
        range = nextrange
    end
    return C
end

vcat(X...) = cat(1, X...)
hcat(X...) = cat(2, X...)

cat{T}(catdim::Integer, A::AbstractArray{T}...) = cat_t(catdim, T, A...)

cat(catdim::Integer, A::AbstractArray...) =
    cat_t(catdim, promote_eltype(A...), A...)

function cat_t(catdim::Integer, typeC, A::AbstractArray...)
    # ndims of all input arrays should be in [d-1, d]

    nargs = length(A)
    dimsA = map(size, A)
    ndimsA = map(ndims, A)
    d_max = maximum(ndimsA)

    if catdim > d_max + 1
        for i=1:nargs
            if dimsA[1] != dimsA[i]
                error("all inputs must have same dimensions when concatenating along a higher dimension");
            end
        end
    elseif nargs >= 2
        for d=1:d_max
            if d == catdim; continue; end
            len = d <= ndimsA[1] ? dimsA[1][d] : 1
            for i = 2:nargs
                if len != (d <= ndimsA[i] ? dimsA[i][d] : 1)
                    error("mismatch in dimension ", d)
                end
            end
        end
    end

    cat_ranges = [ catdim <= ndimsA[i] ? dimsA[i][catdim] : 1 for i=1:nargs ]

    function compute_dims(d)
        if d == catdim
            if catdim <= d_max
                return sum(cat_ranges)
            else
                return nargs
            end
        else
            if d <= ndimsA[1]
                return dimsA[1][d]
            else
                return 1
            end
        end
    end

    ndimsC = max(catdim, d_max)
    dimsC = ntuple(ndimsC, compute_dims)::(Int...)
    C = similar(full(A[1]), typeC, dimsC)

    range = 1
    for k=1:nargs
        nextrange = range+cat_ranges[k]
        cat_one = [ i != catdim ? (1:dimsC[i]) : (range:nextrange-1) for i=1:ndimsC ]
        C[cat_one...] = A[k]
        range = nextrange
    end
    return C
end

vcat(A::AbstractArray...) = cat(1, A...)
hcat(A::AbstractArray...) = cat(2, A...)

# 2d horizontal and vertical concatenation

function hvcat(nbc::Integer, as...)
    # nbc = # of block columns
    n = length(as)
    if mod(n,nbc) != 0
        error("all rows must have the same number of block columns")
    end
    nbr = div(n,nbc)
    hvcat(ntuple(nbr, i->nbc), as...)
end

function hvcat{T}(rows::(Int...), as::AbstractMatrix{T}...)
    nbr = length(rows)  # number of block rows

    nc = 0
    for i=1:rows[1]
        nc += size(as[i],2)
    end

    nr = 0
    a = 1
    for i = 1:nbr
        nr += size(as[a],1)
        a += rows[i]
    end

    out = similar(full(as[1]), T, nr, nc)

    a = 1
    r = 1
    for i = 1:nbr
        c = 1
        szi = size(as[a],1)
        for j = 1:rows[i]
            Aj = as[a+j-1]
            szj = size(Aj,2)
            if size(Aj,1) != szi
                error("mismatched height in block row ", i)
            end
            if c-1+szj > nc
                error("block row ", i, " has mismatched number of columns")
            end
            out[r:r-1+szi, c:c-1+szj] = Aj
            c += szj
        end
        if c != nc+1
            error("block row ", i, " has mismatched number of columns")
        end
        r += szi
        a += rows[i]
    end
    out
end

hvcat(rows::(Int...)) = []

function hvcat{T<:Number}(rows::(Int...), xs::T...)
    nr = length(rows)
    nc = rows[1]

    a = Array(T, nr, nc)
    k = 1
    for i=1:nr
        if nc != rows[i]
            error("row ", i, " has mismatched number of columns")
        end
        for j=1:nc
            a[i,j] = xs[k]
            k += 1
        end
    end
    a
end

function hvcat_fill(a, xs)
    k = 1
    nr, nc = size(a,1), size(a,2)
    for i=1:nr
        for j=1:nc
            a[i,j] = xs[k]
            k += 1
        end
    end
    a
end

function hvcat(rows::(Int...), xs::Number...)
    nr = length(rows)
    nc = rows[1]
    #error check
    for i = 2:nr
        if nc != rows[i]
            error("row ", i, " has mismatched number of columns")
        end
    end
    T = typeof(xs[1])
    for i=2:length(xs)
        T = promote_type(T,typeof(xs[i]))
    end
    hvcat_fill(Array(T, nr, nc), xs)
end

## Reductions and scans ##

function isequal(A::AbstractArray, B::AbstractArray)
    if A === B return true end
    if size(A) != size(B)
        return false
    end
    if isa(A,Range) != isa(B,Range)
        return false
    end
    for i = 1:length(A)
        if !isequal(A[i], B[i])
            return false
        end
    end
    return true
end

function lexcmp(A::AbstractArray, B::AbstractArray)
    nA, nB = length(A), length(B)
    for i = 1:min(nA, nB)
        res = lexcmp(A[i], B[i])
        res == 0 || return res
    end
    return cmp(nA, nB)
end

function (==)(A::AbstractArray, B::AbstractArray)
    if size(A) != size(B)
        return false
    end
    if isa(A,Range) != isa(B,Range)
        return false
    end
    for i = 1:length(A)
        if !(A[i]==B[i])
            return false
        end
    end
    return true
end

# Uses K-B-N summation
function cumsum_kbn{T<:FloatingPoint}(v::AbstractVector{T})
    n = length(v)
    r = similar(v, n)
    if n == 0; return r; end

    s = r[1] = v[1]
    c = zero(T)
    for i=2:n
        vi = v[i]
        t = s + vi
        if abs(s) >= abs(vi)
            c += ((s-t) + vi)
        else
            c += ((vi-t) + s)
        end
        s = t
        r[i] = s+c
    end
    return r
end

# Uses K-B-N summation
function cumsum_kbn{T<:FloatingPoint}(A::AbstractArray{T}, 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)
    C = similar(A)

    for i = 1:length(A)
        if div(i-1, axis_stride) % axis_size == 0
            B[i] = A[i]
            C[i] = zero(T)
        else
            s = B[i-axis_stride]
            Ai = A[i]
            B[i] = t = s + Ai
            if abs(s) >= abs(Ai)
                C[i] = C[i-axis_stride] + ((s-t) + Ai)
            else
                C[i] = C[i-axis_stride] + ((Ai-t) + s)
            end
        end
    end

    return B + C
end

## ipermutedims in terms of permutedims ##

function ipermutedims(A::AbstractArray,perm)
    iperm = Array(Int,length(perm))
    for i = 1:length(perm)
	iperm[perm[i]] = i
    end
    return permutedims(A,iperm)
end

## Other array functions ##

# fallback definition of hvcat in terms of hcat and vcat
function hvcat(rows::(Int...), as...)
    nbr = length(rows)  # number of block rows
    rs = cell(nbr)
    a = 1
    for i = 1:nbr
        rs[i] = hcat(as[a:a-1+rows[i]]...)
        a += rows[i]
    end
    vcat(rs...)
end

function repmat(a::AbstractVecOrMat, m::Int, n::Int=1)
    o, p = size(a,1), size(a,2)
    b = similar(a, o*m, p*n)
    for j=1:n
        d = (j-1)*p+1
        R = d:d+p-1
        for i=1:m
            c = (i-1)*o+1
            b[c:c+o-1, R] = a
        end
    end
    return b
end

function repmat(a::AbstractVector, m::Int)
    o = length(a)
    b = similar(a, o*m)
    for i=1:m
        c = (i-1)*o+1
        b[c:c+o-1] = a
    end
    return b
end

sub2ind(dims) = 1
sub2ind(dims, i::Integer) = int(i)
sub2ind(dims, i::Integer, j::Integer) = sub2ind(dims, int(i), int(j))
sub2ind(dims, i::Int, j::Int) = (j-1)*dims[1] + i
sub2ind(dims, i0::Integer, i1::Integer, i2::Integer) = sub2ind(dims, int(i0),int(i1),int(i2))
sub2ind(dims, i0::Int, i1::Int, i2::Int) =
    i0 + dims[1]*((i1-1) + dims[2]*(i2-1))
sub2ind(dims, i0::Integer, i1::Integer, i2::Integer, i3::Integer) =
    sub2ind(dims, int(i0),int(i1),int(i2),int(i3))
sub2ind(dims, i0::Int, i1::Int, i2::Int, i3::Int) =
    i0 + dims[1]*((i1-1) + dims[2]*((i2-1) + dims[3]*(i3-1)))

function sub2ind(dims, I::Integer...)
    ndims = length(dims)
    index = int(I[1])
    stride = 1
    for k=2:ndims
        stride = stride * dims[k-1]
        index += (int(I[k])-1) * stride
    end
    return index
end

function sub2ind{T<:Integer}(dims::Array{T}, sub::Array{T})
    ndims = length(dims)
    ind = sub[1]
    stride = 1
    for k in 2:ndims
        stride = stride * dims[k - 1]
        ind += (sub[k] - 1) * stride
    end
    return ind
end

sub2ind{T<:Integer}(dims, I::AbstractVector{T}...) =
    [ sub2ind(dims, map(X->X[i], I)...)::Int for i=1:length(I[1]) ]

function ind2sub(dims::(Integer,Integer...), ind::Int)
    ndims = length(dims)
    stride = dims[1]
    for i=2:ndims-1
        stride *= dims[i]
    end

    sub = ()
    for i=(ndims-1):-1:1
        rest = rem(ind-1, stride) + 1
        sub = tuple(div(ind - rest, stride) + 1, sub...)
        ind = rest
        stride = div(stride, dims[i])
    end
    return tuple(ind, sub...)
end

ind2sub(dims::(Integer...), ind::Integer) = ind2sub(dims, int(ind))
ind2sub(dims::(), ind::Integer) = ind==1 ? () : throw(BoundsError())
ind2sub(dims::(Integer,), ind::Int) = (ind,)
ind2sub(dims::(Integer,Integer), ind::Int) =
    (rem(ind-1,dims[1])+1, div(ind-1,dims[1])+1)
ind2sub(dims::(Integer,Integer,Integer), ind::Int) =
    (rem(ind-1,dims[1])+1, div(rem(ind-1,dims[1]*dims[2]), dims[1])+1,
     div(rem(ind-1,dims[1]*dims[2]*dims[3]), dims[1]*dims[2])+1)

function ind2sub{T<:Integer}(dims::(Integer,Integer...), ind::AbstractVector{T})
    n = length(dims)
    l = length(ind)
    t = ntuple(n, x->Array(Int, l))
    for i = 1:l
        s = ind2sub(dims, ind[i])
        for j = 1:n
            t[j][i] = s[j]
        end
    end
    return t
end

function ind2sub!{T<:Integer}(sub::Array{T}, dims::Array{T}, ind::T)
    ndims = length(dims)
    stride = dims[1]
    for i in 2:(ndims - 1)
        stride *= dims[i]
    end
    for i in (ndims - 1):-1:1
        rest = rem1(ind, stride)
        sub[i + 1] = div(ind - rest, stride) + 1
        ind = rest
        stride = div(stride, dims[i])
    end
    sub[1] = ind
    return
end

# Generalized repmat
function repeat{T}(A::Array{T};
                   inner::Array{Int} = ones(Int, ndims(A)),
                   outer::Array{Int} = ones(Int, ndims(A)))
    ndims_in = ndims(A)
    length_inner = length(inner)
    length_outer = length(outer)
    ndims_out = max(ndims_in, length_inner, length_outer)

    if length_inner < ndims_in || length_outer < ndims_in
        msg = "inner/outer repetitions must be set for all input dimensions"
        throw(ArgumentError(msg))
    end

    size_in = Array(Int, ndims_in)
    size_out = Array(Int, ndims_out)
    inner_size_out = Array(Int, ndims_out)

    for i in 1:ndims_in
        size_in[i] = size(A, i)
    end
    for i in 1:ndims_out
        t1 = ndims_in < i ? 1 : size_in[i]
        t2 = length_inner < i ? 1 : inner[i]
        t3 = length_outer < i ? 1 : outer[i]
        size_out[i] = t1 * t2 * t3
        inner_size_out[i] = t1 * t2
    end

    indices_in = Array(Int, ndims_in)
    indices_out = Array(Int, ndims_out)

    length_out = prod(size_out)
    R = Array(T, size_out...)

    for index_out in 1:length_out
        ind2sub!(indices_out, size_out, index_out)
        for t in 1:ndims_in
            # "Project" outer repetitions into inner repetitions
            indices_in[t] = mod1(indices_out[t], inner_size_out[t])
            # Find inner repetitions using flooring division
            if inner[t] != 1
                indices_in[t] = fld1(indices_in[t], inner[t])
            end
        end
        index_in = sub2ind(size_in, indices_in)
        R[index_out] = A[index_in]
    end

    return R
end

## iteration utilities ##

# slow, but useful
function cartesianmap(body, t::(Int...), it...)
    idx = length(t)-length(it)
    if idx == 1
        for i = 1:t[1]
            body(i, it...)
        end
    elseif idx == 2
        for j = 1:t[2]
            for i = 1:t[1]
                body(i, j, it...)
            end
        end
    elseif idx > 1
        for j = 1:t[idx]
            for i = 1:t[idx-1]
                cartesianmap(body, t, i, j, it...)
            end
        end
    else
        throw(ArgumentError("cartesianmap"))
    end
end

cartesianmap(body, t::()) = (body(); nothing)

function cartesianmap(body, t::(Int,))
    for i = 1:t[1]
        body(i)
    end
end

function cartesianmap(body, t::(Int,Int))
    for j = 1:t[2]
        for i = 1:t[1]
            body(i,j)
        end
    end
end

function cartesianmap(body, t::(Int,Int,Int))
    for k = 1:t[3]
        for j = 1:t[2]
            for i = 1:t[1]
                body(i,j,k)
            end
        end
    end
end

##
# generic map on any iterator
function map(f::Callable, iters...)
    result = {}
    len = length(iters)
    states = [start(iters[idx]) for idx in 1:len]
    nxtvals = cell(len)
    cont = true
    for idx in 1:len
        done(iters[idx], states[idx]) && (cont = false; break)
    end
    while cont
        for idx in 1:len
            nxtvals[idx],states[idx] = next(iters[idx], states[idx])
        end
        push!(result, f(nxtvals...))
        for idx in 1:len
            done(iters[idx], states[idx]) && (cont = false; break)
        end
    end
    result
end

## map over arrays ##

## transform any set of dimensions
## dims specifies which dimensions will be transformed. for example
## dims==1:2 will call f on all slices A[:,:,...]
mapslices(f::Function, A::AbstractArray, dims) = mapslices(f, A, [dims...])
function mapslices(f::Function, A::AbstractArray, dims::AbstractVector)
    if isempty(dims)
        return map(f,A)
    end

    dimsA = [size(A)...]
    ndimsA = ndims(A)
    alldims = [1:ndimsA]

    otherdims = setdiff(alldims, dims)

    idx = cell(ndimsA)
    fill!(idx, 1)
    Asliceshape = tuple(dimsA[dims]...)
    itershape   = tuple(dimsA[otherdims]...)
    for d in dims
        idx[d] = 1:size(A,d)
    end

    r1 = f(reshape(A[idx...], Asliceshape))

    # determine result size and allocate
    Rsize = copy(dimsA)
    # TODO: maybe support removing dimensions
    if !isa(r1, AbstractArray) || ndims(r1) == 0
        r1 = [r1]
    end
    Rsize[dims] = [size(r1)...; ones(Int,max(0,length(dims)-ndims(r1)))]
    R = similar(r1, tuple(Rsize...))

    ridx = cell(ndims(R))
    fill!(ridx, 1)
    for d in dims
        ridx[d] = 1:size(R,d)
    end

    R[ridx...] = r1

    first = true
    cartesianmap(itershape) do idxs...
        if first
            first = false
        else
            ia = [idxs...]
            idx[otherdims] = ia
            ridx[otherdims] = ia
            R[ridx...] = f(reshape(A[idx...], Asliceshape))
        end
    end

    return R
end

## 1 argument
function map_to!(f::Callable, first, dest::AbstractArray, A::AbstractArray)
    dest[1] = first
    for i=2:length(A)
        dest[i] = f(A[i])
    end
    return dest
end

function map(f::Callable, A::AbstractArray)
    if isempty(A); return {}; end
    first = f(A[1])
    dest = similar(A, typeof(first))
    return map_to!(f, first, dest, A)
end

## 2 argument
function map_to!(f::Callable, first, dest::AbstractArray, A::AbstractArray, B::AbstractArray)
    dest[1] = first
    for i=2:length(A)
        dest[i] = f(A[i], B[i])
    end
    return dest
end

function map(f::Callable, A::AbstractArray, B::AbstractArray)
    shp = promote_shape(size(A),size(B))
    if isempty(A)
        return similar(A, Any, shp)
    end
    first = f(A[1], B[1])
    dest = similar(A, typeof(first), shp)
    return map_to!(f, first, dest, A, B)
end

## N argument
function map_to!(f::Callable, first, dest::AbstractArray, As::AbstractArray...)
    n = length(As[1])
    i = 1
    ith = a->a[i]
    dest[1] = first
    for i=2:n
        dest[i] = f(map(ith, As)...)
    end
    return dest
end

function map(f::Callable, As::AbstractArray...)
    shape = mapreduce(size, promote_shape, As)
    if prod(shape) == 0
        return similar(As[1], Any, shape)
    end
    first = f(map(a->a[1], As)...)
    dest = similar(As[1], typeof(first), shape)
    return map_to!(f, first, dest, As...)
end

# multi-item push!, unshift! (built on top of type-specific 1-item version)
# (note: must not cause a dispatch loop when 1-item case is not defined)
push!(A) = A
push!(A, a, b) = push!(push!(A, a), b)
push!(A, a, b, c...) = push!(push!(A, a, b), c...)
unshift!(A) = A
unshift!(A, a, b) = unshift!(unshift!(A, b), a)
unshift!(A, a, b, c...) = unshift!(unshift!(A, c...), a, b)

ext: jl

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using Omega
using OmegaRemix

function softtest2(ω = 0.39)
  g(x::Bool) = x & x
  function f(ω)
    if Bool(ω > 0.35)
      baba = ω > 0.5
      # @show typeof(baba)
      g(baba)
    else
      2ω <= 0.2
    end
  end
  softapply(f, ω)
end

function f(ω)
  if Bool(ω > 0.35)
    (baba, err) = ω >ₛ 0.5
    # @show typeof(baba)
    (res, err2) = g(baba)

  else
    2ω <= 0.2
  end
end

softtest2()

function softtest3()
  g(x::Real)::Bool = x > 0.5
  softapply(g, 0.3)
end

softtest3()


text:
import SimpleDiffEq

@inline @inbounds function eom_lorenz_lorenz_bidir(u, p, t)
    c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃ = (p...,)
    x1, x2, x3, y1, y2, y3 = (u...,)

    dx1 = -a₁*(x1 - x2) + c_yx*(y1 - x1)
    dx2 = -x1*x3 + a₂*x1 - x2
    dx3 = x1*x2 - a₃*x3
    dy1 = -b₁*(y1 - y2) + c_xy*(x1 - y1)
    dy2 = -y1*y3 + b₂*y1 - y2
    dy3 = y1*y2 - b₃*y3

    return SVector{6}(dx1, dx2, dx3, dy1, dy2, dy3)
end

"""
    lorenz_lorenz_bidir(; u0 = rand(6),
        c_xy = 0.2, c_yx = 0.2,
        a₁ = 10, a₂ = 28, a₃ = 8/3,
        b₁ = 10, b₂ = 28, b₃ = 9/3) -> ContinuousDynamicalSystem

Initialise a bidirectionally coupled Lorenz-Lorenz system, where each
subsystem is a 3D Lorenz system [1]. Default values for the parameters
`a₁`, `a₂`, `a₃`, `b₁`, `b₂`, `b₃` are as in [1].

## Equations of motion

The dynamics is generated by the following vector field

```math
\\begin{aligned}
\\dot{x_1} &= -a_1 (x_1 - x_2) + c_{yx}(y_1 - x_1) \\\\
\\dot{x_2} &= -x_1 x_3 + a_2 x_1 - x_2 \\\\
\\dot{x_3} &= x_1 x_2 - a_3 x_3 \\\\
\\dot{y_1} &= -b_1 (y_1 - y_2) + c_{xy} (x_1 - y_1) \\\\
\\dot{y_2} &= -y_1 y_3 + b_2 y_1 - y_2 \\\\
\\dot{y_3} &= y_1 y_2 - b_3 y_3
\\end{aligned}
```

## References

1. Amigó, José M., and Yoshito Hirata. "Detecting directional couplings from
    multivariate flows by the joint distance distribution." Chaos: An
    Interdisciplinary Journal of Nonlinear Science 28.7 (2018): 075302.
"""
function lorenz_lorenz_bidir(; u0 = rand(6),
        c_xy = 0.2, c_yx = 0.2,
        a₁ = 10, a₂ = 28, a₃ = 8/3,
        b₁ = 10, b₂ = 28, b₃ = 9/3)
    ContinuousDynamicalSystem(eom_lorenz_lorenz_bidir, u0, [c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃])
end

function lorenzlorenz_bidir_trajectory(npts; sample_dt = 1, Ttr = 1000, dt = 0.1,
    c_xy = 0.1, c_yx = 0.1,
    u0 = rand(6),
    a₁ = 10, a₂ = 28, a₃ = 8/3,
    b₁ = 10, b₂ = 28, b₃ = 9/3)

    s = lorenz_lorenz_bidir(u0 = u0, c_xy = c_xy, c_yx = c_yx, a₁ = a₁, a₂ = a₂, a₃ = a₃, b₁ = b₁, b₂ = b₂, b₃ = b₃)

    # the system is recorded at times t0:dt:T
    T = npts*dt*sample_dt

    o = trajectory(s, T, dt = dt, Ttr = Ttr*dt, alg = SimpleDiffEq.SimpleATsit5())[1:sample_dt:end-1, :]
end

# For some initial conditions, the system wanders off and doesn't settle to an attractor. Create a function that loops until we get a good realization.
function good_lorenzlorenz_bidir_trajectory(npts;
        sample_dt = 1,
        dt = 0.1,
        c_xy = 0.1,
        c_yx = 0.1,
        Da₁ = Uniform(9.5, 10.5),
        Da₂ = Uniform(27, 29),
        Da₃ = Uniform(7.5/3, 8.5/3),
        Db₁ = Uniform(9.5, 10.5),
        Db₂ = Uniform(27, 29),
        Db₃ = Uniform(7.5/3, 8.5/3),
        a₁ = nothing,
        a₂ = nothing,
        a₃ = nothing,
        b₁ = nothing,
        b₂ = nothing,
        b₃ = nothing,
        u0 = rand(6),
        Ttr = 10000,
        n_maxtries = 300)

    n_tries = 0
    while n_tries <= n_maxtries
        a₁ == nothing ? a₁ = rand(Da₁) : nothing
        a₂ == nothing ? a₂ = rand(Da₂) : nothing
        a₃ == nothing ? a₃ = rand(Da₃) : nothing
        b₁ == nothing ? b₁ = rand(Db₁) : nothing
        b₂ == nothing ? b₂ = rand(Db₂) : nothing
        b₃ == nothing ? b₃ = rand(Db₃) : nothing

        pts = lorenzlorenz_bidir_trajectory(npts,
            sample_dt = sample_dt, dt = dt,
            c_xy = c_xy, c_yx = c_yx,
            Ttr = Ttr)

        M = Matrix(pts)

        if all(isfinite.(M)) && all(M .< 1e10) && count(M .≈ 0) < npts*0.1 && count(abs.(M) .< 1e-10) < npts*0.1 &&
            (count(abs.(M) .< 1e-12) < npts*0.1)
            return pts
        end
        println("no attractor found. trying with new initial condition and parameters")
        n_tries += 1
    end
end


text: global PROC; PROC="CPU"

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module SailRoute

include("uncertainty/discretization_error.jl")
include("route/domain.jl")
include("performance/polar.jl")
include("route/shortest_path.jl")
include("weather/load_weather.jl")

end # module


text:
module MatDeforMooneyRivlinADModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D, nstressstrain, nthermstrain
import FinEtoolsDeforLinear.MatDeforModule: AbstractMatDefor, stressvtot!, stressttov!, strainttov!, strainvdet, strainvtr
import ..MatDeforNonlinearModule: AbstractMatDeforNonlinear, totlag2currsymm!
using LinearAlgebra: Transpose, Diagonal, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
using LinearAlgebra: eigen, eigvals, norm, cholesky, cross, dot, log, diagm, det
using ForwardDiff: gradient, hessian

"""
	MatDeforMooneyRivlinAD{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear

Type for triaxial Mooney-Rivlin hyperelastic material.

"""
struct  MatDeforMooneyRivlinAD{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear
	mr::Type{MR} # model reduction type
	mass_density::FFlt # mass density
	E::FFlt # Young's modulus
	nu::FFlt # Poisson ratio
	_lambda::FFlt
	_mu::FFlt
	tangentmoduli!::MTAN # Function to return the tangent moduli matrix
	update!::MUPD # Function to update the material state
end

################################################################################
# 3-D solid model
################################################################################

"""
    MatDeforMooneyRivlinAD(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)

Create triaxial Mooney-Rivlin hyperelastic material.
"""
function MatDeforMooneyRivlinAD(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)
	_mI = diagm(0=>[1, 1, 1, 0.5, 0.5, 0.5]);
	_m1 = vec(FFlt[1 1 1 0 0 0]) ;
	_m1m1 = _m1*_m1';
	_I3 = [1.0 0 0; 0 1.0 0; 0 0 1.0]
	function strainenergy(Cv, MR, lambda, mu)
		trC = strain6vtr(MR, Cv)
		J = sqrt(strain6vdet(MR, Cv))
		lJ = log(J)
		return mu/2*(trC-3) - mu*lJ + lambda/2*(lJ)^2;
	end
	function tangentmoduli3d!(self::MatDeforMooneyRivlinAD, D::FFltMat, statev::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat, label::FInt)
		C = Fn1'*Fn1;
		Cv = fill(0.0, 6)
		strain3x3tto6v!(Cv, C)
		Dtotal = 4 .* hessian(Cv -> strainenergy(Cv, mr, self._lambda, self._mu), Cv);
		return totlag2currsymm!(D, Dtotal, Fn1)
	end
	function update3d!(self::MatDeforMooneyRivlinAD, statev::FFltVec, cauchy::FFltVec, output::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
		@assert length(cauchy) == nstressstrain(self.mr)
		# Note: the code below allocates. This may not be a big deal when the
		# code gets called only a few times, but in explicit dynamics
		# situations  this would be very expensive. The material type needs to
		# provide these temporary arrays as buffers.
		C = Fn1'*Fn1;
		Cv = fill(0.0, 6)
		strain3x3tto6v!(Cv, C)
		Sv = 2 * gradient(Cv -> strainenergy(Cv, mr, self._lambda, self._mu), Cv);
		S = fill(0.0, 3, 3)
		stress6vto3x3t!(S, Sv)
		cauchyt = Fn1*(S/det(Fn1))*Fn1'; # Cauchy stress
		stress3x3tto6v!(cauchy, cauchyt)
		if quantity == :nothing
			#Nothing to be copied to the output array
		elseif quantity == :cauchy || quantity == :Cauchy
			(length(output) >= 6) || (output = zeros(6)) # make sure we can store it
			copyto!(output, cauchy);
		elseif quantity == :pressure || quantity == :Pressure
			output[1]  =  -sum(cauchy[1:3])/3.
		elseif quantity == :princCauchy || quantity == :princcauchy
			ep = eigen(sigma);
			(length(output) >= 3) || (output = zeros(3)) # make sure we can store it
			copyto!(output,  sort(ep.values, rev=true));
		elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
			s1=cauchy[1]; s2=cauchy[2]; s3=cauchy[3];
			s4=cauchy[4]; s5=cauchy[5]; s6=cauchy[6];
			(length(output) >= 1) || (output = zeros(1)) # make sure we can store it
			output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
		end
		return output
	end
	_lambda = E * nu / (1 + nu) / (1 - 2*(nu));
	_mu = E / 2. / (1+nu);
	return MatDeforMooneyRivlinAD(mr, mass_density, E, nu, _lambda, _mu,
		tangentmoduli3d!, update3d!)
end

"""
    MatDeforMooneyRivlinAD(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}

Create Mooney-Rivlin isotropic elastic material.

The mass density is the default value.
"""
function MatDeforMooneyRivlinAD(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}
	mass_density = 1.0
	return MatDeforMooneyRivlinAD(mr, mass_density, E, nu)
end

end


text:
macro checktoplevel(ctx, expr)
    quote
        if !istoplevel($(esc(ctx)).current)
            msg($(esc(ctx)), :E100, "$($(esc(expr))) expression must be at top level")
            return
        end
    end
end

macro checkisa(ctx, var, typ)
    quote
        if !isa($(esc(var)), $(esc(typ)))
            msg($(esc(ctx)), :E101, $(esc(var)), "this expression must be a $($(esc(typ)))")
            return
        end
    end
end

function lintexpr(ex::Symbol, ctx::LintContext)
    registersymboluse(ex, ctx)
end

function lintexpr(ex::QuoteNode, ctx::LintContext)
    if typeof(ex.value) == Expr
        ctx.quoteLvl += 1
        lintexpr(ex.value, ctx)
        ctx.quoteLvl -= 1
    end
end

function lintexpr(ex::Expr, ctx::LintContext)
    # TODO: reenable linthelpers
    if ex.head == :line
        # ignore line numer nodes
        return
    elseif ex.head == :block
        lintblock(ex, ctx)
    elseif ex.head == :quote
        ctx.quoteLvl += 1
        lintexpr(ex.args[1], ctx)
        ctx.quoteLvl -= 1
    elseif ex.head == :if
        lintifexpr(ex, ctx)
    elseif ex.head == :(=) && typeof(ex.args[1])==Expr && ex.args[1].head == :call
        lintfunction(ex, ctx)
    elseif expand_assignment(ex) !== nothing
        ea = expand_assignment(ex)
        lintassignment(Expr(:(=), ea[1], ea[2]), ctx)
    elseif ex.head == :local
        lintlocal(ex, ctx)
    elseif ex.head == :global
        lintglobal(ex, ctx)
    elseif ex.head == :const
        if typeof(ex.args[1]) == Expr && ex.args[1].head == :(=)
            lintassignment(ex.args[1], ctx; isConst = true)
        end
    elseif ex.head == :module
        lintmodule(ex, ctx)
    elseif ex.head == :export
        lintexport(ex, ctx)
    elseif isexpr(ex, [:import, :using, :importall])
        lintimport(ex, ctx)
    elseif ex.head == :comparison # only the odd indices
        for i in 1:2:length(ex.args)
            # comparison like match != 0:-1 is allowed, and shouldn't trigger lint warnings
            if Meta.isexpr(ex.args[i], :(:)) && length(ex.args[i].args) == 2 &&
                typeof(ex.args[i].args[1]) <: Real &&
                typeof(ex.args[i].args[2]) <: Real
                continue
            else
                lintexpr(ex.args[i], ctx)
            end
        end
        lintcomparison(ex, ctx)
    elseif ex.head == :type
        linttype(ex, ctx)
    elseif ex.head == :typealias
        # TODO: deal with X{T} = Y assignments, also const X = Y
        linttypealias(ex, ctx)
    elseif ex.head == :abstract
        lintabstract(ex, ctx)
    elseif ex.head == :bitstype
        lintbitstype(ex, ctx)
    elseif ex.head == :(->)
        lintlambda(ex, ctx)
    elseif ex.head == :($) && ctx.quoteLvl > 0 # an unquoted node inside a quote node
        ctx.quoteLvl -= 1
        lintexpr(ex.args[1], ctx)
        ctx.quoteLvl += 1
    elseif ex.head == :function
        lintfunction(ex, ctx)
    elseif ex.head == :stagedfunction
        lintfunction(ex, ctx, isstaged=true)
    elseif ex.head == :macrocall && ex.args[1] == Symbol("@generated")
        lintfunction(ex.args[2], ctx, isstaged=true)
    elseif ex.head == :macro
        lintmacro(ex, ctx)
    elseif ex.head == :macrocall
        lintmacrocall(ex, ctx)
    elseif ex.head == :call
        lintfunctioncall(ex, ctx)
    elseif ex.head == :(:) # TODO(felipe) check for `Colon()`
        lintrange(ex, ctx)
    elseif ex.head == :(::) # type assert/convert
        lintexpr(ex.args[1], ctx)
    elseif ex.head == :(.) # a.b
        lintexpr(ex.args[1], ctx)
    elseif ex.head == :ref # it could be a ref a[b], or an array Int[1,2]
        lintref(ex, ctx)
    elseif ex.head == :typed_vcat# it could be a ref a[b], or an array Int[1,2]
        linttyped_vcat(ex, ctx)
    elseif ex.head == :vcat
        lintvcat(ex, ctx)
    elseif ex.head == :vect # 0.4
        lintvect(ex, ctx)
    elseif ex.head == :hcat
        linthcat(ex, ctx)
    elseif ex.head == :typed_hcat
        linttyped_hcat(ex, ctx)
    elseif ex.head == :cell1d
        lintcell1d(ex, ctx)
    elseif ex.head == :while
        lintwhile(ex, ctx)
    elseif ex.head == :for
        lintfor(ex, ctx)
    elseif ex.head == :let
        lintlet(ex, ctx)
    elseif ex.head in (:comprehension, :dict_comprehension, :generator)
        lintgenerator(ex, ctx; typed = false)
    elseif ex.head in (:typed_comprehension, :typed_dict_comprehension)
        lintgenerator(ex, ctx; typed = true)
    elseif ex.head == :try
        linttry(ex, ctx)
    elseif ex.head == :curly # e.g. Ptr{T}
        lintcurly(ex, ctx)
    elseif ex.head in [:(&&), :(||)]
        lintboolean(ex.args[1], ctx)
        lintexpr(ex.args[2], ctx) # do not enforce boolean. e.g. b==1 || error("b must be 1!")
    elseif ex.head == :incomplete
        msg(ctx, :E112, ex.args[1])
    else
        for sube in ex.args
            lintexpr(sube, ctx)
        end
    end
end

# no-op fallback for other kinds of expressions (e.g. LineNumberNode) that we
# don’t care to handle
lintexpr(::Any, ::LintContext) = return


text:
using DataDeps

register(DataDep(
  "airfoil-self-noise",
  "https://archive.ics.uci.edu/ml/datasets/Airfoil+Self-Noise",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/00291/airfoil_self_noise.dat",
  "74c75fd71783f1e6b71f8a622b993dc592897a97cd689c5090a07147a1b097b3",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=6,
        feature_indices=1:5,
        delim='\t',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
commands = readlines("input.txt")
x_position = 0
depth = 0
aim = 0

for command in commands
    local direction, dist = split(command)
    local distance = parse(Int64, dist)
    if direction == "forward"
        global x_position += distance
        global depth += distance * aim
    elseif direction == "down"
        global aim += distance
    elseif direction == "up"
        global aim -= distance
    else
        println("Unknown direction! :( ")
    end
end

println("x = ",x_position)
println("depth = ",depth)
println("x*depth = ",x_position*depth)


text:
export set_coords
function set_coords(c::GraphicsContext, x, y, w, h, l, r, t, b)
    Base.depwarn("set_coords is deprecated; use set_coordinates or inner_canvas instead", :set_coords)
    inner_canvas(c::GraphicsContext, x, y, w, h, l, r, t, b)
end
function set_coords(c::GraphicsContext, device::BoundingBox, user::BoundingBox)
    Base.depwarn("set_coords is deprecated; use set_coordinates or inner_canvas instead", :set_coords)
    inner_canvas(c::GraphicsContext, device, user)
end
function set_coords(c::GraphicsContext, user::BoundingBox)
    Base.depwarn("set_coords is deprecated; use set_coordinates instead", :set_coords)
    inner_canvas(c::GraphicsContext, BoundingBox(0, width(c), 0, height(c)), user)
end


text:
using ParticleDA
using LinearAlgebra, Test, HDF5, Random, YAML
using MPI
using StableRNGs
using FFTW
using GaussianRandomFields: Matern

using ParticleDA: FilterParameters

include(joinpath(@__DIR__, "model", "model.jl"))

using .Model, .Model.LLW2d
using .Model: ModelParameters

@testset "LLW2d" begin
    dx = dy = 2e3

    ### set_stations!
    ist = Vector{Int}(undef, 4)
    jst = Vector{Int}(undef, 4)
    LLW2d.set_stations!(ist, jst, 20e3, 20e3, 150e3, 150e3, dx, dy)
    @test ist == [75, 75, 85, 85]
    @test jst == [75, 85, 75, 85]
    ist = rand(Int, 9)
    jst = rand(Int, 9)
    LLW2d.set_stations!(ist, jst, 20e3, 20e3, 150e3, 150e3, dx, dy)

    @test ist == [75, 75, 75, 85, 85, 85, 95, 95, 95]
    @test jst == [75, 85, 95, 75, 85, 95, 75, 85, 95]

    ### initheight!
    nx = 2
    ny = 2
    eta = ones(nx, ny)
    ocean_depth  = ones(nx, ny)
    peak_position = [floor(Int, nx/4) * dx, floor(Int, ny/4) * dy]
    LLW2d.initheight!(eta, ocean_depth, dx, dy, 3e4, 1.0, peak_position)
    @test eta ≈ [0.978266982572228  0.9463188389826958;
                 0.9463188389826958 0.9154140546161575]
    eta = ones(2, 2)
    ocean_depth  = zeros(2, 2)
    LLW2d.initheight!(eta, ocean_depth, dx, dy, 3e4, 1.0, peak_position)
    @test eta ≈ zeros(2, 2)

    # timestep.  TODO: add real tests.  So far we're just making sure code won't
    # crash
    n = 200
    eta1 = rand(n, n)
    mm1 = rand(n, n)
    nn1 = rand(n, n)
    eta0 = rand(n, n)
    mm0 = rand(n, n)
    nn0 = rand(n, n)
    y_averaged_depth = rand(n, n)
    x_averaged_depth = rand(n, n)
    land_filter_m = rand(n, n)
    land_filter_n = rand(n, n)
    land_filter_e = rand(n,n)
    absorbing_boundary = rand(n,n)
    model_matrices = LLW2d.Matrices(absorbing_boundary, ocean_depth,
                                    x_averaged_depth, y_averaged_depth,
                                    land_filter_m, land_filter_n, land_filter_e)
    dxeta = Matrix{Float64}(undef, n, n)
    dyeta = Matrix{Float64}(undef, n, n)
    LLW2d.timestep!(dxeta, dyeta, eta1, mm1, nn1, eta0, mm0, nn0, model_matrices, dx, dy, 1)

    # setup.  TODO: add real tests.  So far we're just making sure code won't
    # crash
    LLW2d.setup(n, n, 3e4, 0.1, 0.015, 10.0)
end

@testset "ParticleDA unit tests" begin
    dx = dy = 2e3

    x = collect(reshape(1.0:9.0, 3, 3, 1))
    # stations at (1,1) (2,2) and (3,3) return diagonal of x[3,3]
    ist = [1,2,3]
    jst = [1,2,3]
    obs = Vector{Float64}(undef, 3)
    Model.get_obs!(obs,x,3,ist,jst)
    @test obs ≈ [1.,5.,9.]

    y = [1.0, 2.0]
    cov_obs = float(I(2))
    weights = Vector{Float64}(undef, 3)
    # model observations with equal distance from true observation return equal weights
    hx = [0.5 0.9 1.5; 2.1 2.5 1.9]
    ParticleDA.get_log_weights!(weights, y, hx, cov_obs)
    @test weights ≈ [-1.9678770664093457, -1.9678770664093457, -1.9678770664093457]
    ParticleDA.normalized_exp!(weights)
    @test weights ≈ ones(3) / 3
    # model observations with decreasing distance from true observation return decreasing weights
    hx = [0.9 0.5 1.5; 2.1 2.5 3.5]
    ParticleDA.get_log_weights!(weights, y, hx, cov_obs)
    ParticleDA.normalized_exp!(weights)
    @test weights[1] > weights[2] > weights[3]

    # multivariate and independent methods give same weights when covariance matrix is diagonal
    weights2 = Vector{Float64}(undef, 3)
    ParticleDA.get_log_weights!(weights2, y, hx, 1.0)
    ParticleDA.normalized_exp!(weights2)
    @test weights2 ≈ weights

    id = zeros(Int, 5)
    # equal weights return the same particles
    w = ones(5) * .2
    ParticleDA.resample!(id,w)
    @test sort(id) == [1,2,3,4,5]
    # weight of 1.0 on first particle returns only copies of that particle
    w = zeros(5)
    w[1] = 1.0
    ParticleDA.resample!(id,w)
    @test id == [1,1,1,1,1]
    # weight of 1.0 on last particle returns only copies of that particle
    w = zeros(5)
    w[end] = 1.0
    ParticleDA.resample!(id,w)
    @test id == [5,5,5,5,5]
    # weights of .4 and .6 on particles 2 and 4 return a 40/60 mix of those particles
    w = zeros(5)
    w[2] = .4
    w[4] = .6
    ParticleDA.resample!(id,w)
    @test sort(id) == [2,2,4,4,4]

    nx = 10
    ny = 10
    dt = 1.0
    nt = 1
    # 0 input gives 0 output
    x0 = x = zeros(nx, ny, 3)
    model_matrices = LLW2d.setup(nx,ny,3.0e4, 0.1, 0.015, 10.0)
    @test size(model_matrices.absorbing_boundary) ==
        size(model_matrices.ocean_depth) ==
        size(model_matrices.x_averaged_depth) ==
        size(model_matrices.land_filter_m) ==
        size(model_matrices.land_filter_n) ==
        size(model_matrices.land_filter_e) == (nx, ny)
    dxeta = Matrix{Float64}(undef, nx, ny)
    dyeta = Matrix{Float64}(undef, nx, ny)
    Model.tsunami_update!(dxeta, dyeta, x, nt, dx, dy, dt, model_matrices)
    @test x ≈ x0

    # Initialise and update a tsunami on a small grid
    s = 4e3
    peak_position = [floor(Int, nx/4) * dx, floor(Int, ny/4) * dy]
    eta = @view(x[:,:,1])
    LLW2d.initheight!(eta, model_matrices, dx, dy, s, 1.0, peak_position)
    @test eta[2,2] ≈ 1.0
    @test sum(eta) ≈ 4.0
    Model.tsunami_update!(dxeta, dyeta, x, nt, dx, dy, dt, model_matrices)
    @test sum(eta, dims=1) ≈ [0.98161253125 1.8161253125 0.98161253125 0.073549875 0.0 0.0 0.0 0.0 0.0 0.0]
    @test sum(eta, dims=2) ≈ [0.98161253125 1.8161253125 0.98161253125 0.073549875 0.0 0.0 0.0 0.0 0.0 0.0]'

    # Test gaussian random field sampling
    x = 1.:2.
    y = 1.:2.
    sigma = lambda = nu = fill(1.0,3)
    grf = Model.init_gaussian_random_field_generator(lambda, nu, sigma,x,y,0,false)
    f = zeros(2, 2)
    rnn = [9.,9.,9.,9.]
    Model.sample_gaussian_random_field!(f,grf[1],rnn)
    @test f ≈ [16.2387054353321 5.115956753643808; 5.115956753643809 2.8210669567042155]

    # Test IO
    params_dict = YAML.load_file(joinpath(@__DIR__, "io_unit_test.yaml"))
    filter_params = ParticleDA.get_params(FilterParameters, params_dict["filter"])
    model_params = ParticleDA.get_params(ModelParameters, params_dict["model"]["llw2d"])
    rm(filter_params.output_filename, force=true)
    data1 = collect(reshape(1.0:(model_params.nx * model_params.ny), model_params.nx, model_params.ny, 1))
    data2 = randn(model_params.nx, model_params.ny, 1)
    tstep = 0
    h5open(filter_params.output_filename, "cw") do file
        Model.write_field(file, @view(data1[:,:,1]), tstep, "m", model_params.title_syn, "height", "unit test", model_params)
        Model.write_field(file, @view(data2[:,:,1]), tstep, "inch", model_params.title_avg,  "height", "unit test", model_params)
    end
    @test h5read(filter_params.output_filename, model_params.state_prefix * "_" * model_params.title_syn * "/t0000/height") ≈ data1
    @test h5read(filter_params.output_filename, model_params.state_prefix * "_" * model_params.title_avg * "/t0000/height") ≈ data2
    attr = h5readattr(filter_params.output_filename, model_params.state_prefix * "_" * model_params.title_syn * "/t0000/height")
    @test attr["Unit"] == "m"
    @test attr["Time step"] == tstep
    attr = h5readattr(filter_params.output_filename, model_params.state_prefix * "_" * model_params.title_avg * "/t0000/height")
    @test attr["Unit"] == "inch"
    @test attr["Time step"] == tstep
    Model.write_params(filter_params.output_filename, model_params)
    attr = h5readattr(filter_params.output_filename, model_params.title_params)
    @test attr["nx"] == model_params.nx
    @test attr["ny"] == model_params.ny
    @test attr["dx"] == model_params.dx
    @test attr["dy"] == model_params.dy
    @test attr["title_avg"] == model_params.title_avg
    @test attr["title_syn"] == model_params.title_syn
    Model.write_grid(filter_params.output_filename, model_params)
    attr = h5readattr(filter_params.output_filename, model_params.title_grid * "/x")
    @test attr["Unit"] == "m"
    attr = h5readattr(filter_params.output_filename, model_params.title_grid * "/y")
    @test attr["Unit"] == "m"

    stations = Model.StationVectors(zeros(Int,4), zeros(Int,4))
    Model.set_stations!(stations, model_params)
    @test stations.ist == [4, 4, 9, 9]
    @test stations.jst == [4, 9, 4, 9]
    Model.write_stations(filter_params.output_filename, stations.ist, stations.jst, model_params)
    @test h5read(filter_params.output_filename, model_params.title_stations * "/x") ≈ stations.ist .* model_params.dx
    @test h5read(filter_params.output_filename, model_params.title_stations * "/y") ≈ stations.jst .* model_params.dy
    attr = h5readattr(filter_params.output_filename, model_params.title_stations * "/x")
    @test attr["Unit"] == "m"
    attr = h5readattr(filter_params.output_filename, model_params.title_stations * "/y")
    @test attr["Unit"] == "m"

    rm(filter_params.output_filename, force=true)

    # Model data and get particles
    input_file = joinpath(dirname(pathof(ParticleDA)), "..", "test", "integration_test_1.yaml")
    model_params_dict = get(ParticleDA.read_input_file(input_file), "model", Dict())
    nprt_per_rank = 1
    my_rank = 0
    model_data = Model.init(model_params_dict, nprt_per_rank, my_rank, Random.TaskLocalRNG())
    # Make sure `get_particles` always returns the same array
    @test pointer(ParticleDA.get_particles(model_data)) == pointer(ParticleDA.get_particles(model_data))
end

@testset "ParticleDA integration tests" begin

    # Test true state with standard parameters
    x_true,x_avg,x_var = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_1.yaml"), BootstrapFilter())
    data_true = h5read(joinpath(@__DIR__, "reference_data.h5"), "integration_test_1")
    @test all(x_true .≈ data_true)

    # Test true state with different parameters
    x_true,x_avg,x_var = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_2.yaml"), BootstrapFilter())
    data_true = h5read(joinpath(@__DIR__, "reference_data.h5"), "integration_test_2")
    @test all(x_true .≈ data_true)

    # Test particle state with ~zero noise
    x_true,x_avg,x_var = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_3.yaml"), BootstrapFilter())
    @test x_true ≈ x_avg
    @test x_var .+ 1.0 ≈ ones(size(x_var))

    if Threads.nthreads() == 1

        avg_ref = h5read(joinpath(@__DIR__, "reference_data.h5"), "integration_test_4")
        # Test particle state with noise
        rng = StableRNG(123)
        x_true,x_avg,x_var = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_4.yaml"), BootstrapFilter(); rng)
        @test all(x_avg .≈ avg_ref)

        # Test that different seed gives different result
        rng = StableRNG(124)
        x_true,x_avg,x_var = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_4.yaml"), BootstrapFilter(); rng)
        @test !all(x_avg .≈ avg_ref)

        # Test optimal filter
        avg_ref = h5read(joinpath(@__DIR__, "reference_data.h5"), "integration_test_6")
        rng = StableRNG(123)
        x_true,x_avg,x_var = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "optimal_filter_test_1.yaml"), OptimalFilter(); rng)
        @test all(x_avg .≈ avg_ref)

    end

end

@testset "MPI -- $(file)" for file in ("mpi.jl", "copy_states.jl", "mean_and_var.jl")
    julia = joinpath(Sys.BINDIR, Base.julia_exename())
    flags = ["--startup-file=no", "-q", "-t$(Base.Threads.nthreads())"]
    script = joinpath(@__DIR__, file)
    mpiexec() do mpiexec
        mktempdir() do dir
            cd(dir) do
                @test success(run(ignorestatus(`$(mpiexec) -n 3 $(julia) $(flags) $(script)`)))
            end
        end
    end
end

@testset "Optimal Filter unit tests" begin

    seed = 123
    Random.seed!(seed)

    # Use default parameters other than smoothness parameter ν for state noise Matern
    # covariance function  which is fixed to 0.5 so that covariance function reduces to
    # r(x, y) = σ^2 exp(-norm([x, y]) / λ)
    model_params = ModelParameters(nu=[0.5, 0.5, 0.5])
    filter_params = FilterParameters()
    grid = ParticleDA.Grid(model_params.nx,
                           model_params.ny,
                           model_params.dx,
                           model_params.dy,
                           model_params.x_length,
                           model_params.y_length)
    grid_ext = ParticleDA.Grid((grid.nx-1)*2,
                               (grid.ny-1)*2,
                               grid.dx,
                               grid.dy,
                               (grid.x_length-grid.dx)*2,
                               (grid.y_length-grid.dy)*2)
    noise_params =  Matern(model_params.lambda[1], model_params.nu[1], σ = model_params.sigma[1])

    # Set station coordinates
    ist = rand(1:model_params.nx, model_params.nobs)
    jst = rand(1:model_params.ny, model_params.nobs)
    stations = (nst = model_params.nobs, ist = ist, jst = jst)
    cov_ext = ParticleDA.extended_covariance(0.0, 0.5 * grid.y_length, grid, noise_params)
    @test cov_ext ≈ noise_params.σ^2 * exp(-norm([0.0, 0.5 * model_params.y_length]) / noise_params.λ)
    @test cov_ext ≈ ParticleDA.extended_covariance(2.0 * grid.x_length, 0.5 * grid.y_length, grid, noise_params)
    @test cov_ext ≈ ParticleDA.extended_covariance(0.0, 1.5 * grid.y_length, grid, noise_params)
    arr = rand(ComplexF64,10,10)
    arr2 = zeros(ComplexF64,10,10)
    arr3 = zeros(ComplexF64,10,10)
    fft_plan, fft_plan! = FFTW.plan_fft(arr), FFTW.plan_fft!(arr)
    ParticleDA.normalized_2d_fft!(arr2, arr, fft_plan, fft_plan!, grid_ext)
    ParticleDA.normalized_2d_fft!(arr3, arr2, fft_plan, fft_plan!, grid_ext, inv)
    @test arr ≈ arr3

    cov_1 = zeros(stations.nst, grid_ext.nx * grid_ext.ny)
    cov_2 = zeros(grid.nx * grid.ny, stations.nst)
    cov_3 = zeros(stations.nst,stations.nst)
    ParticleDA.covariance_stations_extended_grid!(cov_1,grid,grid_ext,stations,noise_params)
    ParticleDA.covariance_grid_stations!(cov_2,grid,stations,noise_params)
    ParticleDA.covariance_stations!(cov_3,grid,stations,noise_params,model_params.obs_noise_std)
    @test all(isfinite, cov_1)
    @test all(isfinite, cov_2)
    @test all(isfinite, cov_3)
    @test cov_3 == Symmetric(cov_3)

    height = rand(grid.nx, grid.ny, filter_params.nprt)
    obs = randn(stations.nst)
    tmp_array = Matrix{ComplexF64}(undef, grid_ext.nx, grid_ext.ny)
    fft_plan, fft_plan! = FFTW.plan_fft(tmp_array), FFTW.plan_fft!(tmp_array)
    mat_off = ParticleDA.init_offline_matrices(grid, grid_ext, stations, noise_params, model_params.obs_noise_std, fft_plan, fft_plan!, Float64)
    mat_on = ParticleDA.init_online_matrices(grid, grid_ext, stations, filter_params.nprt, Float64)
    @test minimum(mat_off.Lambda) > 0.0
    ParticleDA.calculate_mean_height!(mat_on.mean, height, mat_off, obs, stations, grid, grid_ext, fft_plan, fft_plan!, filter_params.nprt, model_params.obs_noise_std)
    @test all(isfinite, mat_on.mean)

    rng = Random.seed!(seed)
    ParticleDA.sample_height_proposal!(height, mat_off, mat_on, obs, stations, grid, grid_ext, fft_plan, fft_plan!, filter_params.nprt, rng, model_params.obs_noise_std)
    @test all(isfinite, mat_on.samples)

end

@testset "Optimal Filter validation" begin

    include(joinpath(@__DIR__,"optimal_filter_validation.jl"));

    params_dict = YAML.load_file(joinpath(@__DIR__, "optimal_filter_validation.yaml"))
    filter_params = ParticleDA.get_params(FilterParameters, params_dict["filter"])
    model_params = ParticleDA.get_params(ModelParameters, params_dict["model"]["llw2d"])

    grid = ParticleDA.Grid(model_params.nx,
                           model_params.ny,
                           model_params.dx,
                           model_params.dy,
                           model_params.x_length,
                           model_params.y_length)
    grid_ext = ParticleDA.Grid((grid.nx-1)*2,
                               (grid.ny-1)*2,
                               grid.dx,
                               grid.dy,
                               (grid.x_length-grid.dx)*2,
                               (grid.y_length-grid.dy)*2)

    noise_params =  Matern(model_params.lambda[1], model_params.nu[1], σ = model_params.sigma[1])

    stations = (nst = model_params.nobs, ist = st.st_ij[:,1].+1, jst = st.st_ij[:,2].+1)

    h(x,y) = sin(x)^2 + cos(y)^2
    height = zeros(model_params.nx, model_params.ny, filter_params.nprt)
    x = (1:model_params.nx) .* 2 * pi / model_params.nx
    y = (1:model_params.ny) .* 4 * pi / model_params.ny
    for i = 1:filter_params.nprt
        height[:,:,i] .= h.(x',y)
    end

    obs = zeros(model_params.nobs)
    for i = 1:model_params.nobs
        obs[i] = height[stations.ist[i], stations.jst[i],1] + rand()
    end

    tmp_array = Matrix{ComplexF64}(undef, grid_ext.nx, grid_ext.ny)
    fft_plan, fft_plan! = FFTW.plan_fft(tmp_array), FFTW.plan_fft!(tmp_array)
    mat_off = ParticleDA.init_offline_matrices(grid, grid_ext, stations, noise_params, model_params.obs_noise_std, fft_plan, fft_plan!, Float64)
    mat_on = ParticleDA.init_online_matrices(grid, grid_ext, stations, filter_params.nprt, Float64)

    rng = Random.seed!(123)
    ParticleDA.sample_height_proposal!(height, mat_off, mat_on, obs, stations, grid, grid_ext, fft_plan, fft_plan!, filter_params.nprt, rng, model_params.obs_noise_std)

    Yobs_t = copy(obs)
    FH_t = copy(reshape(permutedims(height, [3 1 2]), filter_params.nprt, (model_params.nx)*(model_params.ny)))

    th = f_th(noise_params.σ[1], noise_params.λ[1], noise_params.ν[1])
    Mean_height = Calculate_Mean(FH_t, th, st, Yobs_t, Sobs, gr)
    Samples = Sample_Height_Proposal(FH_t, th, st, Yobs_t, Sobs, gr)

    @test mat_on.mean ≈ permutedims(reshape(Mean_height, filter_params.nprt, model_params.nx, model_params.ny), [2 3 1])
    @test mat_on.samples ≈ permutedims(reshape(Samples, filter_params.nprt, model_params.nx, model_params.ny), [2 3 1])

    # print("old mean height: ")
    # @btime Mean_height = Calculate_Mean($FH_t, $th, $st, $Yobs_t, $Sobs, $gr)
    # print("new mean height: ")
    # @btime calculate_mean_height!($mat_on.mean, $height, $mat_off, $obs, $stations, $params)
    # print("old sampling: ")
    # @btime Samples = Sample_Height_Proposal($FH_t, $th, $st, $Yobs_t, $Sobs, $gr)
    # print("new sampling: ")
    # @btime sample_height_proposal!($height, $mat_off, $mat_on, $obs, $stations, $params, $rng)

end


text:
"""
    LinearOperator{T<:VectorSpace,S<:VectorSpace,R<:Union{AbstractMatrix,Factorization}}

Compactly supported operator with effective domain and codomain.

Fields:
- `domain :: T`
- `codomain :: S`
- `coefficients :: R`
"""
struct LinearOperator{T<:VectorSpace,S<:VectorSpace,R<:Union{AbstractMatrix,Factorization}}
    domain :: T
    codomain :: S
    coefficients :: R
    function LinearOperator{T,S,R}(domain::T, codomain::S, coefficients::R) where {T<:VectorSpace,S<:VectorSpace,R<:Union{AbstractMatrix,Factorization}}
        dimension_domain = dimension(domain)
        dimension_codomain = dimension(codomain)
        sz₁, sz₂ = size(coefficients)
        (dimension_codomain == sz₁) & (dimension_domain == sz₂) || return throw(DimensionMismatch("dimensions must match: codomain and domain have dimensions $((dimension_codomain, dimension_domain)), coefficients has size $((sz₁, sz₂))"))
        return new{T,S,R}(domain, codomain, coefficients)
    end
end

LinearOperator(domain::T, codomain::S, coefficients::R) where {T<:VectorSpace,S<:VectorSpace,R<:Union{AbstractMatrix,Factorization}} =
    LinearOperator{T,S,R}(domain, codomain, coefficients)

domain(A::LinearOperator) = A.domain

codomain(A::LinearOperator) = A.codomain

coefficients(A::LinearOperator) = A.coefficients

# order, frequency

order(A::LinearOperator) = (order(domain(A)), order(codomain(A)))
order(A::LinearOperator, i::Int, j::Int) = (order(domain(A), j), order(codomain(A), i))

frequency(A::LinearOperator) = (frequency(domain(A)), frequency(codomain(A)))
frequency(A::LinearOperator, i::Int, j::Int) = (frequency(domain(A), j), frequency(codomain(A), i))

# utilities

function Base.firstindex(A::LinearOperator, i::Int)
    i == 1 && return _firstindex(codomain(A))
    i == 2 && return _firstindex(domain(A))
    return 1
end

function Base.lastindex(A::LinearOperator, i::Int)
    i == 1 && return _lastindex(codomain(A))
    i == 2 && return _lastindex(domain(A))
    return 1
end

Base.length(A::LinearOperator) = length(coefficients(A))

Base.size(A::LinearOperator) = size(coefficients(A))
Base.size(A::LinearOperator, i::Int) = size(coefficients(A), i)

Base.iterate(A::LinearOperator) = iterate(coefficients(A))
Base.iterate(A::LinearOperator, i) = iterate(coefficients(A), i)

Base.eltype(A::LinearOperator) = eltype(coefficients(A))
Base.eltype(::Type{LinearOperator{T,S,R}}) where {T,S,R} = eltype(R)

# getindex, view, setindex!

Base.@propagate_inbounds Base.getindex(A::LinearOperator, α, β) =
    getindex(coefficients(A), _findposition(α, codomain(A)), _findposition(β, domain(A)))

Base.@propagate_inbounds Base.view(A::LinearOperator, α, β) =
    view(coefficients(A), _findposition(α, codomain(A)), _findposition(β, domain(A)))

Base.@propagate_inbounds Base.setindex!(A::LinearOperator, x, α, β) =
    setindex!(coefficients(A), x, _findposition(α, codomain(A)), _findposition(β, domain(A)))

# ==, iszero, isapprox

Base.:(==)(A::LinearOperator, B::LinearOperator) =
    codomain(A) == codomain(B) && domain(A) == domain(B) && coefficients(A) == coefficients(B)

Base.iszero(A::LinearOperator) = iszero(coefficients(A))

Base.isapprox(A::LinearOperator, B::LinearOperator; kwargs...) =
    codomain(A) == codomain(B) && domain(A) == domain(B) && isapprox(coefficients(A), coefficients(B); kwargs...)

# copy, similar

Base.copy(A::LinearOperator) = LinearOperator(domain(A), codomain(A), copy(coefficients(A)))

Base.similar(A::LinearOperator) = LinearOperator(domain(A), codomain(A), similar(coefficients(A)))
Base.similar(A::LinearOperator, ::Type{T}) where {T} = LinearOperator(domain(A), codomain(A), similar(coefficients(A), T))

# zero, one

function Base.zero(A::LinearOperator)
    domain_A = domain(A)
    codomain_A = codomain(A)
    CoefType = eltype(A)
    C = LinearOperator(domain_A, codomain_A, Matrix{CoefType}(undef, dimension(codomain_A), dimension(domain_A)))
    coefficients(C) .= zero(CoefType)
    return C
end

function Base.one(A::LinearOperator)
    CoefType = eltype(A)
    C = zero(A)
    @inbounds for α ∈ indices(domain(A) ∩ codomain(A))
        C[α,α] = one(CoefType)
    end
    return C
end

# float, complex, real, imag, conj, conj!

for f ∈ (:float, :complex, :real, :imag, :conj, :conj!)
    @eval Base.$f(A::LinearOperator) = LinearOperator(domain(A), codomain(A), $f(coefficients(A)))
end

# promotion

Base.convert(::Type{T}, A::T) where {T<:LinearOperator} = A
Base.convert(::Type{LinearOperator{T₁,S₁,R₁}}, A::LinearOperator{T₂,S₂,R₂}) where {T₁,S₁,R₁,T₂,S₂,R₂} =
    LinearOperator{T₁,S₁,R₁}(convert(T₁, domain(A)), convert(S₁, codomain(A)), convert(R₁, coefficients(A)))

Base.promote_rule(::Type{T}, ::Type{T}) where {T<:LinearOperator} = T
Base.promote_rule(::Type{LinearOperator{T₁,S₁,R₁}}, ::Type{LinearOperator{T₂,S₂,R₂}}) where {T₁,S₁,R₁,T₂,S₂,R₂} =
    LinearOperator{promote_type(T₁, T₂), promote_type(S₁, S₂), promote_type(R₁, R₂)}

# Cartesian spaces

Base.eachcol(A::LinearOperator{<:CartesianSpace,<:CartesianSpace}) =
    (@inbounds(component(A, :, j)) for j ∈ Base.OneTo(nb_cartesian_product(domain(A))))
Base.eachrow(A::LinearOperator{<:CartesianSpace,<:CartesianSpace}) =
    (@inbounds(component(A, i, :)) for i ∈ Base.OneTo(nb_cartesian_product(codomain(A))))

eachcomponent(A::LinearOperator{<:CartesianSpace,<:CartesianSpace}) =
    (@inbounds(component(A, i, j)) for i ∈ Base.OneTo(nb_cartesian_product(codomain(A))), j ∈ Base.OneTo(nb_cartesian_product(domain(A))))
eachcomponent(A::LinearOperator{<:CartesianSpace,<:VectorSpace}) =
    (@inbounds(component(A, j)) for i ∈ Base.OneTo(1), j ∈ Base.OneTo(nb_cartesian_product(domain(A))))
eachcomponent(A::LinearOperator{<:VectorSpace,<:CartesianSpace}) =
    (@inbounds(component(A, i)) for i ∈ Base.OneTo(nb_cartesian_product(codomain(A))), j ∈ Base.OneTo(1))

Base.@propagate_inbounds component(A::LinearOperator{<:CartesianSpace,<:CartesianSpace}, i, j) =
    LinearOperator(domain(A)[j], codomain(A)[i], view(coefficients(A), _component_findposition(i, codomain(A)), _component_findposition(j, domain(A))))

Base.@propagate_inbounds component(A::LinearOperator{<:CartesianSpace,<:VectorSpace}, j) =
    LinearOperator(domain(A)[j], codomain(A), view(coefficients(A), :, _component_findposition(j, domain(A))))

Base.@propagate_inbounds component(A::LinearOperator{<:VectorSpace,<:CartesianSpace}, i) =
    LinearOperator(domain(A), codomain(A)[i], view(coefficients(A), _component_findposition(i, codomain(A)), :))

# show

function Base.show(io::IO, ::MIME"text/plain", A::LinearOperator)
    println(io,
        "LinearOperator : " * string_space(domain(A)) * " -> " * string_space(codomain(A)) *
        " with coefficients " * string(typeof(coefficients(A))) * ":"
        )
    Base.print_array(io, coefficients(A))
end


text:
using SimpleTensorNetworks
using Test

@testset "contract" begin
    include("tensorcontract.jl")
end

@testset "simplify" begin
    include("simplify.jl")
end

@testset "contractionorder" begin
    include("contractionorder/contractionorder.jl")
end

@testset "viz" begin
    include("viz.jl")
end

@testset "cuda" begin
    if Base.find_package("CUDA") !== nothing
        include("cuda.jl")
    end
end


text:
"""
Uses `Plots.jl` to plot propagation results for the n-body problem,
and the two-body problem. Uses `Plots.jl`
"""
module OrbitPlots

using Reexport
@reexport using ..AstrodynamicsCore

include("../Misc/DocStringExtensions.jl")
include("../Misc/UnitfulAliases.jl")

using ..TwoBody
using ..ThreeBody
using ..NBody
using ..Propagators

using Plots
using Plots.PlotMeasures

export orbitplot, lagrangeplot

include("PlotTwoBody.jl")
include("PlotNBody.jl")
include("PlotThreeBody.jl")

end


text:
immutable Cauchy <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function Cauchy(l::Real, s::Real)
        s > zero(s) || error("scale must be positive")
	new(float64(l), float64(s))
    end
end

Cauchy(l::Real) = Cauchy(l, 1.0)
Cauchy() = Cauchy(0.0, 1.0)

## Support
@continuous_distr_support Cauchy -Inf Inf

## Properties
mean(d::Cauchy) = NaN
median(d::Cauchy) = d.location
mode(d::Cauchy) = d.location

var(d::Cauchy) = NaN
skewness(d::Cauchy) = NaN
kurtosis(d::Cauchy) = NaN

entropy(d::Cauchy) = log(d.scale) + log(4.0 * pi)

## Functions
pdf(d::Cauchy, x::Real) = 1/(pi*d.scale*(1+((x-d.location)/d.scale)^2))
logpdf(d::Cauchy, x::Real) = -log(pi) - log(d.scale) - log1psq((x-d.location)/d.scale)

cdf(d::Cauchy, x::Real) = atan2(one(x),-(x-d.location)/d.scale)/pi
ccdf(d::Cauchy, x::Real) = atan2(one(x),(x-d.location)/d.scale)/pi
logcdf(d::Cauchy, x::Real) = x <= d.location ? log(cdf(d,x)) : log1p(-ccdf(d,x))
logccdf(d::Cauchy, x::Real) = x <= d.location ? log1p(-cdf(d,x)) : log(ccdf(d,x))

quantile(d::Cauchy, p::Real) = @checkquantile p d.location - d.scale*cospi(p)/sinpi(p)
cquantile(d::Cauchy, p::Real) = @checkquantile p d.location + d.scale*cospi(p)/sinpi(p)
invlogcdf(d::Cauchy, lp::Real) = lp < loghalf ? quantile(d,exp(lp)) : cquantile(d,-expm1(lp))
invlogccdf(d::Cauchy, lp::Real) = lp < loghalf ? cquantile(d,exp(lp)) : quantile(d,-expm1(lp))

mgf(d::Cauchy, t::Real) = t == zero(t) ? 1.0 : NaN
cf(d::Cauchy, t::Real) = exp(im * t * d.location - d.scale * abs(t))

## Fitting
# Note: this is not a Maximum Likelihood estimator
function fit{T <: Real}(::Type{Cauchy}, x::Array{T})
    l, u = iqr(x)
    Cauchy(median(x), (u - l) / 2.0)
end


text:
testobject(val) = Dict("tree" => Dict(
    "val" => val
))

@info "deleteenvkey!"
@test !haskey(ENV, "FOOTEST")
@test_nowarn Configs.deleteenvkey!("FOOTEST")
ENV["FOOTEST"] = "foo"
@test_nowarn Configs.deleteenvkey!("FOOTEST")
@test !haskey(ENV, "FOOTEST")

@info "parseenvkey"
@test Configs.parseenvkey("FOOTEST", "nope") === "nope"
ENV["FOOTEST"] = "yep"
@test Configs.parseenvkey("FOOTEST", "nope") === "yep"
Configs.deleteenvkey!("FOOTEST")

@info "parsecustomenv!"
foo = testobject("FOOTEST")
@test_nowarn Configs.parsecustomenv!(foo)
@test isequal(foo, Dict("tree" => Dict()))
foo = testobject("FOOTEST")
ENV["FOOTEST"] = "yep"
@test_nowarn Configs.parsecustomenv!(foo)
@test isequal(foo, testobject("yep"))
Configs.deleteenvkey!("FOOTEST")

@info "tosymbol"
testdict = @test_nowarn Configs.tosymbol(Dict{String, Any}("foo" => "FOO", "bar" => "BAR"))
for key in keys(testdict)
   @test key isa Symbol
end
@test testdict[Symbol("foo")] === "FOO"
@test Configs.tosymbol("foo") === "foo"

@info "makeimmutable"
isimmutable = @test_nowarn Configs.makeimmutable(testobject(["one", Dict("one" => 1)]))
@test isequal(isimmutable, (tree = (val = ("one", (one = 1,)),),))

@info "override!"
base = Dict()
add = testobject("test")
@test_nowarn Configs.override!(base, add)
@test isequal(base, add)
add = Dict("foo" => ["one", "two"], "bar" => 2, "tree" => Dict("added" => 2))
@test_nowarn Configs.override!(base, add)
@test isequal(base, Dict{Any,Any}("bar" => 2,"tree" => Dict{Any,Any}("val" => "test","added" => 2),"foo" => ["one", "two"]))

@info "getfiles"
@test_throws Configs.Configserror Configs.getfiles("foobar")
@test_nowarn Configs.getfiles("configs")
@test_nowarn Configs.getfiles(joinpath(pwd(), "configs"))

max_stars_repo_name: tisztamo/Configs

max_issues_repo_name: tisztamo/Configs

max_forks_repo_name: tisztamo/Configs


text:
# This runs only on Travis. The goal is to populate the `.julia/compiled/v*` directory
# with some additional files, so that `filter_valid_cachefiles` has to run.
# This is to catch problems like #460.

using Pkg
using Test

Pkg.add(PackageSpec(name="EponymTuples", version="0.2.0"))
using EponymTuples  # force compilation
id = Base.PkgId(EponymTuples)
paths = Base.find_all_in_cache_path(id)
Pkg.rm("EponymTuples") # we don't need it anymore
path = first(paths)
base, ext = splitext(path)
mv(path, base*"blahblah"*ext)
Pkg.add(PackageSpec(name="EponymTuples"))


text:
using Knet,JLD

function main()
    loadmodel()
    loadvocab()
    loadcorpus()
    global data1 = minibatch(corpus, 1)
    global data100 = minibatch(corpus, 100)
    perplexity(model,data100)
end

function loadmodel()
    global model = Dict()
    m = load("turkish_20kmodel.jld")["model"]
    for k in (:fembed,:bembed); model[k] = KnetArray(m[k]); end
    for k in (:forw,:soft,:back); model[k] = map(KnetArray, m[k]); end
end

function loadvocab()
    global vocab = load("turkish_20kvocab.jld")["vocab"]
end

function loadcorpus()
    global corpus = Vector{Int32}[]
    sentence = Int32[]
    unk = vocab["<unk>"]
    for line in eachline("tr-ud-train.conllu")
#    for line in eachline("okirnap/turkish_dev.conllu")
        line = chomp(line)
        fields = split(line,'\t')
        if line == ""
            if !isempty(sentence)
                push!(corpus, sentence)
                sentence = Int32[]
            end
        elseif ismatch(r"^\d+$", fields[1])
            push!(sentence, get(vocab, fields[2], unk))
        elseif line[1] == '#'
            # skip
        elseif ismatch(r"^\d+[-.]\d+$", fields[1])
            # skip
        else
            error("Cannot parse [$line]")
        end
    end
end

function minibatch(corpus, batchsize)
    function batch(a)
        B = length(a)
        T = length(corpus[a[1]])
        s = [ zeros(Int32,B) for i=1:T ]
        for t=1:T, b=1:B
            s[t][b] = corpus[a[b]][t]
        end
        return s
    end
    data = Any[]
    d = Dict{Int,Vector{Int}}()
    for i in 1:length(corpus)
        s = corpus[i]
        l = length(s)
        a = get!(d, l, Int[])
        push!(a, i)
        if length(a) == batchsize
            push!(data, batch(a))
            empty!(a)
        end
    end
    for (k,a) in d
        if !isempty(a)
            push!(data, batch(a))
        end
    end
    return data
end

function blstm(model, sequence; result=nothing)
    # assume the sequence does not have <s> or </s> or padding and add it
    sos = sosinput(model, sequence)
    eos = eosinput(model, sequence)
    sequence = [ sos, sequence..., eos ]
    # @show sequence
    T = length(sequence)
    # sequence[t] is a Vector{Int32} minibatch of B token ids
    B = length(sequence[1])
    # forw[t] and back[t] are the hiddens used to predict sequence[t-1]
    forw = Array(Any,T)
    back = Array(Any,T)
    # zero hidden for first time step
    h0 = zerohidden(model, sequence)
    # run the forward lstm
    hidden = cell = forw[1] = h0
    weight,bias = model[:forw]
    embed = model[:fembed]
    for t=1:T-1
        hidden,cell = lstm(weight,bias,hidden,cell,embed[sequence[t],:])
        forw[t+1] = hidden
    end
    # run the backward lstm
    hidden = cell = back[T] = h0
    weight,bias = model[:back]
    embed = model[:bembed]
    for t=T:-1:2
        hidden,cell = lstm(weight,bias,hidden,cell,embed[sequence[t],:])
        back[t-1] = hidden
    end
    # run the predictions
    total = 0
    for t=1:T
        ypred = hcat(forw[t], back[t]) * model[:soft][1] .+ model[:soft][2]
        total += logprob(sequence[t], ypred)
        # @show tmp = logprob(sequence[t], ypred)
        # total += tmp
    end
    if result != nothing
        result[1] = total
        result[2] = B*T
    end
    # return -total 	    # total loss: longer sequences and larger minibatches have higher loss
    # return -total / (B*T) # per token loss: scale does not depend on sequence length or minibatch
    return -total / B       # per sequence loss: does not depend on minibatch, larger loss for longer seq
end

function blstm1(model, sequence; result=nothing)
    # assume the sequence does not have <s> or </s> or padding
    T = length(sequence)
    # sequence[t] is a Vector{Int32} minibatch of B token ids
    B = length(sequence[1])
    # forw[t] and back[t] are the hiddens used to predict sequence[t]
    forw = Array(Any,T)
    back = Array(Any,T)
    # zero hidden and sos token minibatch for first time step
    h0 = zerohidden(model, sequence)
    x0 = sosinput(model, sequence)
    # run the forward lstm
    hidden,cell,input = h0,h0,x0
    weight,bias = model[:forw]
    embed = model[:fembed]
    for t=1:T
        hidden,cell = lstm(weight,bias,hidden,cell,embed[input,:])
        forw[t] = hidden
        input = sequence[t]
    end
    # run the backward lstm
    hidden,cell,input = h0,h0,x0
    weight,bias = model[:back]
    embed = model[:bembed]
    for t=T:-1:1
        hidden,cell = lstm(weight,bias,hidden,cell,embed[input,:])
        back[t] = hidden
        input = sequence[t]
    end
    # run the predictions
    total = 0
    for t=1:T
        ypred = hcat(forw[t], back[t]) * model[:soft][1] .+ model[:soft][2]
        # total += logprob(sequence[t], ypred)
        @show tmp = logprob(sequence[t], ypred)
        total += tmp
    end
    if result != nothing
        result[1] = total
        result[2] = B*T
    end
    # return -total 	    # total loss: longer sequences and larger minibatches have higher loss
    # return -total / (B*T) # per token loss: scale does not depend on sequence length or minibatch
    return -total / B       # per sequence loss: does not depend on minibatch, larger loss for longer seq
end

function perplexity(model, data)
    result = zeros(2)
    total = zeros(2)
    for sequence in data
        blstm(model, sequence; result=result)
        total += result
    end
    @show total
    return exp(-total[1]/total[2])
end

# row-major lstm
function lstm(weight,bias,hidden,cell,input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

let h0=sos=eos=nothing; global zerohidden, sosinput, eosinput
    function zerohidden(model, sequence)
        B = length(sequence[1])
        W = model[:forw][1]
        H = div(size(W,2),4)
        if h0 == nothing || size(h0) != (B,H) || typeof(h0) != typeof(W)
            h0 = fill!(similar(W, B, H), 0)
        end
        return h0
    end
    function sosinput(model, sequence)
        if sos == nothing || length(sos) != length(sequence[1])
            sos = fill!(similar(sequence[1]), vocab["<s>"])
        end
        return sos
    end
    function eosinput(model, sequence)
        if eos == nothing || length(eos) != length(sequence[1])
            eos = fill!(similar(sequence[1]), vocab["</s>"])
        end
        return eos
    end
end

function logprob(output, ypred)
    nrows,ncols = size(ypred)
    index = similar(output)
    @inbounds for i=1:length(output)
        index[i] = i + (output[i]-1)*nrows
    end
    o1 = logp(ypred,2)
    o2 = o1[index]
    o3 = sum(o2)
    return o3
end


text:
using YaoQDNN
using Yao
using YaoAD
using Flux
using Test

using Zygote

@testset "layer.jl" begin
    ql_1 = QNNL{Float64}(hardware_efficient, 8, 12, 20, "YZ")
    ql_2 = QNNL{Float64}(hardware_efficient, 4, 16, 10, "X")
    m = Chain(ql_1, ql_2)
    loss(x, y) = Flux.mse(m(x), y)
    x = ones(12)
    y = ones(4)

    gs = gradient(()->loss(x,y), params(ql_2))

    _, back = Zygote.pullback(Flux.mse, m(x), y)
    m̄ = back(1)[1]

    @test all(gs[ql_2.params] .≈ YaoQDNN.forwardAD(ql_2, ql_1(x), m̄)[1][1])
    @test all(gs[ql_2.bias] .≈ YaoQDNN.forwardAD(ql_2, ql_1(x), m̄)[1][2])

    # Float32
    ql_1 = QNNL{Float32}(hardware_efficient, 8, 12, 20, "YZ")
    ql_2 = QNNL{Float32}(hardware_efficient, 4, 16, 10, "X")
    m = Chain(ql_1, ql_2)
    loss(x, y) = Flux.mse(m(x), y)
    x = ones(Float32, 12)
    y = ones(Float32, 4)

    _, back = Zygote.pullback(Flux.mse, m(x), y)
    m̄ = back(1)[1]

    gs = gradient(()->loss(x,y), params(ql_2))
    @test all(gs[ql_2.params] .≈ YaoQDNN.forwardAD(ql_2, ql_1(x), m̄)[1][1])
    @test all(gs[ql_2.bias] .≈ YaoQDNN.forwardAD(ql_2, ql_1(x), m̄)[1][2])
end


text:
function I_kron_A(d::Int, A::AbstractMatrix{T}) where {T}
    n, m = size(A)
    B = zeros(T,d*n,d*m)
    @inbounds for k = 1:d
        i = (k-1)*n + 1
        j = (k-1)*m + 1
        B[i:i+n-1, j:j+m-1] .= A
    end
    return B
end

function I_kron_mul!(C,A,B)
    n = checksquare(A)
    n^2 == checksquare(B) == checksquare(C) || throw(DimensionMismatch("invalid dimensions"))
    for j = 1:n
        jj = (j-1)*n+1 : (j-1)*n+n
        for i = 1:n
            ii = (i-1)*n+1 : (i-1)*n+n
            Cij = view(C,ii,jj)
            Bij = view(B,ii,jj)
            mul!(Cij,A,Bij)
        end
    end
    return C
end

# Not too useful:
function At_kron_I_mul_B!(C,A,B)
    n = checksquare(A)
    n^2 == checksquare(B) == checksquare(C) || throw(DimensionMismatch("invalid dimensions"))
    fill!(C,0)
    kk = 1:n
    for j = 1:n
        jj = (j-1)*n+1 : (j-1)*n+n
        for i = 1:n
            ii = (i-1)*n+1 : (i-1)*n+n
            Cij = view(C,ii,jj)
            for k = 1:n
                Cij .+= A[k,i] .* view(B,kk+(k-1)*n,jj)
            end
        end
    end
    return C
end


text:
function prob_func(prob,i,repeat)
    x = 0.3rand(2)
    # remake(prob,p=[p[1:2];x])
    remake(prob)
end

ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sim = solve(ensemble_prob,SRIW1(),trajectories=10^2)
summ = EnsembleSummary(sim,0:0.1:10)


text:
module JotTest1

using Format

function response_func(d::Dict)
  if haskey(d, "add suffix")
    s = d["add suffix"] |> String
    package_root = Base.moduleroot(JotTest1) |> pathof |> splitpath
    rs_path = joinpath(package_root[begin:end-2]..., "response_suffix")
    open(rs_path, "r") do rs
      response_suffix = readchomp(rs)
      return s * response_suffix
    end
  elseif haskey(d, "double")
    n = d["double"] |> Float64
    return n * 2
  else
    error("Input not recognized")
  end
end

end # module


text:
module SimpleSATSolver

export SimpleSAT, iterative_solve

import SatisfiabilityInterface: AbstractSATSolver, SATProblem, solve

struct SimpleSAT <: AbstractSATSolver
end

include("structured_sat_problem.jl")
include("solve.jl")

# SimpleSAT() = SimpleSAT(false)

# TODO: When a new variable is assigned, check that no clauses are violated

# TODO: Unit propagation

end

max_stars_repo_name: dpsanders/SimpleSATSolver

max_issues_repo_name: dpsanders/SimpleSATSolver

max_forks_repo_name: dpsanders/SimpleSATSolver


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "HIP"
version = v"4.0.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/ROCm-Developer-Tools/HIP/archive/rocm-$(version).tar.gz",
                  #"e21c10b62868ece7aa3c8413ec0921245612d16d86d81fe61797bf9a64bc37eb"), # 4.1.0
                  "d7b78d96cec67c55b74ea3811ce861b16d300410bc687d0629e82392e8d7c857"), # 4.0.0
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/HIP*/

# disable broken Pre-Compiled Headers (PCH)
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/disable-pch.patch"

# help hipcc to find BB C++ headers+libs+crt
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/hipcc-bb-paths.patch"

# use rocclr instead of hcc as HIP_PLATFORM
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/rocclr-not-hcc-platform.patch"
#sed -e "s:\$HIP_PLATFORM eq \"hcc\" and \$HIP_COMPILER eq \"clang\":\$HIP_PLATFORM eq \"rocclr\" and \$HIP_COMPILER eq \"clang\":" -i bin/hipcc

# force hipcc to pass paths to clang
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/hipcc-force-paths.patch"

apk add coreutils dateutils

mkdir build && cd build
export HIP_CLANG_PATH=${prefix}/tools
cmake -DCMAKE_INSTALL_PREFIX=${prefix}/hip \
      -DCMAKE_PREFIX_PATH=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DHIP_COMPILER=clang \
      -DHIP_PLATFORM=rocclr \
      -DHIP_RUNTIME=ROCclr \
      -DROCM_PATH=${prefix} \
      -DHSA_PATH=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_clang.cmake \
      -DLLVM_DIR="${libdir}/cmake/llvm" \
      -DClang_DIR="${libdir}/cmake/clang" \
      -DLLD_DIR="${libdir}/cmake/ldd" \
      ..
make -j${nproc}
make install

# link to .hipVersion so clang can find it
ln -s ../hip/bin/.hipVersion ${bindir}/.hipVersion
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx11"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct(["libamdhip64"], :libamdhip64, ["hip/lib"]),
    ExecutableProduct("hipcc", :hipcc, "hip/bin"),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(; name="LLVM_full_jll", version=v"11.0.1")),
    Dependency("hsa_rocr_jll"),
    Dependency("ROCmDeviceLibs_jll"),
    Dependency("ROCmCompilerSupport_jll"),
    Dependency("ROCmOpenCLRuntime_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               preferred_gcc_version=v"8")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil