text:
function Base.promote_array_type{FSA <: FixedArray, N, T<:Number}(
        F, ::Type{T}, ::Type{Array{FSA, N}}
    )
    FSA
end
function Base.promote_array_type{FSA <: FixedArray, N, T<:Number}(
        F, ::Type{Array{FSA, N}}, ::Type{T}
    )
    FSA
end
function Base.promote_array_type{FSA <: FixedArray, T<:Number}(
        F, ::Type{T}, ::Type{FSA}
    )
    FSA
end
function Base.promote_array_type{FSA <: FixedArray, T<:Number}(
        F, ::Type{FSA}, ::Type{T}
    )
    FSA
end

# operations
const unaryOps = (:-, :~, :conj, :abs,
                  :sin, :cos, :tan, :sinh, :cosh, :tanh,
                  :asin, :acos, :atan, :asinh, :acosh, :atanh,
                  :sec, :csc, :cot, :asec, :acsc, :acot,
                  :sech, :csch, :coth, :asech, :acsch, :acoth,
                  :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
                  :sind, :tand, :acosd, :acotd, :acscd, :asecd,
                  :asind, :atand, :rad2deg, :deg2rad,
                  :log, :log2, :log10, :log1p, :exponent, :exp,
                  :exp2, :expm1, :cbrt, :sqrt, :erf,
                  :erfc, :erfcx, :erfi, :dawson, :ceil, :floor,
                  :trunc, :round, :significand, :lgamma,
                  :gamma, :lfact, :frexp, :modf, :airy, :airyai,
                  :airyprime, :airyaiprime, :airybi, :airybiprime,
                  :besselj0, :besselj1, :bessely0, :bessely1,
                  :eta, :zeta, :digamma, :real, :imag)

# vec-vec and vec-scalar
const binaryOps = (:.+, :.-,:.*, :./, :.\, :.^,
                   :.==, :.!=, :.<, :.<=, :.>, :.>=, :+, :-,
                   :min, :max,
                   :div, :fld, :rem, :mod, :mod1, :cmp,
                   :atan2, :besselj, :bessely, :hankelh1, :hankelh2,
                   :besseli, :besselk, :beta, :lbeta)

const reductions = ((:sum,:+), (:prod,:*), (:minimum,:min), (:maximum,:max))

function gen_functor(func::Symbol, unary::Int)
    functor_name  = gensym()
    arguments     = ntuple(i->Symbol("arg$i"), unary)
    functor_expr  = quote
        immutable $functor_name <: Functor{$unary} end
        @compat @inline (::$functor_name)($(arguments...)) = $func($(arguments...))
    end
    return (functor_name, functor_expr)
end

for (callfun, reducefun) in reductions
    functor_name, functor_expr = gen_functor(reducefun, 2)
    eval(quote
        $functor_expr
        @inline $(callfun){T <: FixedArray}(x::T) = reduce($functor_name(), x)
    end)
end

for op in unaryOps
    functor_name, functor_expr = gen_functor(op, 1)
    eval(quote
        $functor_expr
        @inline $(op){T <: FixedArray}(x::T) = map($functor_name(), x)
    end)
end

for op in binaryOps
    functor_name, functor_expr = gen_functor(op, 2)
    eval(quote
        $functor_expr
        @inline $op{F1<:FixedArray, F2<:FixedArray}(x::F1, y::F2) = map($functor_name(), x, y)
        @inline $op{F<:FixedArray}(x::F, y::Number) = map($functor_name(), x, y)
        @inline $op{F<:FixedArray}(x::Number, y::F) = map($functor_name(), x, y)
        @inline $op{T,N}(x::FixedArray{T,N}, y::AbstractArray{T,N}) = map($functor_name(), x, y)
        @inline $op{T,N}(x::AbstractArray{T,N}, y::FixedArray{T,N}) = map($functor_name(), x, y)
    end)
end

# Non-dotted * and / of FSAs and scalars
mul_func_name, mul_func_expr = gen_functor(:*, 2)
div_func_name, div_func_expr = gen_functor(:/, 2)
@eval begin
    $div_func_expr
    $mul_func_expr
    @inline *(x::Number, y::FixedArray) = map($mul_func_name(), x, y)
    @inline *(x::FixedArray, y::Number) = map($mul_func_name(), x, y)
    @inline /(x::FixedArray, y::Number) = map($div_func_name(), x, y)
    function /(x::Number, y::FixedArray)
        Base.depwarn("/(x::Number, y::FixedArray) had unclear semantics and will be removed", :/)
        map($div_func_name(), x, y)
    end
end

@inline function promote{T1 <: FixedArray, T2 <: FixedArray}(a::T1, b::T2)
    T = promote_type(eltype(T1), eltype(T2))
    map(T, a), map(T, b)
end
function promote{T1, T2 <: Number}(a::FixedArray{T1}, b::T2)
    T = promote_type(T1, T2)
    map(T, a), T(b)
end
function promote{T1 <: Number, T2}(a::T1, b::FixedArray{T2})
    T = promote_type(T1, T2)
    T(a), map(T, b)
end

function promote_rule{N, T, X<:Number}(::Type{Vec{N,T}}, ::Type{X})
    Vec{N, promote_type(T, X)}
end

@inline ctranspose{R, C, T}(a::Mat{R, C, T}) = Mat{C,R,T}(ntuple(CRowFunctor(a), Val{R}))
@generated function ctranspose{N,T}(b::Vec{N,T})
    expr = [:(Tuple(b)[$i]',) for i=1:N]
    return quote
        Mat{1,N,T}($(expr...))
    end
end
@inline transpose{R, C, T}(a::Mat{R, C, T}) = Mat(ntuple(RowFunctor(a), Val{R}))
@generated function transpose{N,T}(b::Vec{N,T})
    expr = [:(transpose(Tuple(b)[$i]),) for i=1:N]
    return quote
        Mat{1,N,T}($(expr...))
    end
end
@inline Base.hypot{T}(v::FixedVector{2,T}) = hypot(v[1],v[2])

immutable DotFunctor <: Functor{2} end
@compat (::DotFunctor)(a, b) = a'*b
@inline dot(a::FixedVector, b::FixedVector) = sum(map(DotFunctor(), a, b))

immutable BilinearDotFunctor <: Functor{2} end
@compat (::BilinearDotFunctor)(a, b) = a*b
@inline bilindot(a::FixedVector, b::FixedVector) = sum(map(BilinearDotFunctor(), a, b))

function cross{T}(a::FixedVector{2, T}, b::FixedVector{2, T})
    return a[1]*b[2]-a[2]*b[1]
end
@inline function cross{T1, T2}(a::FixedVector{3, T1}, b::FixedVector{3, T2})
    @inbounds elements = (a[2]*b[3]-a[3]*b[2],
                          a[3]*b[1]-a[1]*b[3],
                          a[1]*b[2]-a[2]*b[1])
    construct_similar(typeof(a), elements)
end

@inline norm{N, T}(a::FixedVector{N, T})     = sqrt(dot(a,a))
@inline normalize{FSA <: FixedArray}(a::FSA) = a / norm(a)

function norm{N, T}(a::FixedVector{N, T}, p)
    isinf(p) && return maxabs(a)
    ret = abs(a[1])^p
    for k = 2:N
        ret += abs(a[k])^p
    end
    ret^(1/p)
end

function Base.isnan(p::FixedArray)
    for elem in p
        isnan(elem) && return true
    end
    false
end

#Matrix
@inline det{T}(A::FixedMatrix{1, 1, T}) = @inbounds return ( A[1] )
@inline det{T}(A::FixedMatrix{2, 2, T}) = @inbounds return ( A[1,1]*A[2,2] - A[1,2]*A[2,1])
@inline det{T}(A::FixedMatrix{3, 3, T}) = @inbounds return (
    A[1,1]*(A[2,2]*A[3,3]-A[2,3]*A[3,2]) -
    A[1,2]*(A[2,1]*A[3,3]-A[2,3]*A[3,1]) +
    A[1,3]*(A[2,1]*A[3,2]-A[2,2]*A[3,1])
)
@inline det{T}(A::FixedMatrix{4, 4, T}) = @inbounds return (
    A[13] * A[10]  * A[7]  * A[4]  - A[9] * A[14] * A[7]  * A[4]   -
    A[13] * A[6]   * A[11] * A[4]  + A[5] * A[14] * A[11] * A[4]   +
    A[9]  * A[6]   * A[15] * A[4]  - A[5] * A[10] * A[15] * A[4]   -
    A[13] * A[10]  * A[3]  * A[8]  + A[9] * A[14] * A[3]  * A[8]   +
    A[13] * A[2]   * A[11] * A[8]  - A[1] * A[14] * A[11] * A[8]   -
    A[9]  * A[2]   * A[15] * A[8]  + A[1] * A[10] * A[15] * A[8]   +
    A[13] * A[6]   * A[3]  * A[12] - A[5] * A[14] * A[3]  * A[12]  -
    A[13] * A[2]   * A[7]  * A[12] + A[1] * A[14] * A[7]  * A[12]  +
    A[5]  * A[2]   * A[15] * A[12] - A[1] * A[6]  * A[15] * A[12]  -
    A[9]  * A[6]   * A[3]  * A[16] + A[5] * A[10] * A[3]  * A[16]  +
    A[9]  * A[2]   * A[7]  * A[16] - A[1] * A[10] * A[7]  * A[16]  -
    A[5]  * A[2]   * A[11] * A[16] + A[1] * A[6]  * A[11] * A[16]
)
det(A::FixedMatrix) = det(Matrix(A))

trace(A::FixedMatrix{1,1}) = A[1,1]
trace(A::FixedMatrix{2,2}) = A[1,1] + A[2,2]
trace(A::FixedMatrix{3,3}) = A[1,1] + A[2,2] + A[3,3]
trace(A::FixedMatrix{4,4}) = A[1,1] + A[2,2] + A[3,3] + A[4,4]
trace(A::FixedMatrix) = trace(Matrix(A))

\{m,n,T1,T2}(mat::Mat{m,n,T1}, v::Vec{n,T2}) = inv(mat)*v

@inline inv{T}(A::Mat{1, 1, T}) = @inbounds return Mat{1, 1, T}(inv(A[1]))
@inline function inv{T}(A::Mat{2, 2, T})
    determinant = det(A)
    @inbounds return Mat{2, 2, T}(
        (A[2,2] /determinant, -A[2,1]/determinant),
        (-A[1,2]/determinant, A[1,1] /determinant)
    )
end
@inline function inv{T}(A::Mat{3, 3, T})
    determinant = det(A)
    @inbounds return  Mat{3, 3, T}(
        ((A[2,2]*A[3,3]-A[2,3]*A[3,2]) /determinant,
        -(A[2,1]*A[3,3]-A[2,3]*A[3,1])/determinant,
        (A[2,1]*A[3,2]-A[2,2]*A[3,1]) /determinant),

        (-(A[1,2]*A[3,3]-A[1,3]*A[3,2])/determinant,
        (A[1,1]*A[3,3]-A[1,3]*A[3,1]) /determinant,
        -(A[1,1]*A[3,2]-A[1,2]*A[3,1])/determinant),

        ((A[1,2]*A[2,3]-A[1,3]*A[2,2]) /determinant,
        -(A[1,1]*A[2,3]-A[1,3]*A[2,1])/determinant,
        (A[1,1]*A[2,2]-A[1,2]*A[2,1]) /determinant)
    )
end

@inline function inv{T}(A::Mat{4, 4, T})
    determinant = det(A)
    @inbounds return Mat{4, 4, T}(
        ((A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) / determinant,
        (A[2,4]*A[3,3]*A[4,1] - A[2,3]*A[3,4]*A[4,1] - A[2,4]*A[3,1]*A[4,3] + A[2,1]*A[3,4]*A[4,3] + A[2,3]*A[3,1]*A[4,4] - A[2,1]*A[3,3]*A[4,4]) / determinant,
        (A[2,2]*A[3,4]*A[4,1] - A[2,4]*A[3,2]*A[4,1] + A[2,4]*A[3,1]*A[4,2] - A[2,1]*A[3,4]*A[4,2] - A[2,2]*A[3,1]*A[4,4] + A[2,1]*A[3,2]*A[4,4]) / determinant,
        (A[2,3]*A[3,2]*A[4,1] - A[2,2]*A[3,3]*A[4,1] - A[2,3]*A[3,1]*A[4,2] + A[2,1]*A[3,3]*A[4,2] + A[2,2]*A[3,1]*A[4,3] - A[2,1]*A[3,2]*A[4,3]) / determinant),

        ((A[1,4]*A[3,3]*A[4,2] - A[1,3]*A[3,4]*A[4,2] - A[1,4]*A[3,2]*A[4,3] + A[1,2]*A[3,4]*A[4,3] + A[1,3]*A[3,2]*A[4,4] - A[1,2]*A[3,3]*A[4,4]) / determinant,
        (A[1,3]*A[3,4]*A[4,1] - A[1,4]*A[3,3]*A[4,1] + A[1,4]*A[3,1]*A[4,3] - A[1,1]*A[3,4]*A[4,3] - A[1,3]*A[3,1]*A[4,4] + A[1,1]*A[3,3]*A[4,4]) / determinant,
        (A[1,4]*A[3,2]*A[4,1] - A[1,2]*A[3,4]*A[4,1] - A[1,4]*A[3,1]*A[4,2] + A[1,1]*A[3,4]*A[4,2] + A[1,2]*A[3,1]*A[4,4] - A[1,1]*A[3,2]*A[4,4]) / determinant,
        (A[1,2]*A[3,3]*A[4,1] - A[1,3]*A[3,2]*A[4,1] + A[1,3]*A[3,1]*A[4,2] - A[1,1]*A[3,3]*A[4,2] - A[1,2]*A[3,1]*A[4,3] + A[1,1]*A[3,2]*A[4,3]) / determinant),

        ((A[1,3]*A[2,4]*A[4,2] - A[1,4]*A[2,3]*A[4,2] + A[1,4]*A[2,2]*A[4,3] - A[1,2]*A[2,4]*A[4,3] - A[1,3]*A[2,2]*A[4,4] + A[1,2]*A[2,3]*A[4,4]) / determinant,
        (A[1,4]*A[2,3]*A[4,1] - A[1,3]*A[2,4]*A[4,1] - A[1,4]*A[2,1]*A[4,3] + A[1,1]*A[2,4]*A[4,3] + A[1,3]*A[2,1]*A[4,4] - A[1,1]*A[2,3]*A[4,4]) / determinant,
        (A[1,2]*A[2,4]*A[4,1] - A[1,4]*A[2,2]*A[4,1] + A[1,4]*A[2,1]*A[4,2] - A[1,1]*A[2,4]*A[4,2] - A[1,2]*A[2,1]*A[4,4] + A[1,1]*A[2,2]*A[4,4]) / determinant,
        (A[1,3]*A[2,2]*A[4,1] - A[1,2]*A[2,3]*A[4,1] - A[1,3]*A[2,1]*A[4,2] + A[1,1]*A[2,3]*A[4,2] + A[1,2]*A[2,1]*A[4,3] - A[1,1]*A[2,2]*A[4,3]) / determinant),

        ((A[1,4]*A[2,3]*A[3,2] - A[1,3]*A[2,4]*A[3,2] - A[1,4]*A[2,2]*A[3,3] + A[1,2]*A[2,4]*A[3,3] + A[1,3]*A[2,2]*A[3,4] - A[1,2]*A[2,3]*A[3,4]) / determinant,
        (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) / determinant,
        (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) / determinant,
        (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) / determinant)
    )
end

inv(A::FixedMatrix) = typeof(A)(inv(Matrix(A)))

lyap{T}(a::Mat{1, 1, T}, c::Mat{1, 1, T}) = Mat{1,1,T}(lyap(a[1,1],c[1,1]))
function lyap{T}(a::Mat{2, 2, T}, c::Mat{2, 2, T})
    d = det(a)
    t = trace(a)
     -(d*c  + (a - t*I)*c*(a-t*I)')/(2*d*t) # http://www.nber.org/papers/w8956.pdf
end
lyap{m,T}(a::Mat{m,m,T},c::Mat{m,m,T}) = Mat(lyap(Matrix(a),Matrix(c)))

chol{T<:Base.LinAlg.BlasFloat}(m::Mat{1, 1, T}) = Mat{1,1}(chol(m[1,1]))
function chol{T<:Base.LinAlg.BlasFloat}(m::Mat{2,2,T})
    m[1,2]==m[2,1]' || error("Matrix not symmetric")
    l11 = chol(m[1,1])
    @inbounds return Mat{2, 2, T}(
        (l11, zero(T)),
        (inv(l11)*m[1,2], chol(m[2,2] - m[2,1]*inv(m[1,1])*m[1,2]))
    )
end
chol{n,T<:Base.LinAlg.BlasFloat}(m::Mat{n,n,T}) = Mat{n,n,T}(full(Base.LinAlg.chol!(Matrix(m))))
chol!(m::Mat, ::Type{UpperTriangular}) = chol(m) # for 0.4
_chol!(m::Mat, ::Type{UpperTriangular}) = chol(m) # for 0.5
chol!(m::Mat, ::Type{Val{:U}}) = chol!(m, UpperTriangular) # for pre-0.5

# Matrix products
# General shape mismatched versions are errors
(*){T1, T2, M, N, O, K}(a::FixedMatrix{M, N, T1}, b::FixedMatrix{O, K, T2}) = throw(DimensionMismatch("$N != $O in $(typeof(a)) and $(typeof(b))"))
(*){T1, T2, M, N, O}(a::FixedMatrix{M, N, T1}, b::FixedVector{O, T2}) = throw(DimensionMismatch("$N != $O in $(typeof(a)) and $(typeof(b))"))
# vector * (row vector)
@generated function *{T1, T2, N, M}(a::FixedVector{N, T1}, b::FixedMatrix{1,M,T2})
    elements = Expr(:tuple, [Expr(:tuple, [:(a[$i] * b[$j]) for i in 1:N]...) for j in 1:M]...)
    :(construct_similar($b, $elements))
end
# matrix * vector
@generated function *{T1, T2, M, N}(a::FixedMatrix{M, N, T1}, b::FixedVector{N, T2})
    total_terms = M*N
    if total_terms <= 64
        # Full unrolling
        elements = Expr(:tuple, [Expr(:call, :+, [:(a[$i,$k]*b[$k]) for k = 1:N]...) for i in 1:M]...)
    else
        # Expand as a bunch of dot products
        elements = Expr(:tuple, [:(bilindot(Vec(row(a,$i)),b)) for i in 1:M]...)
    end
    :(construct_similar($b, $elements))
end
function *(a::AbstractMatrix, b::FixedVector)
    a*Vector(b)
end
@generated function *{T1, T2, M, N}(a::FixedMatrix{M, N, T1}, b::AbstractVector{T2})
    quote
        length(b) == $N || throw(DimensionMismatch("$b is wrong size - expecting vector of length $N"))
        return a*Vec{N,T2}(b)
    end
end
# matrix * matrix
@generated function *{T1, T2, M, N, P}(a::FixedMatrix{M, N, T1}, b::FixedMatrix{N, P, T2})
    total_terms = M*N*P
    if total_terms <= 64  # 4x4 * 4x4
        # Full unrolling
        elements = Expr(:tuple, [Expr(:tuple, [Expr(:call, :+, [:(a[$i,$k]*b[$k,$j]) for k = 1:N]...) for i in 1:M]...) for j in 1:P]...)
    else
        # Expand as a bunch of mat*vec expressions
        elements = Expr(:tuple, [:(Tuple(a*Vec(column(b,$j)))) for j = 1:P]...)
    end
    :(construct_similar($a, $elements))
end
function *(a::AbstractMatrix, b::FixedMatrix)
    a*Matrix(b)
end
function *(a::FixedMatrix, b::AbstractMatrix)
    Matrix(a)*b
end

# Predicates
function (==)(a::FixedVectorNoTuple, b::FixedVectorNoTuple)
    s_a = size(a)
    s_b = size(b)
    s_a == s_b || return false
    @inbounds for i = 1:length(a)
        a[i] == b[i] || return false
    end
    true
end
(==)(a::FixedArray, b::FixedArray) = Tuple(a) == Tuple(b)

(==){R, T, FSA <: FixedVector}(a::FSA, b::Mat{R, 1, T}) = Tuple(a) == column(b, 1)
(==){R, T, FSA <: FixedVector}(a::Mat{R, 1, T}, b::FSA) = column(a,1) == Tuple(b)
function (==)(a::FixedArray, b::AbstractArray)
    s_a = size(a)
    s_b = size(b)
    s_a == s_b || return false
    @inbounds for i = 1:length(a)
        a[i] == b[i] || return false
    end
    true
end

(==)(a::AbstractArray, b::FixedArray) = b == a

isinteger{FSA<:FixedArray}(A::FSA) = all(isinteger, A)
isinteger{T<:Integer}(A::FixedArray{T}) = true

isreal{FSA<:FixedArray}(A::FSA) = all(isreal, A)
isreal{T<:Real}(A::FixedArray{T}) = true

ishermitian{FSA<:FixedMatrix}(A::FSA) = (A == A')

# To support @test_approx_eq
Base.Test.approx_full(a::FixedArray) = a

# UniformScaling

*(J::Base.LinAlg.UniformScaling, A::FixedArray) = J.λ*A
*(A::FixedArray, J::Base.LinAlg.UniformScaling) = A*J.λ
/(A::FixedArray, J::Base.LinAlg.UniformScaling) = A/J.λ

+{m, n, T}(A::Mat{m,n, T}, J::Base.LinAlg.UniformScaling) = A + J.λ*eye(Mat{m,n,T})
+{m, n, T}(J::Base.LinAlg.UniformScaling, A::Mat{m,n, T}) = A + J
-{m, n, T}(A::Mat{m,n, T}, J::Base.LinAlg.UniformScaling) = A + (-J)
-{m, n, T}(J::Base.LinAlg.UniformScaling, A::Mat{m,n, T}) = J.λ*eye(Mat{m,n,T}) - A

# typed floor, ceil, round, trunc
Base.floor{T}(::Type{T}, A::FixedArray) = map(x->floor(T, x), A)
Base.ceil{T}( ::Type{T}, A::FixedArray) = map(x->ceil( T, x), A)
Base.trunc{T}(::Type{T}, A::FixedArray) = map(x->trunc(T, x), A)
Base.round{T}(::Type{T}, A::FixedArray) = map(x->round(T, x), A)

import Base.Cartesian.@nif
shift(t::Tuple) = Base.tail(t)
unshift(t::Tuple, val) = tuple(val, t...)
push(t::Tuple, val) = tuple(t..., val)
@generated function pop{N,T}(t::NTuple{N,T})
    Expr(:tuple, [:(t[$i]) for i in 1:N-1]...)
end

@generated function deleteat{N,T,i}(x::NTuple{N,T}, ::Type{Val{i}})
    (1 <= i <= N) || throw(BoundsError((x,Val{i})))
    Expr(:tuple, [:(x[$j]) for j in deleteat!([1:N...], i)]...)
end
@generated function deleteat{N,T}(x::NTuple{N,T}, i::Int)
    quote
        (1 <= i <= N) || throw(BoundsError(x,i))
        @nif $N d->(i==d) d-> deleteat(x, Val{d})
    end
end

@generated function insert{N,T,i}(t::NTuple{N,T}, ::Type{Val{i}}, item)
    (1 <= i <= N+1) || throw(BoundsError())
    args = Any[:(t[$k]) for k in 1:N]
    insert!(args, i, :item)
    Expr(:tuple, args...)
end
@generated function insert{N,T}(x::NTuple{N,T}, i::Int, item)
    quote
        (1 <= i <= N+1) || throw(BoundsError())
        @nif $(N+1) d->(i==d) d-> insert(x, Val{d}, item)
    end
end
insert(v::FixedVector, i, val) = construct_similar(typeof(v), insert(Tuple(v), i, val))

for f in (:shift, :pop)
    eval(quote
        ($f)(v::FixedVector) = construct_similar(typeof(v), ($f)(Tuple(v)))
    end)
end
for f in (:unshift, :deleteat, :push)
    eval(quote
        ($f)(v::FixedVector, val) = construct_similar(typeof(v), ($f)(Tuple(v), val))
    end)
end

function reverse{T, N}(x::FixedVector{N, T})
    map(typeof(x)) do i
        x[N-i+1]
    end
end


text:
import Distributed: @distributed

"""Returns an array of subsequent values from the Bernoulli map with parameter β, starting
from x over N time steps."""
function bernoulli_map(x::T,β::T,N::Int) where T
    oone = one(T)                   # 1 in format T
    xout = Array{T,1}(undef,N+1)    # preallocate
    xout[1] = x                     # store initial condition
    @inbounds for i in 2:N+1
        x = β*x                     # actual Bernoulli map
        x = x >= oone ? x-oone : x  # x mod 1
        xout[i] = x                 # store iteration
    end
    return xout
end

"""Find the Bernoulli orbit length that is reached starting from x with parameter β.
Recursive function that starts testing for orbits over period length N, but enlarges
that testing period if no orbit was found."""
function orbit_length(  x::T,            # initial condition
                        β::T,            # Bernoulli parameter
                        N::Int) where T  # max period length
    oone = one(T)                        # 1 in format T
    x0 = x                          # new initial condition
    n = 0                           # orbit length (0 = not found yet)
    j = 0                           # iteration counter

    while n == 0 && j < N           # stop when orbit is fousnd (n >0) or when max period is reached
        j += 1
        x = β*x                     # actual Bernoulli map
        x = x >= oone ? x-oone : x  # x mod 1
        n = x0 == x ? j : 0         # check for periodicity
    end

    # recursive call on orbit_length
    if n == 0         # in case no orbit was found in N steps, repeat with 5N steps
        n,x = orbit_length(x,β,5N)
    else
        return n,x    # return orbit length and last x that's on the orbit
    end
end

"""Returns the minimum x on the orbit of length N, starting from x (which is assumed to be on the orbit."""
function orbit_minimum( x::T,               # initial condition on the orbit
                        β::T,               # Bernoulli parameter
                        N::Int) where T     # the period length of the orbit
    oone = one(T)                   # 1 in format T
    xmin = x                        # store the orbit's minimum in xmin
    for i in 2:N+1
        x = β*x
        x = x >= oone ? x-oone : x  # x mod 1
        xmin = x < xmin ? x : xmin  # always store the smallest x on the orbit
    end
    return xmin
end

"""Returns both the length and minimum of the orbit reached from initial condition x."""
function orbit_length_minimum(  x0::T,              # initial condition on the orbit
                                β::Real;            # Bernoulli parameter
                                basin::Symbol=:x0_proportional,
                                N::Int=10) where T  # initial period length to test for
    o,x = orbit_length(x0,β,N)
    x = orbit_minimum(x,β,o)

    # for starting from all floats in [0,1) use x0_proportional basins
    # for starting from random uniform in [0,1) use uniform basins
    b = basin == :x0_proportional ? Float64(eps(x0)) : 1.0
    return Orbit(β,o,x,b)
end

"""For a given number format T, parameter β and n initial conditions, find orbits of the
Bernoulli map and test for their uniqueness."""
function find_orbits(   ::Type{T},                      # Number format
                        β::Real,                        # Bernoulli parameter
                        n::Int=10000;                   # n initial conditions
                        verbose::Bool=false) where T    # feedback on orbits found and time?

    # pre-allocate empty array of orbits
    orbits = Orbit[]
    uint_type = eval(Symbol("UInt"*repr(sizeof(T)*8)))
    Tβ = T(β)

    tic = time()
    orbits = @distributed (reduce_orbits) for i in 0:n-1            # for n ICs calculate orbit lengths & x
        orbit_length_minimum(reinterpret(T,uint_type(i)),β,basin=:x0_proportional)
    end

    orbits = reduce_orbits(orbits)      # convert Orbit to Vector{Orbit} of length 1 otherwise pass
    sort!(orbits)                       # from shortest to longest orbit
    normalise_basins!(orbits)

    if verbose
        toc = time()
        time_elapsed = readable_secs(toc-tic)
        println("Found $(length(orbits)) orbits in $time_elapsed.")
    end
    return orbits
end

"""For a given number format T, parameter β and n initial conditions, find orbits of the
Bernoulli map and test for their uniqueness."""
function find_orbits_rand(  ::Type{T},                      # Number format
                            β::Real,                        # Bernoulli parameter
                            n::Int=10000;                   # n initial conditions
                            verbose::Bool=false) where T    #

    # pre-allocate empty array of orbits
    orbits = Orbit[]
    Tβ = T(β)

    tic = time()
    orbits = @distributed (reduce_orbits) for i in 1:n            # for n ICs calculate orbit lengths & x
        orbit_length_minimum(T(rand()),β,basin=:uniform)
    end

    orbits = reduce_orbits(orbits)      # convert Orbit to Vector{Orbit} of length 1 otherwise pass
    sort!(orbits)                       # from shortest to longest orbit
    normalise_basins!(orbits)

    if verbose
        toc = time()
        time_elapsed = readable_secs(toc-tic)
        println("Found $(length(orbits)) orbits in $time_elapsed.")
    end
    return orbits
end

import StatsBase: fit, Histogram, mean, pweights

"""Convert a vector of orbits to the invariant measure,
a histogram with `nbins` bins in [0,1)."""
function invariant_measure(orbits::Vector{Orbit{T}},nbins::Int=40) where T

    norbits = length(orbits)
    histograms = fill(0.0,norbits,nbins)
    bin_width = 1/nbins
    bin_edges = collect(0:bin_width:1)

    for (i,orbit) in enumerate(orbits)
        X = bernoulli_map(orbit.min,orbit.β,orbit.length-1)
        X = Float64.(X)         # convert to Float64 to avoid type promoting in fit
        H = fit(Histogram,X,bin_edges).weights
        histograms[i,:] = H / sum(H) / bin_width
    end

    # weighted means
    basins = [orbit.b for orbit in orbits]
    return mean(histograms,pweights(basins),dims=1)[1,:]
end

ext: jl

max_stars_repo_name: milankl/BernoulliMap

max_issues_repo_name: milankl/BernoulliMap

max_forks_repo_name: milankl/BernoulliMap


text:
using SparseMatrixLILs
using Base.Test

function test1()
    x0 = sprand(10,15,0.3)
    x = SparseMatrixLIL(x0)
    k = sprand(4,6,0.6)
    x[1:4,1:6] = k
    x0[1:4,1:6] = k
    @test all(x0.==x)
    @test findnz(x0) == findnz(x)
    x[6:9,6:11] = x[6:9,6:11] + k
    x0[6:9,6:11] = x0[6:9,6:11] + k
    @test all(x0.==x)
    @test findnz(x0) == findnz(x)
    x[6:9,6:11] = x[1:4,1:6] + k
    x0[6:9,6:11] = x0[1:4,1:6] + k
    @test all(x0.==x)
    @test findnz(x0) == findnz(x)
    @test vecnorm((x.+sin.(x)) .- (x0.+sin.(x0))) == 0
end

test1()


text:
using DelayEmbeddings, Statistics

"""
    construct_candidate_variables(
        source::Vector{AbstractVector},
        target::Vector{AbstractVector},
        [cond::Vector{AbstractVector}];
        k::Int = 1, include_instantaneous = true,
        τexclude::Union{Int, Nothing} = nothing,
        maxlag::Union{Int, Float64} = 0.05
    ) → ([τs_source, τs_target, τs_cond, ks_targetfuture], [js_source, js_target, js_cond, js_targetfuture])

Construct candidate variables from input time series. `source` is a vector of equal-length time series
assumed to represent the putative source process. `target` and `cond` are the same, but contains time series
of the target process and of the conditional processes, respectively. `k` is the desired prediction lag.

If `include_instantaneous == true`, then the analysis will also consider instantaneous interactions between
the variables.

If `maxlag` is an integer, `maxlag` is taken as the maximum allowed embedding lag. If `maxlag` is a float,
then the maximum embedding lag is taken as `maximum([length.(source); length.(target); length.(cond)])*maxlag`.

If `τexclude` is an integer, all variables whose embedding lag has absolute value equal to `exclude` will be
excluded.
"""
function construct_candidate_variables(source, target, cond;
        k::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "ac_min",
        maxlag::Union{Int, Float64} = 0.05)

    # Ensure all time series are of the same length.
    Ls = [length.(source); length.(target); length.(cond)]
    @assert all(Ls .== maximum(Ls))

    if maxlag isa Int
        τs = 1:maxlag
    else
        τs = 1:ceil(Int, maximum(Ls)*maxlag)
    end

    # Find the maximum allowed embedding lag for each of the candidates.
    τsmax_source = [estimate_delay(s, method_delay, τs) for s in source]
    τsmax_target = [estimate_delay(t, method_delay, τs) for t in target]
    τsmax_cond = [estimate_delay(c, method_delay, τs) for c in cond]

    # Generate candidate set
    startlag = include_instantaneous ? 0 : -1

    τs_source = [[startlag:-1:-τ...,] for τ in τsmax_source]
    τs_target = [[startlag:-1:-τ...,] for τ in τsmax_target]
    τs_cond = [[startlag:-1:-τ...,] for τ in τsmax_cond]

    ks_targetfuture = [k for i in 1:length(target)]
    js_targetfuture = [i for i in length(τs_source)+1:length(τs_source)+length(τs_target)]
    τs = [τs_source..., τs_target..., τs_cond...]
    js = [[i for x in 1:length(τs[i])] for i = 1:length(τs)]

    # Variable filtering, if desired
    if τexclude isa Int
        τs = [filtered_τs(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
        js = [filtered_js(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
    end
    return [τs..., ks_targetfuture], [js..., js_targetfuture]
end

# Usaully, we use all lags from startlag:-\tau_max to construct variables. In some situations,
# we may want to exclude som of those variables.
function filtered_τs(τs::AbstractVector{Int}, js::AbstractVector{Int}, τexclude::Int)
    [τ for τ in τs if abs(τ) != abs.(τexclude)]
end

function filtered_js(τs::AbstractVector{Int}, js::AbstractVector{Int}, τexclude::Int)
    [j for (τ, j) in zip(τs, js) if abs(τ) != abs.(τexclude)]
end

# source & target variant
function construct_candidate_variables(source, target;
        k::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "ac_min",
        maxlag::Union{Int, Float64} = 0.05)

    # Ensure all time series are of the same length.
    Ls = [length.(source); length.(target)]
    @assert all(Ls .== maximum(Ls))

    if maxlag isa Int
        τs = 1:maxlag
    else
        τs = 1:ceil(Int, maximum(Ls)*maxlag)
    end

    # Find the maximum allowed embedding lag for each of the candidates.
    τsmax_source = [estimate_delay(s, method_delay, τs) for s in source]
    τsmax_target = [estimate_delay(t, method_delay, τs) for t in target]

    # Generate candidate set
    startlag = include_instantaneous ? 0 : -1
    τs_source = [[startlag:-1:-τ...,] for τ in τsmax_source]
    τs_target = [[startlag:-1:-τ...,] for τ in τsmax_target]

    ks_targetfuture = [k for i in 1:length(target)]
    js_targetfuture = [i for i in length(τs_source)+1:length(τs_source)+length(τs_target)]
    τs = [τs_source..., τs_target...,]
    js = [[i for x in 1:length(τs[i])] for i = 1:length(τs)]

    # Variable filtering, if desired
    if τexclude isa Int
        τs = [filtered_τs(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
        js = [filtered_js(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
    end

    return [τs..., ks_targetfuture], [js..., js_targetfuture]
end

# source, target & cond variant
function embed_candidate_variables(source, target, cond;
        η::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "mi_min",
        maxlag::Union{Int, Float64} = 0.05)

    τs, js = construct_candidate_variables(source, target, cond, k = η, τexclude = τexclude)

    # TODO: This is more efficient if not using datasets. Re-do manually.
    data = Dataset([source..., target..., cond...,]...,)
    ℰ = genembed(data, ((τs...)...,), ((js...)...,))

    # Get all variables except the target future (which are the last columns of ℰ)
    n_timeseries = size(ℰ, 2)
    n_timeseries_target = length(target)
    Ω = [ℰ[:, i] for i = 1:n_timeseries - n_timeseries_target]
    Y⁺ = ℰ[:, n_timeseries - n_timeseries_target+1:end]

    # We need to keep track of which variables are from the source, because
    # when computing the final TE, we need a marginal which is 𝒮 \ 𝒮_source.
    # Hence, we need to know which indices in `js` correspond to the source.
    idxs_source = 1:length(source)
    idxs_target = length(source)+1:length(source)+length(target)
    idxs_cond = length(source)+length(target)+1:length(source)+length(target)+length(cond)

    return Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond
end

# source & target variant
function embed_candidate_variables(source, target;
        η::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "mi_min",
        maxlag::Union{Int, Float64} = 0.05)

    τs, js = construct_candidate_variables(source, target, k = η, τexclude = τexclude)

    # TODO: This is more efficient if not using datasets. Re-do manually.
    data = Dataset([source..., target...,]...,)
    ℰ = genembed(data, ((τs...)...,), ((js...)...,))

    # Get all variables except the target future (which are the last columns of ℰ)
    n_timeseries = size(ℰ, 2)
    n_timeseries_target = length(target)
    Ω = [ℰ[:, i] for i = 1:n_timeseries - n_timeseries_target]

    Y⁺ = ℰ[:, n_timeseries - n_timeseries_target+1:end]
    idxs_source = 1:length(source)
    idxs_target = length(source)+1:length(source)+length(target)
    idxs_cond = Int[]

    return Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond
end

function optim_te(Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond, est;
        uq = 0.95, nsurr = 100, q = 1, base = 2)

    τs_comb = [(τs...)...,]
    js_comb = [(js...)...,]

    npts = length(Y⁺)
    n_candidate_variables = length(Ω)

    𝒮 = Vector{Vector{Float64}}(undef, 0)
    𝒮_τs = Vector{Int}(undef, 0)
    𝒮_js = Vector{Int}(undef, 0)

    k = 1
    while k <= n_candidate_variables
        n_remaining_candidates = length(Ω)
        CMIs_between_Y⁺_and_candidates = zeros(n_remaining_candidates)

        # At first iteration, only loop through source variable. If no source variable is found that
        # yields significant TE, terminate.
        for i = 1:n_remaining_candidates
            if k == 1 || length(𝒮) == 0
                Cᵢ = Ω[i]
                CMI_Y⁺_Cᵢ =
                    genentropy(Dataset(Y⁺, Dataset(Cᵢ)), est, q = q, base = base) -
                    genentropy(Dataset(Cᵢ), est, q = q, base = base)
            else
                Cᵢ = [Ω[i], 𝒮...]
                CMI_Y⁺_Cᵢ =
                    genentropy(Dataset(Y⁺, Dataset(Cᵢ...,)), est, q = q, base = base) -
                    genentropy(Dataset(Cᵢ...,), est, q = q, base = base)
            end
            CMIs_between_Y⁺_and_candidates[i] = CMI_Y⁺_Cᵢ
        end

        idx = findfirst(x -> x == minimum(CMIs_between_Y⁺_and_candidates), CMIs_between_Y⁺_and_candidates)
        Wₖ = Ω[idx]

        # Test significance of this candidate by using a random permutation test
        CMI_permutations = zeros(nsurr)

        # A circular shift surrogate generator, to exclude effects of autocorrelation
        s = surrogenerator(Wₖ, CircShift(collect(1:npts - 1)))
        #s = surrogenerator(Wₖ, RandomShuffle())

        if k == 1
            cmiₖ = CMIs_between_Y⁺_and_candidates[idx]

            for i = 1:nsurr
                surr_wₖ = s() # Surrogate version of Wₖ
                CMI_permutations[i] = mutualinfo(Y⁺, surr_wₖ, est)
            end
        else
            # Precompute terms that do not change during surrogate loop
            H_Y⁺_𝒮 = genentropy(Dataset(Y⁺, Dataset(𝒮...,)), est, q = q, base = base)

            # ORIGIANL TE
            H_𝒮 = genentropy(Dataset(𝒮...), est, q = q, base = base)
            cmiₖ = H_Y⁺_𝒮 +
                    genentropy(Dataset([Wₖ, 𝒮...,]...,), est, q = q, base = base) -
                    genentropy(Dataset(Y⁺, Dataset([Wₖ, 𝒮...,]...,)), est, q = q, base = base) -
                    H_𝒮

            for i = 1:nsurr
                surr_wₖ = s() # Surrogate version of Wₖ
                CMI_permutations[i] = H_Y⁺_𝒮 +
                    genentropy(Dataset([surr_wₖ, 𝒮...]...,), est, q = q, base = base) -
                    genentropy(Dataset(Y⁺, Dataset([surr_wₖ, 𝒮...]...,)), est, q = q, base = base) -
                    H_𝒮
            end

        end
       # If the candidate passes the significance test
        if cmiₖ > quantile(CMI_permutations, uq)
            # Add the candidate to list of selected candidates
            push!(𝒮, Wₖ)
            push!(𝒮_τs, τs_comb[idx])
            push!(𝒮_js, js_comb[idx])

            # Delete the candidate from the list of remaining candidates
            deleteat!(Ω, idx)
            deleteat!(τs_comb, idx)
            deleteat!(js_comb, idx)

            k = k + 1
        else
            k = n_candidate_variables + 1
        end
    end

    # No variables were selected
    if length(𝒮) == 0
        return 0.0, Int[], Int[], idxs_source, idxs_target, idxs_cond
    end

    # No variables were selected from the source process
    n_source_vars_picked = count(x -> x ∈ idxs_source, 𝒮_js)
    if n_source_vars_picked == 0
        return 0.0, Int[], Int[], idxs_source, idxs_target, idxs_cond
    end

    # No variables were selected from the target or conditional processes.
    𝒮_nonX = [ts for (ts, j) in zip(𝒮, 𝒮_js) if j ∉ idxs_source]
    if length(𝒮_nonX) == 0
        return 0.0, Int[], Int[], idxs_source, idxs_target, idxs_cond
    end

    CE2 = genentropy(Dataset(Y⁺, Dataset(𝒮...,)), est, base = base, q = q) -
        genentropy(Dataset(𝒮...,), est, base = base, q = q)

    CE1 = genentropy(Dataset(Y⁺, Dataset(𝒮_nonX...,)), est, base = base, q = q) -
        genentropy(Dataset(𝒮_nonX...,), est, base = base, q = q)

    CMI = CE1 - CE2
    return CMI, 𝒮_js, 𝒮_τs, idxs_source, idxs_target, idxs_cond

end

process_input(ts::Vector{T}) where T <: Real = [ts]
process_input(ts::AbstractVector{Vector{T}}) where T <: Real = ts


text:
# -*- coding: utf-8 -*-
import Revise
import GLRenderer
import PoseComposition
import Rotations
import FileIO

R = Rotations
P = PoseComposition
GL = GLRenderer

obj_path = joinpath(@__DIR__, "035_power_drill/textured_simple.obj")
texture_path = joinpath(@__DIR__, "035_power_drill/texture_map.png")

camera_intrinsics = GL.CameraIntrinsics(
    640

    , 480,
    1000.0, 1000.0,
    320.0, 240.0,
    0.01, 5.0
)
camera_intrinsics = GL.scale_down_camera(camera_intrinsics, 4)

renderer = GL.setup_renderer(camera_intrinsics, GL.DepthMode())

mesh_data = GL.get_mesh_data_from_obj_file(obj_path)
GL.load_object!(renderer, mesh_data)

depth_image_base = GL.gl_render(renderer, [1], [P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))], P.IDENTITY_POSE)
GL.view_depth_image(depth_image_base)

for _ in 1:1000
    rand_pose = P.Pose(rand(3) * 10.0, R.RotXYZ(rand(3)*2*pi...))
    rand_pose
    depth_image = GL.gl_render(renderer, [1], [rand_pose * P.Pose([0.0, 0.0, 1.0], R.RotXYZ(0.1, 0.4, 0.9))], rand_pose)
    total_mismatches = sum(abs.(depth_image .- depth_image_base) .> 0.01)
    @assert total_mismatches < 3
end


text:
#exports
export PeriDynTest, Test, run, NameParam

Solvers = Dict(:qs => :quasi_static!, :vv => :velocity_verlet!, :min => :minimize!)

abstract type PeriDynTest end

struct NameParam
    name
    args
    kwargs
end

function NameParam(name, args::Tuple{Any})
    NameParam(name, args, Dict())
end

function Base.length(x::NameParam) 1 end
function Base.first(x::NameParam) x end
function Base.length(x::Function) 1 end
function Base.first(x::Function) x end

struct Test <: PeriDynTest
    solver
    geom
    gen_material
    spc_material
    bc
    RM
    f
    blocknames
    function Test(args...; names=nothing)
        if names === nothing
            names = ["Default name" for i in 1:length(args[3])]
        end
        new([length(i)==1 ? [first(i)] : i for i in args]..., names)
    end
end

function run(test::T) where T <: PeriDynTest
    error("No run defined for test of type **$(typeof(test))**")
end

function getproperty_(a, b)
    if typeof(b)==Symbol
        return getproperty(a, b)
    elseif typeof(b) <: Function
        return b
    else
        error("It should be either a Symbol or a Function.")
    end
end

function stage!(test::Test)
    solver = getproperty(PeriDyn, Solvers[first(test.solver).name])
    skwargs = first(test.solver).kwargs
    sargs = first(test.solver).args

    geoms = [PDMesh.create(getproperty_(PDBenchmark, x.name)(x.args...); x.kwargs...) for x in test.geom]

    gen_mat = [getproperty_(PeriDyn, x.name)(y..., x.args...; x.kwargs...)  for (x, y) in zip(test.gen_material, geoms)]

    spc_mat = [getproperty_(PeriDyn, x.name)(x.args...; x.kwargs...) for x in test.spc_material]

    block = [getproperty_(PeriDyn, :PeridynamicsMaterial)(x, y; name=name) for (x,y, name) in zip(gen_mat, spc_mat, test.blocknames)]

    RM = []
    for x in test.RM
        bks = pop!(x.kwargs, :blocks)
        y = [block[i] for i in bks]
        push!(RM, getproperty_(PeriDyn, x.name)(x.args..., y...; x.kwargs...))
    end

    env =  PeriDyn.Env(1, block, RM, Any[], 1.0)

    BC = [getproperty_(PeriDyn, x.name)(x.args[1](env)...; x.kwargs...) for x in test.bc]

    for bc in BC
        push!(env.boundary_conditions, bc)
    end

    for ff in test.f
        ff(env)
    end

    println("\n\n")
    println("=========================")
    println("     Pre Test Report     ")
    println("=========================")
    println("\n=========Solver==========")
    println(solver)
    println("\tArgs: ", sargs)
    println("\tKwargs: ", skwargs)

    println("\n====General Materials====")
    for i in 1:length(gen_mat)
        println("$i.)")
        println("Block name: $(test.blocknames[i])")
        show(gen_mat[i])
    end

    println("\n===Specific Materials ===")
    for i in 1:length(spc_mat)
        println("$i.)")
        println("Block name: $(test.blocknames[i])")
        show(spc_mat[i])
    end

    println("\n===Boundary Conditions===")
    for i in 1:length(BC)
        println("$i.)")
        show(BC[i])
    end

    println("\n====Repulsive Models ====")
    for i in 1:length(RM)
        println("$i.)")
        show(RM[i])
    end

    return env, (env_) -> solver([env_], sargs...; skwargs...)

end

function run!(test::Test; pseudorun=false)

    env, env_solve! = stage!(test)
    if ~pseudorun
        env_solve!(env)
    end

    return env
end

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark


text:
module ReinforcementLearning
export RL
const RL = ReinforcementLearning

using ReinforcementLearningEnvironments

include("extensions/extensions.jl")

using Reexport
include("Utils/Utils.jl")
include("components/components.jl")
include("glue/glue.jl")
end


text:
@testset "Base" begin

@testset "view" begin
    xs = copyto!(ROCVector{Int}(undef, 4), 1, collect(1:4), 1, 4)
    a = view(xs, 1:2)
    b = view(xs, 3:4)
    @test a isa ROCVector{Int}
    @test b isa ROCVector{Int}
    @test collect(a)::Vector{Int} == 1:2
    @test collect(b)::Vector{Int} == 3:4
    @allowscalar begin
        @test a[[1, 2]] == 1:2
        @test b[[1, 2]] == 3:4
    end

    b_cpu = Vector{Int}(undef, 2)
    @test copyto!(b_cpu, 1, b, 1, 2) === b_cpu
    @test b_cpu == 3:4

    c = ROCVector{Int}(undef, 2)
    @test copyto!(c, 1, b, 1, 2) === c
    @test collect(c) == 3:4
end

end


text:
@testset "Push, pop, shift, unshift, etc" begin
    v = @SVector [1, 2, 3]

    @test push(v, 4) === @SVector [1, 2, 3, 4]
    @test pop(v) === @SVector [1, 2]

    @test unshift(v, 0) === @SVector [0, 1, 2, 3]
    @test shift(v) === @SVector [2, 3]

    @test (@inferred insert(v, 2, -2)) === @SVector [1, -2, 2, 3]
    @test (@inferred deleteat(v, 2)) === @SVector [1, 3]
end


text:
using Distributions, Plots; pyplot()

triangDist = TriangularDist(0,2,1)
xGrid = 0:0.1:2
N = 10^6
inverseSampledData = quantile.(triangDist,rand(N))

histogram( inverseSampledData, bins=30, normed=true,
	ylims=(0,1.1), label="Inverse transform data")
plot!( xGrid, pdf.(triangDist,xGrid), c=:red, lw=4,
	xlabel="x", label="PDF", ylabel = "Density", legend=:topright)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Test, PlotReferenceImages

@testset "PlotReferenceImages" begin
    @test !isdir(reference_path(:notabackend, v"1"))
    @test isfile(reference_file(:gr, 1, v"1"))
end


text:
using Test, NumberTheoreticTransforms, Random

Random.seed!(42)

include("ntt.jl")
include("fnt.jl")


text:
Tables.istable(::Type{<:File}) = true
Tables.columnaccess(::Type{<:File}) = true
Tables.schema(f::File)  = Tables.Schema(getnames(f), _eltype.(gettypes(f)))
Tables.columns(f::File) = f
Base.propertynames(f::File) = getnames(f)

function Base.getproperty(f::File, col::Symbol)
    lookup = getfield(f, :lookup)
    return get(lookup, col) do
        getfield(f, col)
    end
end

function Base.copy(c::Column{T}) where {T}
    len = length(c)
    A = Vector{T}(undef, len)
    @simd for i = 1:len
        @inbounds A[i] = c[i]
    end
    return A
end

function Base.copy(c::Column{T, T}) where {T <: Union{String, Union{String, Missing}}}
    len = length(c)
    A = StringVector{T}(undef, len)
    @simd for i = 1:len
        @inbounds A[i] = c[i]
    end
    return A
end

function Base.copy(c::Column{T, S}) where {T <: Union{String, Union{String, Missing}}, S <: Union{PooledString, Union{PooledString, Missing}}}
    len = length(c)
    catg = c.catg
    tape = c.tape
    crefs = c.refs
    if S === PooledString
        refs = Dict{String, UInt32}()
        foreach(x->setindex!(refs, UInt32(x[1]), x[2]), enumerate(crefs))
        values = Vector{UInt32}(undef, len)
        @simd for i = 1:len
            @inbounds values[i] = ref(tape[i])
        end
    else
        if catg
            refs = Dict{String, UInt32}()
            foreach(x->setindex!(refs, UInt32(x[1]), x[2]), enumerate(crefs))
            missingref = UInt32(0)
            values = Vector{UInt32}(undef, len)
        else # Union{PooledString, Missing}
            refs = Dict{Union{String, Missing}, UInt32}()
            foreach(x->setindex!(refs, UInt32(x[1]), x[2]), enumerate(crefs))
            missingref = UInt32(length(refs) + 1)
            refs[missing] = missingref
            values = Vector{UInt32}(undef, len)
        end
        @simd for i = 1:len
            @inbounds v = ref(tape[i])
            @inbounds values[i] = ifelse(v == UInt32(0), missingref, v)
        end
    end
    if catg
        pool = CategoricalPool(refs)
        levels!(pool, sort(levels(pool)))
        A = CategoricalArray{T, 1}(values, pool)
    else
        A = PooledArray(PooledArrays.RefArray(values), refs)
    end
    return A
end

Base.@propagate_inbounds function Base.getindex(c::Column{Missing}, row::Int)
    @boundscheck checkbounds(c, row)
    return missing
end

Base.@propagate_inbounds function Base.getindex(c::Column{T, S}, row::Int) where {T, S}
    @boundscheck checkbounds(c, row)
    @inbounds x = reinterp_func(T)(c.tape[row])
    return x
end

Base.@propagate_inbounds function Base.getindex(c::Column{Union{T, Missing}, S}, row::Int) where {T, S}
    @boundscheck checkbounds(c, row)
    @inbounds x = c.tape[row]
    return ifelse(x === c.sentinel, missing, reinterp_func(T)(x))
end

Base.@propagate_inbounds function Base.getindex(c::Column{T, PooledString}, row::Int) where {T}
    @boundscheck checkbounds(c, row)
    @inbounds x = c.tape[row]
    @inbounds str = c.refs[x]
    return str
end

Base.@propagate_inbounds function Base.getindex(c::Column{T, Union{PooledString, Missing}}, row::Int) where {T}
    @boundscheck checkbounds(c, row)
    @inbounds x = c.tape[row]
    if x == 0
        return missing
    else
        @inbounds y = c.refs[x]
        return y
    end
end

Base.@propagate_inbounds function Base.getindex(c::Column{T, String}, row::Int) where {T}
    @boundscheck checkbounds(c, row)
    @inbounds offlen = c.tape[row]
    s = PointerString(pointer(c.buf, getpos(offlen)), getlen(offlen))
    return escapedvalue(offlen) ? unescape(s, c.e) : String(s)
end

Base.@propagate_inbounds function Base.getindex(c::Column{T, Union{String, Missing}}, row::Int) where {T}
    @boundscheck checkbounds(c, row)
    @inbounds offlen = c.tape[row]
    if missingvalue(offlen)
        return missing
    else
        s = PointerString(pointer(c.buf, getpos(offlen)), getlen(offlen))
        return escapedvalue(offlen) ? unescape(s, c.e) : String(s)
    end
end


text:
"""
    const ByteData = Union{Vector{UInt8}, Base.CodeUnits{UInt8, <: AbstractString}}

Either a Vector of `UInt8` or a `Base.CodeUnit{UInt8}` object. Compressors
should be able to compress both of these types.
"""
const ByteData = Union{Vector{UInt8}, Base.CodeUnits{UInt8, <: AbstractString}}

# ---------------------------
#    AbstractCompressor
# ---------------------------

"""
    AbstractCompressor

A compressor interface type that represent string compressors.

# Mandatory methods
- compressed_length( <: AbstractCompressor, ::InformationDistances.ByteData)

# Optional methods
- compressed_lengths( <: AbstractCompressor, iter)
"""
abstract type AbstractCompressor end

"""
    compressed_length(compressor, s)

The number of resulting bytes when `s` is compressed with `compressor`.

When implementing a subtype `Compressor <: AbstractCompressor` one should implement
`compressed_length(compressor::Compressor, s::InformationDistances.ByteData)

# Examples
```jldoctest
julia> compressed_length(LibDeflateCompressor(), "hello")
10
```
"""
function compressed_length(compressor::AbstractCompressor, s::AbstractString)

    return compressed_length(compressor, collect(codeunits(s)))
end

"""
    compressed_lengths(compressor, iter)

Calculate for each `s` in `iter` the number of resulting bytes when `s` is compressed with `compressor`.

Implementing this method for a specific subtype of `AbstractCompressor` might lead
to some performance improvements as some compressors need to allocate some resources before
compressing, therefore batch processing might lead to performance improvements as
the resources have to be allocated only once.

It is recommended but not necessary to implement this method for a custom subtype
`Compressor <: AbstractCompressor`. The method signature in that case should be
`compressed_lengths(compressor::Compressor, iter)`.

As Julia does not allow one to specify the eltype of an iterator, one should make at least
sure, that the elements of `iter` can be of type `InformationDistances.ByteData` and optionally
could also be of type `AbstractString`.

# Examples
```jldoctest
julia> compressed_lengths(LibDeflateCompressor(), ["hello", "world", "!"])
3-element Array{Int64,1}:
 10
 10
  6
```
"""
function compressed_lengths(compressor, iter)

    return map(s -> compressed_length(compressor, s), iter)
end

# ---------------------------
#    LibDeflateCompressor
# ---------------------------

"""
    LibDeflateCompressor <: AbstractCompressor

A compressor that uses a `LibDeflate.jl` for compressing.

-----

    LibDeflateCompressor(;compresslevel=12)

Create a `LibDeflateCompressor` with compression level `compresslevel`.

# Examples
```jldoctest
julia> LibDeflateCompressor()
LibDeflateCompressor(12)

julia> LibDeflateCompressor(;compresslevel=8)
LibDeflateCompressor(8)
```
"""
struct LibDeflateCompressor <: AbstractCompressor

    compresslevel::Int

    function LibDeflateCompressor(;compresslevel=12)

        return new(compresslevel)
    end
end

function compressed_length(compressor::LibDeflateCompressor, s::ByteData)

    compressor = LibDeflate.Compressor(compressor.compresslevel)

    outvector = Vector{UInt8}(undef, _deflate_maxoutlen(length(s)))

    @GC.preserve outvector s return unsafe_compress!(
        compressor, pointer(outvector), length(outvector), pointer(s), length(s))
end

function compressed_lengths(compressor::LibDeflateCompressor, iter)

    compressor = LibDeflate.Compressor(compressor.compresslevel)

    outvector = UInt8[]

    return map(iter) do s
        if !(s isa ByteData)
            s = codeunits(s)
        end
        resize!(outvector, max(length(outvector), _deflate_maxoutlen(length(s))))

        @GC.preserve outvector s return unsafe_compress!(
            compressor, pointer(outvector), length(outvector), pointer(s), length(s))
    end
end

function _deflate_maxoutlen(inlen)

    # TODO this calculation is taken from libdeflate, it would maybe be better to call the
    # libdeflate_deflate_compress_bound function directly
    max_num_blocks = max(div(inlen, 10000, RoundUp), 1)
    return (5 * max_num_blocks) + inlen + 9
end

# ---------------------------
#   CodecCompressor
# ---------------------------

"""
    CodecCompressor{ <: TranscodingStreams.Codec} <: AbstractCompressor

A compressor that uses a `TranscodingStreams.Codec` for compressing.

-----

    CodecCompressor{C <: TranscodingStreams.Codec}(;kwargs...)

Create a `CodecCompressor` for the codec `C` with a additional keyword arguments passed
to the constructor of that codec.

# Examples
```jldoctest
julia> using CodecXz: XzCompressor

julia> CodecCompressor{XzCompressor}(; level=6)
CodecCompressor{XzCompressor}(Base.Iterators.Pairs(:level => 6))
```
"""
struct CodecCompressor{C <: TranscodingStreams.Codec} <: AbstractCompressor

    kwargs::Base.Iterators.Pairs

    function CodecCompressor{C}(;kwargs...) where {C}

        return new{C}(kwargs)
    end
end

function compressed_lengths(compressor::CodecCompressor{C}, iter) where {C}

    codec = C(;compressor.kwargs...)
    TranscodingStreams.initialize(codec)

    try
        return map(iter) do s
            if !(s isa ByteData)
                s = codeunits(s)
            end
            length(transcode(codec, s))
        end
    catch
        rethrow()
    finally
        TranscodingStreams.finalize(codec)
    end
end

function compressed_length(compressor::CodecCompressor{C}, s::ByteData) where {C}

    codec = C(;compressor.kwargs...)
    TranscodingStreams.initialize(codec)
    try
        return length(transcode(codec, s))
    catch
        rethrow()
    finally
        TranscodingStreams.finalize(codec)
    end
end


text:
include("../../../src/common.jl")

struct Board <: AbstractImage
    board::Matrix{Int8}
end

Base.size(board::Board) = size(board.board)

function Base.show(io::IO, x::Board)
    show(io, x.board)
    println(io, "\n")
end


text:
### operators

# infelicities:
# * LowRankOperator makes all its arguments 2d; does it need to?

import Base: size, *, Ac_mul_B, show, Array, getindex, dot

export IndexingOperator,
       AbstractLowRankOperator, LowRankOperator,
       IndexedLowRankOperator,
       *, Ac_mul_B, size, show, Array, getindex,
       thin_update!

abstract Operator{T} <: AbstractMatrix{T}

function show(io::IO, o::Operator)
  print(io, "$(Array(o))")
end

size{T}(op::Operator{T}, i::Int) = size(op)[i]

##### Indexing operators

type IndexingOperator{T<:Number}<:Operator{T}
	m
	n
	iobs
end
IndexingOperator(m,n,iobs) = IndexingOperator{Float64}(m,n,iobs)

size{T}(op::IndexingOperator{T}) = (length(op.iobs), (op.m, op.n))

*{T}(op::IndexingOperator{T}, X::Array{T}) = X[op.iobs]
Ac_mul_B{T}(op::IndexingOperator{T}, y::AbstractVector) = reshape(sparsevec(op.iobs, y, op.m*op.n), (op.m, op.n))
function Ac_mul_B{T}(op::IndexingOperator{T}, y::AbstractMatrix)
  size(y,2) == 1 || error("multiplication by IndexingOperators is only defined on vectors; got $(typeof(y))")
  reshape(sparsevec(op.iobs, vec(y), op.m*op.n), (op.m, op.n))
end

##### Low rank operators

mult_map = Dict{Symbol, Function}(:N => *, :T => Ac_mul_B)
tmult_map = Dict{Symbol, Function}(:T => *, :N => Ac_mul_B)
t_map = Dict{Symbol, Function}(:T => ctranspose, :N => identity)

abstract AbstractLowRankOperator{T<:Number}<:Operator{T}
type LowRankOperator{T<:Number}<:AbstractLowRankOperator{T}
	factors::Tuple{Vararg{AbstractMatrix{T}}}
	transpose::Tuple{Vararg{Symbol}}
end
function LowRankOperator(a::AbstractArray...; transpose = tuple(fill(:N, length(a))...))
  t = tuple([make_2d(ai) for ai in a]...)
  l = LowRankOperator(t, transpose)
  return l
end
function LowRankOperator{T<:Number}(t::Tuple{Vararg{AbstractMatrix{T}}}; transpose = tuple(fill(:N, length(a))...))
  return LowRankOperator(t,transpose)
end
# make_2d turns column vectors into matrices with 1 column
function make_2d{T}(a::AbstractArray{T,2})
	return a
end
function make_2d{T}(a::AbstractArray{T,1})
	return reshape(a, (length(a), 1))
end
copy(l::AbstractLowRankOperator) = AbstractLowRankOperator(deepcopy(l.factors), deepcopy(l.transpose))

# defining for specific cases just to remove ambiguity w/Base method ugh
# could fix using invoke???
# i can't seem to make a small problematic example, see attempt in test/ambiguous.jl
function *{T}(l::LowRankOperator{T}, x::AbstractVector{T})
	for i in length(l.factors):-1:1
		x = mult_map[l.transpose[i]](l.factors[i], x)
	end
	return x
end
function *{T}(l::LowRankOperator{T}, x::AbstractMatrix{T})
  for i in length(l.factors):-1:1
		x = mult_map[l.transpose[i]](l.factors[i], x)
	end
	return x
end
# function *{T}(x, l::LowRankOperator{T})
# 	xl = copy(l)
# 	unshift!(xl.factors, x)
# 	unshift!(xl.transpose, :N)
# 	return xl
# end

function Ac_mul_B{T}(l::LowRankOperator{T}, x::AbstractArray{T,1})
	for i in 1:length(l.factors)
		x = tmult_map[l.transpose[i]](l.factors[i], x)
	end
	return x
end
function Ac_mul_B{T}(l::LowRankOperator{T}, x::AbstractArray{T,2})
	for i in 1:length(l.factors)
		x = tmult_map[l.transpose[i]](l.factors[i], x)
	end
	return x
end
function size{T}(l::LowRankOperator{T})
	(size(l.factors[1], 1),
	 size(l.factors[end], t_map[l.transpose[end]] == :N ? 2 : 1)) # 2nd dimension if not transposed; 1st if transposed
end

function Array(l::AbstractLowRankOperator)
	a = t_map[l.transpose[end]](l.factors[end])
	for i in length(l.factors)-1:-1:1
		a = mult_map[l.transpose[i]](l.factors[i], a)
	end
	return a
end

function is_svd(X::LowRankOperator)
  if (length(X.factors)==3) && (X.transpose == (:N, :N, :T))
    # X.factors[2] is diagonal
    # X.factors[1] and [3] are orthogonal
    return true
  else
    return false
  end
end
function is_rank_1(X::LowRankOperator)
  if (length(X.factors)==2) && (X.transpose == (:N, :T))
    return true
  else
    return false
  end
end

# compute alpha*X + (1-alpha)*Delta
# update the factorization of X using idea from
# Brand 2006 "Fast Low-Rank Modifications of the Thin SVD"
function thin_update!{T<:Number}(X::LowRankOperator{T}, Delta::LowRankOperator{T}, alpha::T=1.0, beta::T=1.0)
  # assume that the factors of X form an SVD factorization
  @assert is_svd(X)
  # and that Delta is rank 1
  @assert is_rank_1(Delta)

  # scale X and Delta so we just need to compute U * Sigma * V' + a * b'
  U, Sigma, V = X.factors
  Sigma *= alpha
  r = size(Sigma,1)
  a,b = Delta.factors
  a *= beta

  # eqn 6
  m = U'*a
  p = a-U*m
  ra = norm(p)
  P = p/ra

  # eqn 7
  n = V'*b
  q = b-V*n
  rb = norm(q)
  Q = q/rb

  # eqn 8 - can be made more efficient exploiting diag + rank 1
  K = [Sigma zeros(r); zeros(1,r) 0] + [m; ra]*[n' rb]
  Uk,Sk,Vk = svd(K)

  Up = [U P]*Uk
  Vp = [V Q]*Vk
  X.factors = (Up, spdiagm(Sk), Vp)
  return X
end

#### now mix IndexingOperator and LowRankOperator

type IndexedLowRankOperator{T<:Number}<:AbstractLowRankOperator{T}
	factors::Tuple{IndexingOperator{T},AbstractVector{T}}
  transpose::Tuple{Vararg{Symbol}}
end
IndexedLowRankOperator(a...) = IndexedLowRankOperator{Float64}((a...), (:T, :N))
function A_mul_B!{T}(u::AbstractArray{T,1}, l::IndexedLowRankOperator{T}, v::AbstractVector{T})
  @assert(length(u)==l.factors[1].m)
  @assert(length(v)==l.factors[1].n)
	iobs = l.factors[1].iobs
  js = round(Int,floor((iobs-1)/l.factors[1].m)+1)
  is = (iobs-1)%l.factors[1].m+1
  for ii in length(iobs)
    u[is[ii]] += l.factors[2][ii]*v[js[ii]]
  end
	return u
end
*{T}(l::IndexedLowRankOperator{T}, v::AbstractVector{T}) =
    A_mul_B!(zeros(l.factors[1].m), l, v)
function A_mul_B!{T}(v::AbstractMatrix{T}, l::IndexedLowRankOperator{T}, u::AbstractMatrix{T})
    for i=1:size(u,2)
      A_mul_B!(view(v,:,i), l, view(u,:,i))
    end
    v
end
*{T}(l::IndexedLowRankOperator{T}, v::AbstractMatrix{T}) =
    A_mul_B!(zeros(l.factors[1].m, size(v,2)), l, v)
function Ac_mul_B!{T}(v::AbstractArray{T,1}, l::IndexedLowRankOperator{T}, u::AbstractArray{T,1})
  @assert(length(u)==l.factors[1].m)
  @assert(length(v)==l.factors[1].n)
	iobs = l.factors[1].iobs
  js = round(Int,floor((iobs-1)/l.factors[1].m)+1)
  is = (iobs-1)%l.factors[1].m+1
  v = zeros(l.factors[1].n)
  for ii in length(iobs)
    v[js[ii]] += l.factors[2][ii]*u[is[ii]]
  end
	return v
end
Ac_mul_B{T}(l::IndexedLowRankOperator{T}, u::AbstractVector{T}) =
    Ac_mul_B!(zeros(l.factors[1].n), l, u)
function Ac_mul_B!{T}(v::AbstractMatrix{T}, l::IndexedLowRankOperator{T}, u::AbstractMatrix{T})
    for i=1:size(u,2)
      Ac_mul_B!(view(v,:,i), l, view(u,:,i))
    end
    v
end
Ac_mul_B{T}(l::IndexedLowRankOperator{T}, u::AbstractMatrix{T}) =
    Ac_mul_B!(zeros(l.factors[1].n, size(u,2)), l, u)
function size{T}(l::IndexedLowRankOperator{T})
  size(l.factors[1])[2]
end
# not sure why this is needed for svds
function getindex(l::IndexedLowRankOperator,i::Int,j::Int)
  idx = findin(l.factors[1].iobs, i+l.factors[1].m*(j-1))
  if length(idx)==0
    return 0
  else
    return sum(l.factors[2][idx])
  end
end
function Array(l::IndexedLowRankOperator)
  iobs = l.factors[1].iobs
  js = round(Int,floor((iobs-1)/l.factors[1].m)+1)
  is = (iobs-1)%l.factors[1].m+1
  return sparse(is, js, l.factors[2], l.factors[1].m, l.factors[1].n)
end

function *{T}(iop::IndexingOperator{T}, lrop::LowRankOperator{T})
	# we only know how to do this for rank 1 LowRankOperators, for now
	if length(lrop.factors) == 2 && size(lrop.factors[1])[2] == 1 & size(lrop.factors[2])[1] == 1
		u = lrop.factors[1]
		v = lrop.factors[2]
		y = zeros(length(iop.iobs))
		for idx in 1:length(iop.iobs)
			i = (iop.iobs[idx]-1) % iop.m + 1
			j = floor(Int, (iop.iobs[idx]-1) / iop.m + 1)
			y[idx] = u[i]*v[j]
		end
		return y
	else
		warn("materializing LowRankOperator...")
		return iop*Array(lrop)
		# error("We only know how to multiply IndexingOperators by LowRankOperators for rank 1 LowRankOperators, for now")
	end
end

function getindex(op::LowRankOperator, i::Int, j::Int)
	if length(op.factors)==2
		if op.transpose[1]==:N
			xi = op.factors[1][i,:]
		else
			xi = op.factors[1][:,i]
		end
		if op.transpose[2]==:N
    	yj = op.factors[2][:,j]
		else
			yj = op.factors[2][j,:]
		end
		return dot(vec(xi), vec(yj))
	elseif is_svd(op)
    return dot(op.factors[1][i,:], op.factors[2]*op.factors[3][j,:])
  else
		error("indexing not defined for LowRankOperator with $(length(op.factors)) factors")
	end
end

*(a::Number, o::AbstractLowRankOperator) = (scale!(o.factors[end], a); o)
*(o::AbstractLowRankOperator, a::Number) = *(a,o)

function dot(A::SparseMatrixCSC, lrop::LowRankOperator)
  # @assert is_svd(lrop)
  # @assert size(A) == size(lrop)
  # construct U, V so that U*V = lrop
  U = lrop.factors[1]
  n,r = size(lrop.factors[3])
  V = Array(Float64,r,n)
  for i=1:r
    V[i,:] = lrop.factors[2].nzval[i] * lrop.factors[3][:,i]
  end
  out = 0
  m,n = size(A)
  for j=1:n
    for iobs = A.colptr[j]:(A.colptr[j+1]-1)
      i = A.rowval[iobs]
      out += dot(U[i,:],V[:,j])*A.nzval[iobs]
    end
  end
  return out
end

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
struct GeometricBrownianMotion{T1,T2}
  μ::T1
  σ::T2
end

function (p::GeometricBrownianMotion)(W,dt,rng) #dist
  drift = p.μ-(1/2)*p.σ^2
  if typeof(W.dW) <: AbstractArray
    rand_val = wiener_randn(rng,W.dW)
  else
    rand_val = wiener_randn(rng,typeof(W.dW))
  end
  new_val = @.. exp(drift*dt + p.σ*sqrt(dt)*rand_val)
  return W[end]*(new_val-1)
end

#=
u = q*h
t = h
s = 0

t-u = h-qh = (1-q)h
t-s = h
u - s = q*h
https://math.stackexchange.com/questions/412470/conditional-distribution-in-brownian-motion

=#
function gbm_bridge(gbm,W,W0,Wh,q,h,rng)
  if typeof(W.dW) <: AbstractArray
    return gbm.σ*sqrt((1-q)*q*abs(h))*wiener_randn(rng,W.dW)+q*Wh
  else
    return gbm.σ*sqrt((1-q)*q*abs(h))*wiener_randn(rng,typeof(W.dW))+q*Wh
  end
end
function gbm_bridge!(rand_vec,gbm,W,W0,Wh,q,h,rng)
  wiener_randn!(rng,rand_vec)
  @.. rand_vec = gbm.σ*sqrt((1-q)*q*abs(h))*rand_vec+q*Wh
end

function GeometricBrownianMotionProcess(μ,σ,t0,W0,Z0=nothing;kwargs...)
  gbm = GeometricBrownianMotion(μ,σ)
  NoiseProcess(t0,W0,Z0,gbm,(W,W0,Wh,q,h,rng)->gbm_bridge(gbm,W,W0,Wh,q,h,rng);kwargs...)
end

struct GeometricBrownianMotion!{T1,T2}
  μ::T1
  σ::T2
end
function (p::GeometricBrownianMotion!)(rand_vec,W,dt,rng) #dist!
  wiener_randn!(rng,rand_vec)
  @.. rand_vec = W[end]*(exp(p.μ-(1/2)*p.σ*dt + p.σ*sqrt(dt)*rand_vec)-1)
end
function GeometricBrownianMotionProcess!(μ,σ,t0,W0,Z0=nothing;kwargs...)
  gbm = GeometricBrownianMotion!(μ,σ)
  NoiseProcess(t0,W0,Z0,gbm,(rand_vec,W,W0,Wh,q,h,rng)->gbm_bridge!(rand_vec,gbm,W,W0,Wh,q,h,rng);kwargs...)
end


text:
using Fermi.HartreeFock

export RCCSD

abstract type RCCSDAlgorithm end

function get_rccsd_alg()
    implemented = [RCCSDa()]
    N = Options.get("cc_alg")
    try
        return implemented[N]
    catch BoundsError
        throw(FermiException("implementation number $N not available for RCCSD."))
    end
end

"""
    Fermi.CoupledCluster.RCCSD
    TODO

Fermi struct that holds information about RCCSD wavefunctions

# Fields

    CorrelationEnergy   CCSD correlation energy
    T1                  T1 amplitudes
    T2                  T2 amplitudes

_struct tree:_

**RCCSD** <: AbstractCCWavefunction <: AbstractCorrelatedWavefunction <: AbstractWavefunction
"""
struct RCCSD{T} <: AbstractCCWavefunction
    guessenergy::Float64
    energy::Float64
    correlation::T
    T1::AbstractArray{T,2}
    T2::AbstractArray{T,4}
    δE::T
    residue::T
end

function RCCSD()
    aoints = IntegralHelper{Float64}()
    #aoints = IntegralHelper{Float64}(eri_type=JKFIT())
    rhf = RHF(aoints)
    moints = IntegralHelper(orbitals=rhf.orbitals)
    RCCSD(moints, aoints)
end

function RCCSD(moints::IntegralHelper{T1,E1,O}, aoints::IntegralHelper{T2,E2,AtomicOrbitals}) where {T1<:AbstractFloat,
                                                                                T2<:AbstractFloat,O<:AbstractRestrictedOrbitals,
                                                                                E1<:AbstractERI,E2<:AbstractERI}
    mo_from_ao!(moints, aoints, "Fd","OOOO", "OOOV", "OOVV", "OVOV", "OVVV", "VVVV")
    RCCSD(moints)
end

function RCCSD(moints::IntegralHelper{T1,E1,O}, aoints::IntegralHelper{T2,E2,AtomicOrbitals}) where {T1<:AbstractFloat,T2<:AbstractFloat,
                                                                                E1<:AbstractDFERI,E2<:AbstractDFERI,O<:AbstractRestrictedOrbitals}
    mo_from_ao!(moints, aoints, "Fd","BOO", "BOV", "BVV")
    RCCSD(moints)
end

function RCCSD{Float64}(mol = Molecule(), ints = IntegralHelper{Float64}())

    # Compute Restricted Hartree-Fock
    refwfn = RHF(mol, ints)

    # Delete integrals that are not gonna be used anymore
    delete!(ints, "S", "T", "V", "JKERI")
    RCCSD{Float64}(refwfn, ints)
end

function RCCSD{Float32}(mol = Molecule(), ints = IntegralHelper{Float32}())

    # Note that using this method a new IntegralHelper object
    # is created within the RHF call. This is necessary becasue
    # RHF needs a double precision integral helper

    # Compute Restricted Hartree-Fock
    refwfn = RHF(mol)

    RCCSD{Float32}(refwfn, ints)
end

function RCCSD(moints::IntegralHelper{T,E,O}) where {T<:AbstractFloat,E<:AbstractERI,O<:AbstractRestrictedOrbitals}

    # Create zeroed guesses for amplitudes

    o = moints.molecule.Nα - Options.get("drop_occ")
    v = size(moints.orbitals.C,1) - Options.get("drop_vir") - moints.molecule.Nα

    output("Using MP2 guess")
    T1guess = moints["Fia"]
    T2guess = permutedims(moints["OVOV"], (1,3,2,4))

    # Orbital energies line
    if haskey(moints.cache, "D1")
        d = moints["D1"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        d = FermiMDArray([ϵo[i]-ϵv[a] for i=eachindex(ϵo), a=eachindex(ϵv)])
        moints["D1"] = d
    end

    if haskey(moints.cache, "D2")
        D = moints["D2"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        D = FermiMDArray([ϵo[i]+ϵo[j]-ϵv[a]-ϵv[b] for i=eachindex(ϵo), j=eachindex(ϵo), a=eachindex(ϵv), b=eachindex(ϵv)])
        moints["D2"] = D
    end

    T1guess ./= d
    T2guess ./= D

    RCCSD(moints, T1guess, T2guess, get_rccsd_alg())
end

# For each implementation a singleton type must be create
struct RCCSDa <: RCCSDAlgorithm end
include("RCCSDa.jl")


text:
function test_uniform_grid_type(T, N, H)
  grid = UniformGrid(T, N, H)
  T == typeof(grid.H) && T == eltype(grid.zc) && T == eltype(grid.zf)
end

function test_uniform_grid_spacing(T, N, H)
  Δ = T(H/N)
  grid = UniformGrid(T, N, H)
  Δ == grid.Δc && Δ == grid.Δf
end

function test_uniform_grid_limits_zc_left(T, N, H)
  Δ = H/N
  grid = UniformGrid(T, N, H)
  grid.zc[1] ≈ -H+0.5Δ
end

function test_uniform_grid_limits_zc_right(T, N, H)
  Δ = H/N
  grid = UniformGrid(T, N, H)
  grid.zc[N] ≈ -0.5Δ
end

function test_uniform_grid_limits_zf_left(T, N, H)
  Δ = H/N
  grid = UniformGrid(T, N, H)
  grid.zf[1] ≈ -H
end

function test_uniform_grid_limits_zf_right(T, N, H)
  Δ = H/N
  grid = UniformGrid(T, N, H)
  grid.zf[N+1] ≈ 0
end

@testset "Grids" begin
    N, H = 3, 4.2
    for T in (Float64, Float32, Float16)
        @test test_uniform_grid_type(T, N, H)
        @test test_uniform_grid_spacing(T, N, H)
        @test test_uniform_grid_limits_zc_right(T, N, H)
        @test test_uniform_grid_limits_zf_right(T, N, H)
        @test test_uniform_grid_limits_zc_left(T, N, H)
        @test test_uniform_grid_limits_zf_left(T, N, H)
    end
end


text:
using MinFEM
using LinearAlgebra

function semilinear(mesh::Mesh, L::AbstractMatrix, M::AbstractMatrix,
                    s::AbstractVector, BoundaryIndices::Set{Int64}=[], tol=1e-10)

  y = zeros(mesh.nnodes)

  pde = PDESystem(A=L, b=M*s, bc=zeros(mesh.nnodes), DI=BoundaryIndices)

  res = Inf
  while res > tol
    pde.A = L + asmCubicDerivativeMatrix(mesh, y)
    pde.b = -L*y + M*s - asmCubicTerm(mesh, y)
    refresh(pde)
    solve(pde)

    y += pde.state
    res = norm(pde.state)
    println(res)
  end
  return y
end

mesh = import_mesh("../meshes/semilinear.msh")

L = asmLaplacian(mesh)
M = asmMassMatrix(mesh)

# y = 3*sin(x[1]*pi)*sin(x[2]*pi)
f(x) = 3*2*pi^2*sin(x[1]*pi)*sin(x[2]*pi) + (3*sin(x[1]*pi)*sin(x[2]*pi))^3
s = evaluateMeshFunction(mesh, f)

boundary = union(mesh.Boundaries[1001].Nodes,
                 mesh.Boundaries[1002].Nodes,
                 mesh.Boundaries[1003].Nodes,
                 mesh.Boundaries[1004].Nodes);

y = semilinear(mesh, L, M, s, boundary);

vtkfile = open_vtk_file(mesh, "semilinear.vtu")
write_point_data(vtkfile, y, "y")
write_point_data(vtkfile, s, "s")
save_vtk_file(vtkfile)


text:
"""
    StateActionReward(m::Union{MDP,POMDP})

Robustly create a reward function that depends only on the state and action.

If `reward(m, s, a)` is implemented, that will be used, otherwise the mean of `reward(m, s, a, sp)` for MDPs or `reward(m, s, a, sp, o)` for POMDPs will be used.

# Example
```jldoctest
using POMDPs
using POMDPModels
using POMDPModelTools

m = BabyPOMDP()

rm = StateActionReward(m)

rm(true, true)

# output

-15.0
```
"""
abstract type StateActionReward end

function StateActionReward(m)
    if hasmethod(reward, Tuple{typeof(m), statetype(m), actiontype(m)})
        return FunctionSAR(m)
    else
        return LazyCachedSAR(m)
    end
end

struct FunctionSAR{M} <: StateActionReward
    m::M
end

function (sar::FunctionSAR)(s, a)
    if isterminal(sar.m, s)
        return 0.0
    else
        return reward(sar.m, s, a)
    end
end

struct LazyCachedSAR{M} <: StateActionReward
    m::M
    cache::Matrix{Union{Missing,Float64}}
end

function LazyCachedSAR(m)
    ns = length(states(m))
    na = length(actions(m))
    return LazyCachedSAR(m, Matrix{Union{Missing, Float64}}(missing, ns, na))
end

function (sar::LazyCachedSAR)(s, a)::Float64
    si = stateindex(sar.m, s)
    ai = actionindex(sar.m, a)
    entry = sar.cache[si, ai]
    if ismissing(entry)
        r = mean_reward(sar.m, s, a)
        sar.cache[si, ai] = r
    else
        r = entry
    end
    return r
end

function mean_reward(m::MDP, s, a)
    if isterminal(m, s)
        return 0.0
    else
        td = transition(m, s, a)
        rsum = 0.0
        wsum = 0.0
        for (sp, w) in weighted_iterator(td)
            rsum += w*reward(m, s, a, sp)
            wsum += w
        end
        return rsum/wsum
    end
end

function mean_reward(m::POMDP, s, a)
    if isterminal(m, s)
        return 0.0
    else
        td = transition(m, s, a)
        rsum = 0.0
        wsum = 0.0
        for (sp, w) in weighted_iterator(td)
            od = observation(m, s, a, sp)
            for (o, ow) in weighted_iterator(od)
                rsum += ow*w*reward(m, s, a, sp, o)
                wsum += ow*w
            end
        end
        return rsum/wsum
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using HDF5
using LightXML
importall LightXML
importall Base

function new_element(name::AbstractString, attrs::Dict)
    x = new_element(name)
    for (k, v) in attrs
        x[k] = v
    end
    return x
end

function setindex!(x::XMLElement, content::Any, attr_name::AbstractString)
    set_attribute(x, attr_name, content)
end

function haskey(x::XMLElement, key::AbstractString)
    return has_child(x, key) || has_attribute(x, key)
end

function has_child(x::XMLElement, child_name::AbstractString)
    return get_child(x, child_name) != nothing
end

function get_attribute(x::XMLElement, attr_name::AbstractString)
    attr = attribute(x, attr_name)
    numeric = tryparse(Int64, attr)
    isnull(numeric) && (numeric = tryparse(Float64, attr))
    isnull(numeric) && return attr
    return get(numeric)
end

function new_child(xparent::XMLElement, name::AbstractString, attrs::Dict)
    x = new_child(xparent, name)
    for (k, v) in attrs
        x[k] = v
    end
    return x
end

function new_child(xparent::XMLElement, name::AbstractString, attrs::Pair...)
    x = new_child(xparent, name)
    for (k, v) in attrs
        x[k] = v
    end
    return x
end

""" Basic traverse support, so that it's possible to find data from xml using
path syntax e.g. /foo/bar[2]/baz[@Name=Frame 1]/DataItem. If several elements
with same name exists in tree, pick first by default and next ones can be picked
using [] syntax or [@attr=value] syntax, see [1] for details. For last item use
[end].

[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format
"""
function get_child(x::XMLElement, child_name::AbstractString)
    '/' in child_name && return nothing
    m = match(r"(\w+)\[(.+)\]", child_name)
    if m != nothing
        child_name = m[1]
    end
    childs = []
    for child in child_elements(x)
        if name(child) == child_name
            push!(childs, child)
        end
    end
    length(childs) == 0 && return nothing
    m == nothing && return first(childs)
    j = tryparse(Int, m[2])
    isnull(j) || return childs[get(j)]
    m[2] == "end" && return childs[end]
    m2 = match(r"@(.+)=(.+)", m[2])
    m2 == nothing && throw("Unable to parse: $(m[2])")
    attr_name = m2[1]
    attr_value = m2[2]
    for child in childs
        has_attribute(child, attr_name) || continue
        if get_attribute(child, attr_name) == attr_value
            return child
        end
    end
    throw("Unable to parse: $(m[2])")
end

function getindex(x::XMLElement, attr_name::AbstractString)
    attr_name = strip(attr_name, '/')
    child = get_child(x, attr_name)
    child == nothing || return child
    has_attribute(x, attr_name) && return get_attribute(x, attr_name)
    if '/' in attr_name
        items = split(attr_name, '/')
        attr_name = first(items)
        length(items) > 1 || throw(KeyError(attr_name))
        haskey(x, attr_name) || throw(KeyError(attr_name))
        path = join(items[2:end], '/')
        new_item = getindex(x, attr_name)
        return new_item[path]
    else
        throw(KeyError(attr_name))
    end
end

type Xdmf
    name :: AbstractString
    xml :: XMLElement
    hdf :: HDF5File
end

function new_child(xdmf::Xdmf, args...; kwargs...)
    new_child(xdmf.xml, args...; kwargs...)
end

function read(xdmf::Xdmf, path::AbstractString)
    result = getindex(xdmf.xml, path)
    if endswith(path, "DataItem")
        format = get_attribute(result, "Format")
        @assert format == "HDF"
        h5file, path = map(ASCIIString, split(content(result), ':'))
        h5file = dirname(xdmf.name) * "/" * h5file
        isfile(h5file) || throw("Xdmf: h5 file $h5file not found!")
        return read(xdmf.hdf, path)
    else
        return result
    end
end

function Xdmf()
    return Xdmf(tempname())
end

function h5file(xdmf::Xdmf)
    return xdmf.name*".h5"
end

function xmffile(xdmf::Xdmf)
    return xdmf.name*".xmf"
end

function Xdmf(name::AbstractString)
    xdmf = new_element("Xdmf")
    set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude")
    set_attribute(xdmf, "Version", "2.1")
    h5file = "$name.h5"
    flag = isfile(h5file) ? "r+" : "w"
    hdf = h5open(h5file, flag)
    return Xdmf(name, xdmf, hdf)
end

function save!(xdmf::Xdmf)
    doc = XMLDocument()
    set_root(doc, xdmf.xml)
    save_file(doc, xmffile(xdmf))
end

function new_dataitem{T,N}(xdmf::Xdmf, path::AbstractString, data::Array{T,N}; format="HDF")
    dataitem = new_element("DataItem")
    datatype = replace("$T", "64", "")
    dimensions = join(size(data), " ")
    set_attribute(dataitem, "DataType", datatype)
    set_attribute(dataitem, "Dimensions", dimensions)
    set_attribute(dataitem, "Format", format)
    if format == "HDF"
        hdf = basename(h5file(xdmf))
        if !exists(xdmf.hdf, path)
            write(xdmf.hdf, path, data)
        end
        add_text(dataitem, "$hdf:$path")
    elseif format == "XML"
        add_text(dataitem, strip(string(data), ['[', ']']))
    end
    return dataitem
end


text:
using ModelingToolkitStandardLibrary.Blocks
using ModelingToolkit, OrdinaryDiffEq

@parameters t

@testset "Gain" begin
    @named c = Constant(; k=1)
    @named gain = Gain(1;)
    @named int = Integrator(; k=1)
    @named model = ODESystem([connect(c.output, gain.input), connect(gain.output, int.input)], t, systems=[int, gain, c])
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>1.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())

    @test sol[int.output.u][end] ≈ 2
end

@testset "Feedback loop" begin
    @named c = Constant(; k=2)
    @named gain = Gain(1;)
    @named int = Integrator(; k=1)
    @named fb = Feedback(;)
    @named model = ODESystem(
        [
            connect(c.output, fb.input1),
            connect(fb.input2, int.output),
            connect(fb.output, gain.input),
            connect(gain.output, int.input),
        ],
        t,
        systems=[int, gain, c, fb]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 100.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 2
end

@testset "Add" begin
    @named c1 = Constant(; k=1)
    @named c2 = Constant(; k=2)
    @named add = Add(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, add.input1),
            connect(c2.output, add.input2),
            connect(add.output, int.input),
        ],
        t,
        systems=[int, add, c1, c2]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 3
end

@testset "Product" begin
    @named c1 = Constant(; k=1)
    @named c2 = Constant(; k=2)
    @named prod = Product(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, prod.input1),
            connect(c2.output, prod.input2),
            connect(prod.output, int.input),
        ],
        t,
        systems=[int, prod, c1, c2]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 2
end

@testset "Division" begin
    @named c1 = Constant(; k=1)
    @named c2 = Constant(; k=2)
    @named div = Division(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, div.input1),
            connect(c2.output, div.input2),
            connect(div.output, int.input),
        ],
        t,
        systems=[int, div, c1, c2]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 1/2
end

@testset "Abs" begin
    @named c = Constant(; k=-1)
    @named abs = Abs(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c.output, abs.input),
            connect(abs.output, int.input),
        ],
        t,
        systems=[int, abs, c]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 1
end

@testset "Sqrt" begin
    @named c = Constant(; k=4)
    @named sqr = Sqrt(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c.output, sqr.input),
            connect(sqr.output, int.input),
        ],
        t,
        systems=[int, sqr, c]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 2
end

@testset "Sign" begin
    @named c = Constant(; k=3)
    @named sig = Sign(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c.output, sig.input),
            connect(sig.output, int.input),
        ],
        t,
        systems=[int, sig, c]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ 1
end

@testset "MatrixGain" begin
    K = [1 2; 3 4]
    @named gain = MatrixGain(K;)
    # TODO:
end

@testset "Sum" begin
    @named s = Sum(2;)
    # TODO:
end

@testset "Math" begin
    for (block, func) in [(Abs, abs), (Sin, sin), (Cos, cos), (Tan, tan), (Asin, asin), (Acos, acos), (Atan, atan), (Sinh, sinh), (Cosh, cosh), (Tanh, tanh), (Exp, exp)]
        @named source = Sine(frequency=1)
        @named b = block()
        @named int = Integrator()
        @named model = ODESystem([connect(source.output, b.input), connect(b.output, int.input)], t, systems=[int, b, source])
        sys = structural_simplify(model)

        prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

        sol = solve(prob, Rodas4())
        @test sol[b.output.u] ≈ func.(sol[source.output.u])
    end

    # input must be positive
    for (block, func) in [(Sqrt, sqrt), (Log, log), (Log10, log10)]
        @named source = Sine(; frequency=1, offset=2)
        @named b = block()
        @named int = Integrator()
        @named model = ODESystem([connect(source.output, b.input), connect(b.output, int.input)], t, systems=[int, b, source])
        sys = structural_simplify(model)

        prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

        sol = solve(prob, Rodas4())
        @test sol[b.output.u] ≈ func.(sol[source.output.u])
    end
end

@testset "Atan2" begin
    @named c1 = Constant(; k=1)
    @named c2 = Constant(; k=2)
    @named b = Atan2(;)
    @named int = Integrator(; k=1)
    @named model = ODESystem(
        [
            connect(c1.output, b.input1),
            connect(c2.output, b.input2),
            connect(b.output, int.input),
        ],
        t,
        systems=[int, b, c1, c2]
    )
    sys = structural_simplify(model)

    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 1.0))

    sol = solve(prob, Rodas4())
    @test sol[int.output.u][end] ≈ atan(1, 2)
end


text:
# Note: some definitions in this file come from PETSc.jl

const PetscErrorCode = Cint
const comm_type = MPI.Comm

const Mat = Ptr{Nothing}

PetscViewer = Ptr{Nothing}

PetscInt = Int32
PetscScalar = Float64
PetscBool = UInt32

const PETSC_FALSE = (UInt32)(0)
const PETSC_TRUE = (UInt32)(1)
const PETSC_NULL = C_NULL
const PETSC_IGNORE = C_NULL
const PETSC_DECIDE = (Int32)(-1)
const PETSC_DETERMINE = PETSC_DECIDE
const PETSC_DEFAULT = (Int32)(-2)

const MatType = Ptr{UInt8}
const MATAIJ = "aij"

const MatAssemblyType = UInt32
const MAT_FINAL_ASSEMBLY = (UInt32)(0)

const InsertMode = UInt32
const NOT_SET_VALUES = (UInt32)(0)
const INSERT_VALUES = (UInt32)(1)
const ADD_VALUES = (UInt32)(2)
const MAX_VALUES = (UInt32)(3)
const INSERT_ALL_VALUES = (UInt32)(4)
const ADD_ALL_VALUES = (UInt32)(5)
const INSERT_BC_VALUES = (UInt32)(6)
const ADD_BC_VALUES = (UInt32)(7)

const MatOption = Cint
const MAT_OPTION_MIN = (Int32)(-5)
const MAT_NEW_NONZERO_LOCATION_ERR = (Int32)(-4)
const MAT_UNUSED_NONZERO_LOCATION_ERR = (Int32)(-3)
const MAT_NEW_NONZERO_ALLOCATION_ERR = (Int32)(-2)
const MAT_ROW_ORIENTED = (Int32)(-1)
const MAT_SYMMETRIC = (Int32)(1)
const MAT_STRUCTURALLY_SYMMETRIC = (Int32)(2)
const MAT_NEW_DIAGONALS = (Int32)(3)
const MAT_IGNORE_OFF_PROC_ENTRIES = (Int32)(4)
const MAT_USE_HASH_TABLE = (Int32)(5)
const MAT_KEEP_NONZERO_PATTERN = (Int32)(6)
const MAT_IGNORE_ZERO_ENTRIES = (Int32)(7)
const MAT_USE_INODES = (Int32)(8)
const MAT_HERMITIAN = (Int32)(9)
const MAT_SYMMETRY_ETERNAL = (Int32)(10)
const MAT_DUMMY = (Int32)(11)
const MAT_IGNORE_LOWER_TRIANGULAR = (Int32)(12)
const MAT_ERROR_LOWER_TRIANGULAR = (Int32)(13)
const MAT_GETROW_UPPERTRIANGULAR = (Int32)(14)
const MAT_SPD = (Int32)(15)
const MAT_NO_OFF_PROC_ZERO_ROWS = (Int32)(16)
const MAT_NO_OFF_PROC_ENTRIES = (Int32)(17)
const MAT_NEW_NONZERO_LOCATIONS = (Int32)(18)
const MAT_OPTION_MAX = (Int32)(19)

const Vec = Ptr{Nothing}
const VecType = Ptr{UInt8}
const VECMPI = "mpi"

const ScatterMode = UInt32
const SCATTER_FORWARD = (UInt32)(0)
const SCATTER_REVERSE = (UInt32)(1)
const SCATTER_FORWARD_LOCAL = (UInt32)(2)
const SCATTER_REVERSE_LOCAL = (UInt32)(3)
const SCATTER_LOCAL = (UInt32)(2)

const KSP = Ptr{Nothing}

const KSPType = Ptr{UInt8}
const KSPRICHARDSON = "richardson"
const KSPCHEBYSHEV = "chebyshev"
const KSPCG = "cg"
const KSPGROPPCG = "groppcg"
const KSPPIPECG = "pipecg"
const KSPCGNE = "cgne"
const KSPNASH = "nash"
const KSPSTCG = "stcg"
const KSPGLTR = "gltr"
const KSPFCG = "fcg"
const KSPGMRES = "gmres"
const KSPFGMRES = "fgmres"
const KSPLGMRES = "lgmres"
const KSPDGMRES = "dgmres"
const KSPPGMRES = "pgmres"
const KSPTCQMR = "tcqmr"
const KSPBCGS = "bcgs"
const KSPIBCGS = "ibcgs"
const KSPFBCGS = "fbcgs"
const KSPFBCGSR = "fbcgsr"
const KSPBCGSL = "bcgsl"
const KSPCGS = "cgs"
const KSPTFQMR = "tfqmr"
const KSPCR = "cr"
const KSPPIPECR = "pipecr"
const KSPLSQR = "lsqr"
const KSPPREONLY = "preonly"
const KSPQCG = "qcg"
const KSPBICG = "bicg"
const KSPMINRES = "minres"
const KSPSYMMLQ = "symmlq"
const KSPLCD = "lcd"
const KSPPYTHON = "python"
const KSPGCR = "gcr"

const ISLocalToGlobalMapping = Ptr{Nothing}

const VecScatter = Ptr{Nothing}

const ScatterMode = UInt32

const SCATTER_FORWARD = (UInt32)(0)
const SCATTER_REVERSE = (UInt32)(1)
const SCATTER_FORWARD_LOCAL = (UInt32)(2)
const SCATTER_REVERSE_LOCAL = (UInt32)(3)
const SCATTER_LOCAL = (UInt32)(2)

const NormType = UInt32

const NORM_1 = (UInt32)(0)
const NORM_2 = (UInt32)(1)
const NORM_FROBENIUS = (UInt32)(2)
const NORM_INFINITY = (UInt32)(3)
const NORM_1_AND_2 = (UInt32)(4)


text:
import Base.hash

abstract type AbstractBlock end

struct BlockHeader <: AbstractBlock
    version::UInt32
    prev_block::Vector{UInt8}
    merkle_root::Vector{UInt8}
    timestamp::UInt32
    bits::UInt32
    nonce::UInt32
    BlockHeader(version, prev_block, merkle_root, timestamp, bits, nonce) = new(version, prev_block, merkle_root, timestamp, bits, nonce)
end

function show(io::IO, z::BlockHeader)
    print(io, "Block\n--------\nVersion : ", z.version,
            "\nPrevious Block : ", bytes2hex(z.prev_block),
            "\nMerkle Root : ", bytes2hex(z.merkle_root),
            "\nTime Stamp : ", unix2datetime(z.timestamp),
            "\nBits : ", z.bits,
            "\nNonce : ", z.nonce)
end

"""
Takes a byte stream and parses a block. Returns a Block object
"""
function io2blockheader(s::IOBuffer)
    version = ltoh(reinterpret(UInt32, read(s, 4))[1])
    prev_block = read(s, 32)
    merkle_root = read(s, 32)
    timestamp = ltoh(reinterpret(UInt32, read(s, 4))[1])
    bits = ltoh(reinterpret(UInt32, read(s, 4))[1])
    nonce = ltoh(reinterpret(UInt32, read(s, 4))[1])
    return BlockHeader(version, prev_block, merkle_root, timestamp, bits, nonce)
end

"""
Returns the 80 byte block header
"""
function serialize(block::BlockHeader)
    result = Vector(reinterpret(UInt8, [htol(block.version)]))
    prev_block = copy(block.prev_block)
    append!(result, prev_block)
    append!(result, block.merkle_root)
    append!(result, Vector(reinterpret(UInt8, [htol(block.timestamp)])))
    append!(result, Vector(reinterpret(UInt8, [htol(block.bits)])))
    append!(result, Vector(reinterpret(UInt8, [htol(block.nonce)])))
end

"""
Returns the hash256 interpreted little endian of the block
"""
function hash(block::BlockHeader)
    s = serialize(block)
    h256 = hash256(s)
    return reverse!(h256)
end

"""
Human-readable hexadecimal of the block hash
"""
function id(block::BlockHeader)
    return bytes2hex(hash(block))
end

"""
Returns whether this block is signaling readiness for BIP9

    BIP9 is signalled if the top 3 bits are 001
    remember version is 32 bytes so right shift 29 (>> 29) and see if
    that is 001
"""
function bip9(block::BlockHeader)
    return block.version >> 29 == 0b001
end

"""
Returns whether this block is signaling readiness for BIP91

    BIP91 is signalled if the 5th bit from the right is 1
    shift 4 bits to the right and see if the last bit is 1
"""
function bip91(block::BlockHeader)
    return block.version >> 4 & 1 == 1
end

"""
Returns whether this block is signaling readiness for BIP141

    BIP91 is signalled if the 2nd bit from the right is 1
    shift 1 bit to the right and see if the last bit is 1
"""
function bip141(block::BlockHeader)
    return block.version >> 1 & 1 == 1
end

"""
Returns the proof-of-work target based on the bits

    last byte is exponent
    the first three bytes are the coefficient in little endian
    the formula is: coefficient * 256**(exponent-3)
"""
function target(block::BlockHeader)
    exponent = block.bits >> 24
    coefficient = block.bits & 0x00ffffff
    return coefficient * big(256)^(exponent - 3)
end

"""
Returns the block difficulty based on the bits

    difficulty is (target of lowest difficulty) / (block's target)
    lowest difficulty has bits that equal 0xffff001d
"""
function difficulty(block::BlockHeader)
    lowest = 0xffff * big(256)^(0x1d - 3)
    return div(lowest, target(block))
end

"""
Returns whether this block satisfies proof of work

    get the hash256 of the serialization of this block
    interpret this hash as a little-endian number
    return whether this integer is less than the target
"""
function check_pow(block::BlockHeader)
    block_hash = hash(block)
    proof = to_int(block_hash, little_endian=true)
    return proof < target(block)
end

"""
    validate_merkle_root(block::BlockHeader, hashes::Vector{Vector{UInt8}})
    -> Bool

Gets the merkle root of the hashes and checks that it's
the same as the merkle root of this block header.
"""
function validate_merkle_root(block::BlockHeader, hashes::Vector{Vector{UInt8}})
    hashes = [reverse!(copy(h)) for h in hashes]
    root = merkle_root(hashes)
    merkle_root(hashes) == block.merkle_root
end

struct Block <: AbstractBlock
    magic::UInt32
    size::UInt32
    header::BlockHeader
    tx_counter
    tx::Vector{UInt8}
end

# TODO Add Block serialization function


text:
# using ShaleDrillingModel
using JLD
using Plots
gr()

include("makeStateSpace.jl")

path_prices = joinpath(Pkg.dir("ShaleDrillingData"), "data", "price-transitions.jld")
@load path_prices pspace Πp

jldpath = joinpath(Pkg.dir("ShaleDrillingData"), "results", "results-full-exp-M2500.jld")
@load jldpath coefmat_full_OPG

# some primitives
β = (1.02 / 1.125) ^ (1./12.)  # real discount rate
royalty_rates = [0.125, 1./6., 0.1875, 0.20, 0.225, 0.25]
geology_types = 1:10
royalty_types = 1:length(royalty_rates)

# initial parameters
θfull = coefmat_full_OPG[2:9,1]
θt = θfull[1:end-1]
σ = θfull[end]

# problem sizes
Πz = Πp
nψ, nz =  51, size(Πz,1)
wp = well_problem(3,4,10)
zspace, ψspace = (pspace,1:2,), range(-4.0, stop=4.0, length=nψ)
prim = dcdp_primitives(:exp, β, wp, zspace, Πp, ψspace)
tmpv = dcdp_tmpvars(prim)
evs = dcdp_Emax(prim)

# check sizes of models
itype = (1./6., 8,)
ShaleDrillingModel.check_size(prim, evs)
solve_vf_all!(evs, tmpv, prim, θt, σ, itype, Val{false})

fillflows!(tmpv, prim, θt, σ, itype...)

ubv = tmpv.uin[:,:,1:2,1] .+ β .* evs.EV[:,:,end-1:end]
q2 = similar(ubv)
softmax3!(q2, tmpv.lse, tmpv.tmp, ubv)

EU = similar(evs.EV)
ubUfull = similar(tmpv.ubVfull)
q = similar(tmpv.ubVfull)
rin = similar(tmpv.uex)
rex = similar(tmpv.uex)
P = zeros(Float64,nz,nψ,dmax(wp)+1,length(wp))
pdct = makepdct(prim, Val{:u})
ShaleDrillingModel.fillflowrevs!(Val{:exp}, ShaleDrillingModel.flowrev, rin, rex, θt, σ, pdct, itype...)

size(P)
size(evs.EV)

prDrill_infill!(
    evs.EV, tmpv.uin, tmpv.uex, tmpv.ubVfull,
    EU, rin, rex, ubUfull,
    P, q,
    tmpv.lse, tmpv.tmp, tmpv.IminusTEVp,
    prim.wp, prim.Πz,
    tmpv.Πψtmp, prim.ψspace, σ, prim.β
    )

tmpc = dcdp_tmpcntrfact(prim)
prDrill_infill!(evs, tmpv, prim, EU, P, tmpc, θfull, (1./6., 8,) )
size(P)

prdrill = 1.0 .- P[:,:,1,:]
plot(ψspace, prdrill[17,:,:])

# let prim = dcdp_primitives(:exp, β, wp, zspace, Πz, ψspace),
#     tmpv = dcdp_tmpvars(prim),
#     sev = SharedEV([1,], prim, [1.0/8.0], 1:1),
#     zdims = length.(prim.zspace),
#     nψ = _nψ(prim),
#     nd = _nd(prim),
#     sEU = Array{Float64}(size(sev.EV)),
#     sP = Array{Float64}(zdims..., nψ, nd, length(prim.wp), length.(sev.itypes)...)
#
#     serial_solve_vf_all!(sev, tmpv, prim, θfull, Val{false})
#     serial_counterfact_all!(sev, tmpv, prim, sEU, sP, tmpc, θfull)
# end

# rmprocs(workers())
# pids = IN_SLURM ? addprocs_slurm(parse(Int, ENV["SLURM_CPUS_PER_TASK"])) : addprocs()

# @everywhere @show pwd()
# @everywhere using ShaleDrillingModel
#
# sev = SharedEV(pids, prim, royalty_rates, 1:1)
# @eval @everywhere begin
#     set_g_dcdp_primitives($prim)
#     set_g_dcdp_tmpvars($tmpv)
#     set_g_SharedEV($sev)
#     set_g_ItpSharedEV(ItpSharedEV(get_g_SharedEV(), get_g_dcdp_primitives(), $σv))
# end
# s = remotecall_fetch(get_g_SharedEV, pids[2])
# fetch(s)
#
# zero!(sev)
# @show parallel_solve_vf_all!(sev, vcat(θt,σv), Val{true})


text:
abstract type AbstractSherlog <: AbstractFloat end

struct Sherlog64{T<:AbstractFloat,i} <: AbstractSherlog    # T is the bitpattern type for logging, i the logbook id
    val::Float64                                           # the value is always Float64
end

# conversions back from Sherlog64
Base.Float64(x::Sherlog64) = x.val
Base.Float32(x::Sherlog64) = Float32(x.val)
Base.Float16(x::Sherlog64) = Float16(x.val)
Base.Int64(x::Sherlog64) = Int64(x.val)
Base.Int32(x::Sherlog64) = Int32(x.val)
Base.Int16(x::Sherlog64) = Int16(x.val)

# otherwise conversion to itself is ambigious
Sherlog64{T,i}(x::Sherlog64) where {T<:AbstractFloat,i} = Sherlog64{T,i}(x.val)
Sherlog64(x::Sherlog64{T,i}) where {T<:AbstractFloat,i} = Sherlog64{T,i}(x.val)

# generator functions from Float, Int and Bool
Sherlog64(x::AbstractFloat) = Sherlog64{Float16,1}(x)   # use Float16 and logbook #1 as standard
Sherlog64{T}(x::AbstractFloat) where T = Sherlog64{T,1}(x)      # use logbook #1 if no logbook provided
Sherlog64(x::Integer) = Sherlog64{Float16,1}(x)           # same for integers
Sherlog64{T}(x::Integer) where T = Sherlog64{T,1}(x)
Sherlog64{T}(x::Bool) where T = if x Sherlog64{T,1}(1) else Sherlog64{T,1}(0) end

Base.promote_rule(::Type{Int64},::Type{Sherlog64{T,i}}) where {T,i} = Sherlog64{T,i}
Base.promote_rule(::Type{Int32},::Type{Sherlog64{T,i}}) where {T,i} = Sherlog64{T,i}
Base.promote_rule(::Type{Int16},::Type{Sherlog64{T,i}}) where {T,i} = Sherlog64{T,i}

Base.promote_rule(::Type{Float64},::Type{Sherlog64{T,i}}) where {T,i} = Sherlog64{T,i}
Base.promote_rule(::Type{Float32},::Type{Sherlog64{T,i}}) where {T,i} = Sherlog64{T,i}
Base.promote_rule(::Type{Float16},::Type{Sherlog64{T,i}}) where {T,i} = Sherlog64{T,i}

Base.bitstring(x::Sherlog64) = bitstring(x.val)
Base.show(io::IO,x::Sherlog64) = print(io,"Sherlog64(",string(x.val),")")

Base.eps(::Type{Sherlog64{T,i}}) where {T,i} = eps(Float64)
Base.eps(x::Sherlog64) = eps(x.val)

Base.typemin(::Sherlog64{T,i}) where {T,i} = Sherlog64{T,i}(typemin(Float64))
Base.typemax(::Sherlog64{T,i}) where {T,i} = Sherlog64{T,i}(typemax(Float64))
Base.floatmin(::Sherlog64{T,i}) where {T,i} = Sherlog64{T,i}(floatmin(Float64))
Base.floatmax(::Sherlog64{T,i}) where {T,i} = Sherlog64{T,i}(floatmax(Float64))
Base.precision(::Sherlog64{T,i}) where {T,i} = Sherlog64{T,i}(precision(Float64))

function +(x::Sherlog64{T,i},y::Sherlog64{T,i}) where {T,i}
    r = x.val + y.val
    log_it(T,r,i)
    return Sherlog64{T,i}(r)
end

function -(x::Sherlog64{T,i},y::Sherlog64{T,i}) where {T,i}
    r = x.val - y.val
    log_it(T,r,i)
    return Sherlog64{T,i}(r)
end

function *(x::Sherlog64{T,i},y::Sherlog64{T,i}) where {T,i}
    r = x.val * y.val
    log_it(T,r,i)
    return Sherlog64{T,i}(r)
end

function /(x::Sherlog64{T,i},y::Sherlog64{T,i}) where {T,i}
    r = x.val / y.val
    log_it(T,r,i)
    return Sherlog64{T,i}(r)
end

function ^(x::Sherlog64{T,i},y::Sherlog64{T,i}) where {T,i}
    r = x.val ^ y.val
    log_it(T,r,i)
    return Sherlog64{T,i}(r)
end

for O in ( :(-), :(+),
           :sign,
           :prevfloat, :nextfloat,
           :round, :trunc, :ceil, :floor,
           :inv, :abs, :sqrt, :cbrt,
           :exp, :expm1, :exp2, :exp10,
           :log, :log1p, :log2, :log10,
           :rad2deg, :deg2rad, :mod2pi, :rem2pi,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
           :sinc, :sinpi, :cospi,
           :sind, :cosd, :tand, :cscd, :secd, :cotd,
           :asind, :acosd, :atand, :acscd, :asecd, :acotd
          )
    @eval begin
        function Base.$O(x::Sherlog64{T,i}) where {T,i}
            r = $O(x.val)
            log_it(T,r,i)
            return Sherlog64{T,i}(r)
        end
    end
end

for O in ( :(<), :(<=))
    @eval begin
        Base.$O(x::Sherlog64, y::Sherlog64) = $O(x.val, y.val)
    end
end


text:
module LightGBM

using Dates
import Base
import Libdl
import StatsBase
import Libdl

struct LibraryNotFoundError <: Exception
    msg::String
end

function find_library(library_name::String, custom_paths::Vector{String})

    # Search system filedirs first, returns empty string if not found
    output = Libdl.find_library(library_name)

    if output == ""
        # try specified paths
        @info("$(library_name) not found in system dirs, trying fallback")
        output = Libdl.find_library(library_name, custom_paths)
    else
        @info("$(library_name) found in system dirs!")
    end

    if output == ""
        throw(LibraryNotFoundError("$(library_name) not found. Please ensure this library is either in system dirs or the dedicated paths: $(custom_paths)"))
    end

    return output

end

const LGBM_library = find_library("lib_lightgbm", [@__DIR__])

include("wrapper.jl")
include("estimators.jl")
include("utils.jl")
include("fit.jl")
include("predict.jl")
include("cv.jl")
include("search_cv.jl")
include(joinpath(@__DIR__, "MLJInterface.jl"))

export fit!, predict, predict_classes, cv, search_cv, savemodel, loadmodel!
export LGBMEstimator, LGBMRegression, LGBMClassification

end # module LightGBM


text:
module ColorBitstring
export printbits, printlnbits

printred(io::IO, x) = print(io, "\x1b[31m" * x * "\x1b[0m")
printgreen(io::IO, x) = print(io, "\x1b[32m" * x * "\x1b[0m")
printblue(io::IO, x) = print(io, "\x1b[34m" * x * "\x1b[0m")

function printbits(io::IO, x::Float16)
       bts = bitstring(x)
       printred(io, bts[1:1])
       printgreen(io, bts[2:6])
       printblue(io, bts[7:end])
end

function printbits(io::IO, x::Float32)
   bts = bitstring(x)
    printred(io, bts[1:1])
    printgreen(io, bts[2:2+8-1])
    printblue(io, bts[2+8:end])
end

function printbits(io::IO, x::Float64)
    bts = bitstring(x)
    printred(io, bts[1:1])
    printgreen(io, bts[2:2+11-1])
    printblue(io, bts[2+11:end])
end

function printbits(io::IO, x::Signed)
    bts = bitstring(x)
    printred(io, bts[1:1])
    printblue(io, bts[2:end])
end

printbits(io::IO, x::Unsigned) = printblue(io, bitstring(x))

printbits(x) = printbits(stdout, x)
function printlnbits(x)
    printbits(x)
    println()
end
end