text:
module RNGTest

    using Random

    import Base: convert, getindex, pointer

    import Distributed: pmap

    include("../deps/deps.jl")

    swrite = Ref{Ptr{Ptr{Bool}}}()

    function __init__()
        swrite[] = cglobal(("swrite_Basic", libtestu01), Ptr{Bool})
        unsafe_store!(swrite[], 0, 1)
    end

    # WrappedRNG

    # TestU01 expects a standard function as input (created via
    # cfunction here). When one wants to test the random stream of
    # an AbstractRNG, a little more work has to be done before
    # passing it to TestU01. The type WrappedRNG wraps an
    # AbstractRNG into an object which knows which type of random
    # numbers to produce, and also whether the scalar or array API
    # should be used (this is useful when different algorithms are
    # used for each case, which results in different streams which
    # should be tested separately). Such an object can then be
    # passed to the Unif01 constructor.

    const TestableNumbers = Union{Int8, UInt8, Int16, UInt16, Int32, UInt32,
        Int64, UInt64, Int128, UInt128, Float16, Float32, Float64}

    mutable struct WrappedRNG{T<:TestableNumbers, RNG<:AbstractRNG}
        rng::RNG
        cache::Vector{T}
        fillarray::Bool
        vals::Vector{UInt32}
        idx::Int
    end

    function WrappedRNG(rng::RNG, ::Type{T}, fillarray = true, cache_size = 3*2^11 ÷ sizeof(T)) where {RNG, T}
        if T <: Integer && cache_size*sizeof(T) % sizeof(UInt32) != 0
            error("cache_size must be a multiple of $(Int(4/sizeof(T))) (for type $T)")
        elseif T === Float16 && cache_size % 6 != 0 || T === Float32 && cache_size % 3 != 0
            error("cache_size must be a multiple of 3 (resp. 6) for Float32 (resp. Float16)")
        end
        cache = Vector{T}(undef, cache_size)
        fillcache(WrappedRNG{T, RNG}(rng, cache, fillarray,
            unsafe_wrap(Array, convert(Ptr{UInt32}, pointer(cache)), sizeof(cache)÷sizeof(UInt32)),
            0)) # 0 is a dummy value, which will be set correctly by fillcache
    end

    # The ability to play with the cache size and the fillarray option is for advanced uses,
    # when one wants to test different code path of the particular RNG implementations, like
    # MersenneTwister from Base.
    # For now let's document only the type parameter in the wrap function:
    wrap(rng::AbstractRNG, ::Type{T}) where {T<:TestableNumbers} = WrappedRNG(rng, T)

    function fillcache(g::WrappedRNG{T}) where T
        if g.fillarray
            rand!(g.rng, g.cache)
        else
            for i = 1:length(g.cache)
                @inbounds g.cache[i] = rand(g.rng, T)
            end
        end
        g.idx = 0
        return g
    end

    function (g::WrappedRNG{T})() where T<:Integer
        g.idx+1 > length(g.vals) && fillcache(g)
        @inbounds return g.vals[g.idx+=1]
    end

    function (g::WrappedRNG{Float64})()
        g.idx+1 > length(g.cache) && fillcache(g)
        @inbounds return g.cache[g.idx+=1]
    end

    function (g::WrappedRNG{Float32})()
        g.idx+3 > length(g.cache) && fillcache(g)
        @inbounds begin
            f = Float64(g.cache[g.idx+1])
            # a Float32 has 24 bits of precision, but only 23 bit of entropy
            f += Float64(g.cache[g.idx+2])/exp2(23)
            f += Float64(g.cache[g.idx+=3])/exp2(46)
            return f % 1.0
        end
    end

    function (g::WrappedRNG{Float16})()
        g.idx+6 > length(g.cache) && fillcache(g)
        @inbounds begin
            f = Float64(g.cache[g.idx+1])
            # a Float16 has 10 bits of entropy
            f += Float64(g.cache[g.idx+2])/exp2(10)
            f += Float64(g.cache[g.idx+3])/exp2(20)
            f += Float64(g.cache[g.idx+4])/exp2(30)
            f += Float64(g.cache[g.idx+5])/exp2(40)
            f += Float64(g.cache[g.idx+=6])/exp2(50)
            return f % 1.0
        end
    end

    # RNGGenerator struct
    mutable struct Unif01
        ptr::Ptr{Cvoid}
        gentype::Type
        name::String
        function Unif01(f::Function, genname)
            for i in 1:100
                tmp = f()
                if typeof(tmp) != Float64 error("Function must return Float64") end
                if tmp < 0 || tmp > 1 error("Function must return values on [0,1]") end
            end
            cf = @cfunction($f, Float64, ())
            # TestU01 crashed if two unif01 object are generated. The only safe thing is to explicitly delete the object when used instead of using finalizers
            return new(ccall((:unif01_CreateExternGen01, libtestu01), Ptr{Cvoid}, (Ptr{UInt8}, Ptr{Cvoid}), genname, cf), Float64, genname)
        end

        function Unif01(g::WrappedRNG{T}, genname) where {T<:AbstractFloat}
            # we assume that g being created out of an AbstractRNG, it produces Floats in the interval [0,1)
            cf = @cfunction($g, Float64, ())
            return new(ccall((:unif01_CreateExternGen01, libtestu01), Ptr{Cvoid}, (Ptr{UInt8}, Ptr{Cvoid}), genname, cf), Float64, genname)
        end

        function Unif01(g::WrappedRNG{T}, genname) where {T<:Integer}
            @assert Cuint === UInt32
            cf = @cfunction($g, UInt32, ())
            return new(ccall((:unif01_CreateExternGenBits, libtestu01), Ptr{Cvoid}, (Ptr{UInt8}, Ptr{Cvoid}), genname, cf), UInt32, genname)
        end
    end
    function delete(obj::Unif01)
            if obj.gentype === Float64
                ccall((:unif01_DeleteExternGen01, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
            else
                ccall((:unif01_DeleteExternGenBits, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
            end
    end

    const RNGGenerator = Union{Function, WrappedRNG}

    # Result types

    ## gofw
    struct Gotw_TestArray
        data::Vector{Float64}
        Gotw_TestArray() = new(Vector{Float64}(undef, 11))
    end
    function getindex(obj::Gotw_TestArray, i::Symbol)
        i == :KSP && return obj.data[1]
        i == :KSM && return obj.data[2]
        i == :KS && return obj.data[3]
        i == :AD && return obj.data[4]
        i == :CM && return obj.data[5]
        i == :WG && return obj.data[6]
        i == :WU  && return obj.data[7]
        i == :Mean && return obj.data[8]
        i == :Var && return obj.data[9]
        i == :Cor && return obj.data[10]
        i == :Sum && return obj.data[11]
        throw(BoundsError())
    end

    for (t, sCreate, sDelete, sPval) in
        ((:ResPoisson, :sres_CreatePoisson, :sres_DeletePoisson, :getPValPoisson),
         (:MarsaRes, :smarsa_CreateRes, :smarsa_DeleteRes, :getPValSmarsa),
         (:KnuthRes2, :sknuth_CreateRes2, :sknuth_DeleteRes2, :getPValRes2))
        @eval begin
            # The types
            mutable struct $t
                ptr::Ptr{Cvoid}
                function $(t)()
                    res = new(ccall(($(string(sCreate)), libtestu01), Ptr{Cvoid}, (), ))
                    finalizer(delete, res)
                    return res
                end
            end
            # Finalizers
            function delete(obj::$t)
                ccall(($(string(sDelete)), libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
            end
            # pvalue extractors
            pvalue(obj::$t) = ccall(($(string(sPval)), libtestu01extractors), Float64, (Ptr{Cvoid},), obj.ptr)
        end
    end

    # Sres
    # Basic
    ## Type
    mutable struct ResBasic
        ptr::Ptr{Cvoid}
        function ResBasic()
            res = new(ccall((:sres_CreateBasic, libtestu01), Ptr{Cvoid}, (), ))
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::ResBasic)
        ccall((:sres_DeleteBasic, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::ResBasic)
        res = Gotw_TestArray()
        ccall((:getPValBasic, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return res
    end

    # Chi2
    ## Type
    mutable struct ResChi2
        ptr::Ptr{Cvoid}
        N::Int
        function ResChi2(N::Integer)
            res = new(ccall((:sres_CreateChi2, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::ResChi2)
        ccall((:sres_DeleteChi2, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::ResChi2)
        res = Gotw_TestArray()
        ccall((:getPValChi2, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return res[obj.N == 1 ? :Mean : :Sum]
    end

    # sknuth
    ## Type
    mutable struct KnuthRes1
        ptr::Ptr{Cvoid}
        N::Int
        function KnuthRes1(N::Integer)
            res = new(ccall((:sknuth_CreateRes1, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::KnuthRes1)
        ccall((:sknuth_DeleteRes1, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::KnuthRes1)
        chi = Gotw_TestArray()
        bas = Gotw_TestArray()
        ccall((:getPValRes1, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}), obj.ptr, chi.data, bas.data)
        return chi[obj.N == 1 ? :Mean : :Sum], bas[obj.N == 1 ? :Mean : :AD]
    end

    # smarsa
    ## Type
    mutable struct MarsaRes2
        ptr::Ptr{Cvoid}
        N::Int
        function MarsaRes2(N::Integer)
            res = new(ccall((:smarsa_CreateRes2, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::MarsaRes2)
        ccall((:smarsa_DeleteRes2, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::MarsaRes2)
        res = Gotw_TestArray()
        ccall((:getPValSmarsa2, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return obj.N == 1 ? res[:Mean] : res[:Sum]
    end

    # Walk
    ## Type
    mutable struct WalkRes
        ptr::Ptr{Cvoid}
        N::Int
        function WalkRes(N::Integer)
            res = new(ccall((:swalk_CreateRes, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::WalkRes)
        ccall((:swalk_DeleteRes, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::WalkRes)
        pH = Gotw_TestArray()
        pM = Gotw_TestArray()
        pJ = Gotw_TestArray()
        pR = Gotw_TestArray()
        pC = Gotw_TestArray()
        ccall((:getPVal_Walk, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}), obj.ptr, pH.data, pM.data, pJ.data, pR.data, pC.data)
        return obj.N == 1 ? (pH[:Mean], pM[:Mean], pJ[:Mean], pR[:Mean], pC[:Mean]) : (pH[:Sum], pM[:Sum], pJ[:Sum], pR[:Sum], pC[:Sum])
    end

    # Npairs
    struct Snpair_StatArray
        data::Vector{Float64}
        Snpair_StatArray() = new(Vector{Float64}(undef, 11))
    end
    function getindex(obj::Snpair_StatArray, i::Symbol)
        i == :NP && return obj.data[1]
        i == :NPS && return obj.data[2]
        i == :NPPR && return obj.data[3]
        i == :mNP && return obj.data[4]
        i == :mNP1 && return obj.data[5]
        i == :mNP1S && return obj.data[6]
        i == :mNP2  && return obj.data[7]
        i == :mNP2S && return obj.data[8]
        i == :NJumps && return obj.data[9]
        i == :BB && return obj.data[10]
        i == :BM && return obj.data[11]
        throw(BoundsError())
    end
    ## Type
    mutable struct NpairRes
        ptr::Ptr{Cvoid}
        N::Int
        function NpairRes(N::Integer)
            res = new(ccall((:snpair_CreateRes, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::NpairRes)
        ccall((:snpair_DeleteRes, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::NpairRes)
        res = Snpair_StatArray()
        ccall((:getPVal_Npairs, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return obj.N == 1 ? (res[:NP], res[:mNP]) : (res[:NP], res[:mNP1], res[:mNP2], res[:NJumps])
    end

    # scomp
    ## Type
    mutable struct CompRes
        ptr::Ptr{Cvoid}
        N::Int
        function CompRes(N::Integer)
            res = new(ccall((:scomp_CreateRes, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::CompRes)
        ccall((:scomp_DeleteRes, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::CompRes)
        num = Gotw_TestArray()
        size = Gotw_TestArray()
        ccall((:getPValScomp, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}), obj.ptr, num.data, size.data)
        return num[obj.N == 1 ? :Mean : :Sum], size[obj.N == 1 ? :Mean : :Sum]
    end

    # sspectral
    ## Type
    mutable struct SpectralRes
        ptr::Ptr{Cvoid}
        function SpectralRes()
            res = new(ccall((:sspectral_CreateRes, libtestu01), Ptr{Cvoid}, (), ))
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::SpectralRes)
        ccall((:sspectral_DeleteRes, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::SpectralRes)
        res = Gotw_TestArray()
        ccall((:getPValSspectral, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return res[:AD]
    end

    # sstring
    ## Type
    mutable struct StringRes
        ptr::Ptr{Cvoid}
        function StringRes()
            res = new(ccall((:sstring_CreateRes, libtestu01), Ptr{Cvoid}, (), ))
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::StringRes)
        ccall((:sstring_DeleteRes, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::StringRes)
        res = Gotw_TestArray()
        ccall((:getPValStringRes, libtestu01extractors), Cvoid, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return res
    end

    ## Type
    mutable struct StringRes2
        ptr::Ptr{Cvoid}
        N::Int
        function StringRes2(N::Integer)
            res = new(ccall((:sstring_CreateRes2, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::StringRes2)
        ccall((:sstring_DeleteRes2, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::StringRes2)
        res = Gotw_TestArray()
        res2 = ccall((:getPValStringRes2, libtestu01extractors), Float64, (Ptr{Cvoid}, Ptr{Float64}), obj.ptr, res.data)
        return res[obj.N == 1 ? :Mean : :Sum], res2
    end

    ## Type
    mutable struct StringRes3
        ptr::Ptr{Cvoid}
        N::Int
        function StringRes3(N::Integer)
            res = new(ccall((:sstring_CreateRes3, libtestu01), Ptr{Cvoid}, (), ), N)
            finalizer(delete, res)
            return res
        end
    end
    ## Finalizers
    function delete(obj::StringRes3)
        ccall((:sstring_DeleteRes3, libtestu01), Cvoid, (Ptr{Cvoid},), obj.ptr)
    end
    ## pvalue extractors
    function pvalue(obj::StringRes3)
        res1 = Gotw_TestArray()
        res2 = Gotw_TestArray()
        ccall((:getPValStringRes3, libtestu01extractors), Float64, (Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}), obj.ptr, res1.data, res2.data)
        return res1[obj.N == 1 ? :Mean : :Sum], res2[obj.N == 1 ? :Mean : :Sum]
    end

    #########
    # Tests #
    #########
    ## smarsa
    function smarsa_BirthdaySpacings(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer, t::Integer, p::Integer)
        unif01 = Unif01(gen, "")
        sres = ResPoisson()
        try
            ccall((:smarsa_BirthdaySpacings, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Clong, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, d, t, p)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function smarsa_GCD(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer)
        unif01 = Unif01(gen, "")
        sres = MarsaRes2(N)
        try
            ccall((:smarsa_GCD, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function smarsa_CollisionOver(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = MarsaRes()
        try
            ccall((:smarsa_CollisionOver, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Clong, Cint),
             unif01.ptr, sres.ptr, N, n,
             r, d, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function smarsa_Savir2(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, m::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:smarsa_Savir2, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Clong, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, m, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function smarsa_SerialOver(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:smarsa_SerialOver, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Clong, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, d, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function smarsa_MatrixRank(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, L::Integer, k::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:smarsa_MatrixRank, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s, L, k)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ## sknuth
    function sknuth_Collision(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = KnuthRes2()
        try
            ccall((:sknuth_Collision, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Clong, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, d, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_CollisionPermut(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = KnuthRes2()
        try
            ccall((:sknuth_CollisionPermut, libtestu01), Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                Cint, Cint),
            unif01.ptr, sres.ptr, N, n,
            r, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_CouponCollector(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:sknuth_CouponCollector, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, d)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_Gap(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, Alpha::Real, Beta::Real)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:sknuth_Gap, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Float64, Float64),
                unif01.ptr, sres.ptr, N, n,
                r, Alpha, Beta)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_MaxOft(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = KnuthRes1(N)
        try
            ccall((:sknuth_MaxOft, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, d, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_Permutation(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:sknuth_Permutation, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_Run(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, up::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:sknuth_Run, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, up)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sknuth_SimpPoker(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, d::Integer, k::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:sknuth_SimpPoker, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, d, k)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ## svaria
    function svaria_AppearanceSpacings(gen::RNGGenerator, N::Integer, Q::Integer, K::Integer, r::Integer, s::Integer, L::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:svaria_AppearanceSpacings, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Clong, Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, Q,
                K, r, s, L)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function svaria_SampleProd(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, t::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:svaria_SampleProd, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, t)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function svaria_SampleMean(gen::RNGGenerator, N::Integer, n::Integer, r::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:svaria_SampleMean, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint),
                unif01.ptr, sres.ptr, N, n,
                r)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function svaria_SampleCorr(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, k::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:svaria_SampleCorr, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, k)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function svaria_SumCollector(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, g::Float64)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:svaria_SumCollector, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cdouble),
                unif01.ptr, sres.ptr, N, n,
                r, g)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function svaria_WeightDistrib(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, k::Integer, alpha::Real, beta::Real)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:svaria_WeightDistrib, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Clong, Float64, Float64),
                unif01.ptr, sres.ptr, N, n,
                r, k, alpha, beta)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ## sstring
    function sstring_AutoCor(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, d::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:sstring_AutoCor, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s, d)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sstring_HammingCorr(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, L::Integer)
        unif01 = Unif01(gen, "")
        sres = StringRes()
        try
            ccall((:sstring_HammingCorr, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s, L)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sstring_HammingIndep(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, L::Integer, d::Integer)
        unif01 = Unif01(gen, "")
        sres = StringRes()
        try
            ccall((:sstring_HammingIndep, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s, L, d)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sstring_HammingWeight2(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, L::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:sstring_HammingWeight2, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Clong),
                unif01.ptr, sres.ptr, N, n,
                r, s, L)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sstring_LongestHeadRun(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, L::Integer)
        unif01 = Unif01(gen, "")
        sres = StringRes2(N)
        try
            ccall((:sstring_LongestHeadRun, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Clong),
                unif01.ptr, sres.ptr, N, n,
                r, s, L)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sstring_PeriodsInStrings(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer)
        unif01 = Unif01(gen, "")
        sres = ResChi2(N)
        try
            ccall((:sstring_PeriodsInStrings, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function sstring_Run(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer)
        unif01 = Unif01(gen, "")
        sres = StringRes3(N)
        try
            ccall((:sstring_Run, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ## swalk
    function swalk_RandomWalk1(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer, L0::Integer, L1::Integer)
        unif01 = Unif01(gen, "")
        sres = WalkRes(N)
        try
            ccall((:swalk_RandomWalk1, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Clong, Clong),
                unif01.ptr, sres.ptr, N, n,
                r, s, L0, L1)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ## snpair
    function snpair_ClosePairs(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, t::Integer, p::Integer, m::Integer)
        unif01 = Unif01(gen, "")
        sres = NpairRes(N)
        try
            ccall((:snpair_ClosePairs, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint, Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, t, p, m)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ## scomp
    function scomp_LempelZiv(gen::RNGGenerator, N::Integer, k::Integer, r::Integer, s::Integer)
        unif01 = Unif01(gen, "")
        sres = ResBasic()
        try
            ccall((:scomp_LempelZiv, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Cint,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, k,
                r, s)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end
    function scomp_LinearComp(gen::RNGGenerator, N::Integer, n::Integer, r::Integer, s::Integer)
        unif01 = Unif01(gen, "")
        sres = CompRes(N)
        try
            ccall((:scomp_LinearComp, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, n,
                r, s)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    # sspectral
    function sspectral_Fourier3(gen::RNGGenerator, N::Integer, k::Integer, r::Integer, s::Integer)
        unif01 = Unif01(gen, "")
        sres = SpectralRes()
        try
            ccall((:sspectral_Fourier3, libtestu01), Cvoid,
                (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Cint,
                    Cint, Cint),
                unif01.ptr, sres.ptr, N, k,
                r, s)
        finally
            delete(unif01)
        end
        return pvalue(sres)
    end

    ##################
    # Test Batteries #
    ##################
    for (snm, fnm) in ((:SmallCrush, :smallcrushTestU01), (:Crush, :crushTestU01), (:BigCrush, :bigcrushTestU01), (:pseudoDIEHARD, :diehardTestU01), (:FIPS_140_2, :fips_140_2TestU01))
        @eval begin
            function $(fnm)(f::RNGGenerator, fname::String)
                unif01 = Unif01(f, fname)
                try
                    ccall(($(string("bbattery_", snm)), libtestu01),
                        Cvoid, (Ptr{Cvoid},), unif01.ptr)
                finally
                    delete(unif01)
                end
            end
            $(fnm)(f::RNGGenerator) = $(fnm)(f::RNGGenerator, "")
        end
    end
     function smallcrushJulia(f::RNGGenerator)
         # @everywhere f = ()->g()
         testnames = [g->smarsa_BirthdaySpacings(g, 1, 5000000, 0, 1073741824, 2, 1),
                      g->sknuth_Collision(g, 1, 5000000, 0, 65536, 2),
                      g->sknuth_Gap(g, 1, 200000, 22, 0.0, .00390625),
                      g->sknuth_SimpPoker(g, 1, 400000, 24, 64, 64),
                      g->sknuth_CouponCollector(g, 1, 500000, 26, 16),
                      g->sknuth_MaxOft(g, 1, 2000000, 0, 100000, 6),
                      g->svaria_WeightDistrib(g, 1, 200000, 27, 256, 0.0, 0.125),
                      g->smarsa_MatrixRank(g, 1, 20000, 20, 10, 60, 60),
                      g->sstring_HammingIndep(g, 1, 500000, 20, 10, 300, 0)[:Mean],
                      g->swalk_RandomWalk1(g, 1, 1000000, 0, 30, 150, 150)]
         return pmap(t->t(f), testnames)
     end
     function bigcrushJulia(f::RNGGenerator)
         testnames = [g->smarsa_SerialOver(g, 1, 10^9, 0, 2^8, 3)[:Mean],
                      g->smarsa_SerialOver(g, 1, 10^9, 22, 2^8, 3)[:Mean],
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 0, 2^21, 2),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 9, 2^21, 2),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 0, 2^14, 3),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 16, 2^14, 3),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 0, 64, 7),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 24, 64, 7),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 0, 8, 14),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 27, 8, 14),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 0, 4, 21),
                      g->smarsa_CollisionOver(g, 30, 2*10^7, 28, 4, 21),
                      g->smarsa_BirthdaySpacings(g, 100, 10^7, 0, 2^31, 2, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 2*10^7, 0, 2^21, 3, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 3*10^7, 14, 2^16, 4, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 2*10^7, 0, 2^9, 7, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 2*10^7, 7, 2^9, 7, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 3*10^7, 14, 2^8, 8, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 3*10^7, 22, 2^8, 8, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 3*10^7, 0, 2^4, 16, 1),
                      g->smarsa_BirthdaySpacings(g, 20, 3*10^7, 26, 2^4, 16, 1),
                      g->snpair_ClosePairs(g, 30, 6*10^6, 0, 3, 0, 30),
                      g->snpair_ClosePairs(g, 20, 4*10^6, 0, 5, 0, 30),
                      g->snpair_ClosePairs(g, 10, 3*10^6, 0, 9, 0, 30),
                      g->snpair_ClosePairs(g, 5, 2*10^6, 0, 16, 0, 30),
                      g->sknuth_SimpPoker(g, 1, 4*10^8, 0, 8, 8),
                      g->sknuth_SimpPoker(g, 1, 4*10^8, 27, 8, 8),
                      g->sknuth_SimpPoker(g, 1, 10^8, 0, 32, 32),
                      g->sknuth_SimpPoker(g, 1, 10^8, 25, 32, 32),
                      g->sknuth_CouponCollector(g, 1, 2*10^8, 0, 8),
                      g->sknuth_CouponCollector(g, 1, 2*10^8, 10, 8),
                      g->sknuth_CouponCollector(g, 1, 2*10^8, 20, 8),
                      g->sknuth_CouponCollector(g, 1, 2*10^8, 27, 8),
                      g->sknuth_Gap(g, 1, 5*10^8, 0, 0.0, 1/16),
                      g->sknuth_Gap(g, 1, 3*10^8, 25, 0.0, 1/32),
                      g->sknuth_Gap(g, 1, 10^8, 0, 0.0, 1/128),
                      g->sknuth_Gap(g, 1, 10^7, 20, 0.0, 1/1024),
                      g->sknuth_Run(g, 5, 10^9, 0, 0),
                      g->sknuth_Run(g, 5, 10^9, 15, 1),
                      g->sknuth_Permutation(g, 1, 10^9, 0, 3),
                      g->sknuth_Permutation(g, 1, 10^9, 0, 5),
                      g->sknuth_Permutation(g, 1, 5*10^8, 0, 7),
                      g->sknuth_Permutation(g, 1, 5*10^8, 10, 10),
                      g->sknuth_CollisionPermut(g, 20, 2*10^7, 0, 14),
                      g->sknuth_CollisionPermut(g, 20, 2*10^7, 10, 14),
                      g->sknuth_MaxOft(g, 40, 10^7, 0, 10^5, 8),
                      g->sknuth_MaxOft(g, 30, 10^7, 0, 10^5, 16),
                      g->sknuth_MaxOft(g, 20, 10^7, 0, 10^5, 24),
                      g->sknuth_MaxOft(g, 20, 10^7, 0, 10^5, 32),
                      g->svaria_SampleProd(g, 40, 10^7, 0, 8)[:AD],
                      g->svaria_SampleProd(g, 20, 10^7, 0, 16)[:AD],
                      g->svaria_SampleProd(g, 20, 10^7, 0, 24)[:AD],
                      g->svaria_SampleMean(g, 2*10^7, 30, 0)[:AD],
                      g->svaria_SampleMean(g, 2*10^7, 30, 10)[:AD],
                      g->svaria_SampleCorr(g, 1, 2*10^9, 0, 1)[:Mean],
                      g->svaria_SampleCorr(g, 1, 2*10^9, 0, 2)[:Mean],
                      g->svaria_AppearanceSpacings(g, 1, 10^7, 10^9, 0, 3, 15)[:Mean],
                      g->svaria_AppearanceSpacings(g, 1, 10^7, 10^9, 27, 3, 15)[:Mean],
                      g->svaria_WeightDistrib(g, 1, 2*10^7, 0, 256, 0.0, 1/4),
                      g->svaria_WeightDistrib(g, 1, 2*10^7, 20, 256, 0.0, 1/4),
                      g->svaria_WeightDistrib(g, 1, 2*10^7, 28, 256, 0.0, 1/4),
                      g->svaria_WeightDistrib(g, 1, 2*10^7, 0, 256, 0.0, 1/16),
                      g->svaria_WeightDistrib(g, 1, 2*10^7, 10, 256, 0.0, 1/16),
                      g->svaria_WeightDistrib(g, 1, 2*10^7, 26, 256, 0.0, 1/16),
                      g->svaria_SumCollector(g, 1, 5*10^8, 0, 10.0),
                      g->smarsa_MatrixRank(g, 10, 10^6, 0, 5, 30, 30),
                      g->smarsa_MatrixRank(g, 10, 10^6, 25, 5, 30, 30),
                      g->smarsa_MatrixRank(g, 1, 5000, 0, 4, 1000, 1000),
                      g->smarsa_MatrixRank(g, 1, 5000, 26, 4, 1000, 1000),
                      g->smarsa_MatrixRank(g, 1, 80, 15, 15, 5000, 5000),
                      g->smarsa_MatrixRank(g, 1, 80, 0, 30, 5000, 5000),
                      g->smarsa_Savir2(g, 10, 10^7, 10, 2^30, 30),
                      g->smarsa_GCD(g, 10, 5*10^7, 0, 30),
                      g->swalk_RandomWalk1(g, 1, 10^8, 0, 5, 50, 50),
                      g->swalk_RandomWalk1(g, 1, 10^8, 25, 5, 50, 50),
                      g->swalk_RandomWalk1(g, 1, 10^7, 0, 10, 1000, 1000),
                      g->swalk_RandomWalk1(g, 1, 10^7, 20, 10, 1000, 1000),
                      g->swalk_RandomWalk1(g, 1, 10^6, 0, 15, 1000, 1000),
                      g->swalk_RandomWalk1(g, 1, 10^6, 15, 15, 10000, 10000),
                      g->scomp_LinearComp(g, 1, 400000, 0, 1),
                      g->scomp_LinearComp(g, 1, 400000, 29, 1),
                      g->scomp_LempelZiv(g, 10, 27, 0, 30)[:Sum],
                      g->scomp_LempelZiv(g, 10, 27, 15, 15)[:Sum],
                      g->sspectral_Fourier3(g, 100000, 14, 0, 3),
                      g->sspectral_Fourier3(g, 100000, 14, 27, 3),
                      g->sstring_LongestHeadRun(g, 1, 1000, 0, 3, 10^7),
                      g->sstring_LongestHeadRun(g, 1, 1000, 27, 3, 10^7),
                      g->sstring_PeriodsInStrings(g, 10, 5*10^8, 0, 10),
                      g->sstring_PeriodsInStrings(g, 10, 5*10^8, 20, 10),
                      g->sstring_HammingWeight2(g, 10, 10^9, 0, 3, 10^6)[:Sum],
                      g->sstring_HammingWeight2(g, 10, 10^9, 27, 3, 10^6)[:Sum],
                      g->sstring_HammingCorr(g, 1, 10^9, 10, 10, 30)[:Mean],
                      g->sstring_HammingCorr(g, 1, 10^8, 10, 10, 300)[:Mean],
                      g->sstring_HammingCorr(g, 1, 10^8, 10, 10, 1200)[:Mean],
                      g->sstring_HammingIndep(g, 10, 3*10^7, 0, 3, 30, 0)[:Sum],
                      g->sstring_HammingIndep(g, 10, 3*10^7, 27, 3, 30, 0)[:Sum],
                      g->sstring_HammingIndep(g, 1, 3*10^7, 0, 4, 300, 0)[:Mean],
                      g->sstring_HammingIndep(g, 1, 3*10^7, 26, 4, 300, 0)[:Mean],
                      g->sstring_HammingIndep(g, 1, 10^7, 0, 5, 1200, 0)[:Mean],
                      g->sstring_HammingIndep(g, 1, 10^7, 25, 5, 1200, 0)[:Mean],
                      g->sstring_Run(g, 1, 2*10^9, 0, 3),
                      g->sstring_Run(g, 1, 2*10^9, 27, 3),
                      g->sstring_AutoCor(g, 10, 10^9, 0, 3, 1)[:Sum],
                      g->sstring_AutoCor(g, 10, 10^9, 0, 3, 3)[:Sum],
                      g->sstring_AutoCor(g, 10, 10^9, 27, 3, 1)[:Sum],
                      g->sstring_AutoCor(g, 10, 10^9, 27, 3, 3)[:Sum]]
         return pmap(t->t(f), testnames)
     end
end


text:
using DynamicGrids, DimensionalData, Dates, Test
using DynamicGrids: isshowable, frameindex, storeframe!, SimData, stoppedframe

# Mostly outputs are tested in integration.jl
@testset "Output construction" begin
    init = [10.0 11.0
            0.0   5.0]

    output = ArrayOutput(init; tspan=Date(2001):Year(1):Date(2010))
    ruleset = Ruleset(Life())

    @test length(output) == 10
    @test size(output) == (10,)
    @test step(output) == Year(1)
    @test stoppedframe(output) == 10
    @test timestep(output) == Year(1)
    @test_throws ArgumentError DynamicGrids.ruleset(output)
    @test frameindex(output, 5) == 5
    @test isshowable(output, 5) == false
    @test output[1] == output[Ti(1)] == init

    @testset "DimensionalData interface" begin
        @test output isa AbstractDimArray{<:Array,1,<:Tuple{<:Ti}}
        @test dims(output) isa Tuple{<:Ti}
        @test DimensionalData.name(output) == NoName()
        @test metadata(output) == NoMetadata()
        da = output[Ti(Between(Date(2002), Date(2003)))]
        @test da isa DimArray{<:Array,1,<:Tuple{<:Ti}}
        @test index(da) == (Date(2002):Year(1):Date(2002),)
    end

    @testset "errors" begin
        @test_throws UndefKeywordError ArrayOutput(ones(5, 5))
        @test_throws ArgumentError ArrayOutput((a=ones(5, 5), b=ones(4, 4)); tspan=1:10)
        @test_throws ArgumentError ArrayOutput(ones(5, 5); mask=ones(2, 2), tspan=1:10)
    end
end

ext: jl


text:
@doc raw"""
    hankel_image(v, channel_size, kernel_size)

Creates a Hankel array from `v` with channel and kernel dimensions
`channel_size` and `kernel_size`, respectively. The returned array `V` satisfies:

```math
V_{c,j,k,n} = V_{c,j+k-1,n}
```

where `c,j,k,n` are multi-indices, with `c` traversing channels,
`j` traversing the kernel, `k` the convolution output, and `n` batches.
"""
function hankel_image(
    img::AbstractArray, channel_sz::NTuple{C,Int}, kernel_sz::NTuple{N,Int},
    stride::NTuple{N,Int}, pad::NTuple{N2,Int}, dilation::NTuple{N,Int}
) where {C,N,N2}
    input_sz = size(img)[(C + 1):(C + N)]
    batch_sz = size(img)[(C + N + 1):end]
    @assert size(img) == (channel_sz..., input_sz..., batch_sz...)
    out_sz = output_size(kernel_sz, input_sz; stride, pad, dilation)
    pad_lo = CartesianIndex(splitpad(pad).lo)
    X = zeros(eltype(img), channel_sz..., kernel_sz..., out_sz..., batch_sz...)
    for k in CartesianIndices(out_sz), j in CartesianIndices(kernel_sz)
        i = out2in(j, k; dilation, stride) - pad_lo
        all(1 .≤ Tuple(i) .≤ input_sz) || continue
        for b in CartesianIndices(batch_sz), c in CartesianIndices(channel_sz)
            X[c,j,k,b] = img[c,i,b]
        end
    end
    return X
end

function hankel_image(
    img::AbstractArray, channel_sz::NTuple{C,Int}, kernel_sz::NTuple{N,Int};
    stride = 1, pad = 0, dilation = 1
) where {C,N}
    return hankel_image(
        img, channel_sz, kernel_sz,
        expand_tuple(Val(N), stride),
        parsepad(kernel_sz, dilation, pad),
        expand_tuple(Val(N), dilation)
    )
end

function hankel_image(img::AbstractArray, C::Int, kernel_size::TupleN{Int}; kwargs...)
    channel_size = size(img)[1:C]
    return hankel_image(img, channel_size, kernel_size; kwargs...)
end

@doc raw"""
    hankel_weight(w, channel_size, input_size)

Creates a Hankel array from `v` with channel and kernel dimensions
`channel_size` and `kernel_size`, respectively. The returned array `V` satisfies:

```math
V_{c,j,k,n} = V_{c,j+k-1,n}
```

where `c,j,k,n` are multi-indices, with `c` traversing channels,
`j` traversing the kernel, `k` the convolution output, and `n` batches.
"""
function hankel_weight(
    w::AbstractArray, channel_sz::NTuple{C,Int}, in_sz::NTuple{N,Int};
    stride = 1, pad = 0, dilation = 1
) where {C,N}
    kernel_sz = size(w)[(C + 1):(C + N)]
    return hankel_weight(
        w, channel_sz, in_sz,
        expand_tuple(Val(N), stride),
        parsepad(kernel_sz, dilation, pad),
        expand_tuple(Val(N), dilation)
    )
end

function hankel_weight(
    w::AbstractArray, channel_sz::NTuple{C,Int}, in_sz::NTuple{N,Int},
    stride::NTuple{N,Int}, pad::NTuple{N2,Int}, dilation::NTuple{N,Int}
) where {C,N,N2}
    kernel_sz = size(w)[(C + 1):(C + N)]
    hidden_sz = size(w)[(C + N + 1):end]
    @assert size(w) == (channel_sz..., kernel_sz..., hidden_sz...)
    out_sz = output_size(kernel_sz, in_sz; stride, pad, dilation)
    pad_lo = CartesianIndex(splitpad(pad).lo)
    A = zeros(eltype(w), channel_sz..., in_sz..., hidden_sz..., out_sz...)
    for k in CartesianIndices(out_sz), j in CartesianIndices(kernel_sz)
        i = out2in(j, k; dilation, stride) - pad_lo
        all(1 .≤ Tuple(i) .≤ in_sz) || continue
        for μ in CartesianIndices(hidden_sz), c in CartesianIndices(channel_sz)
            A[c,i,μ,k] = w[c,j,μ]
        end
    end
    return A
end

function hankel_weight(w::AbstractArray, C::Int, in_sz::TupleN{Int}; kwargs...)
    channel_size = size(w)[1:C]
    return hankel_weight(w, channel_size, in_sz; kwargs...)
end

"""
    hankel(convrbm, input_size)

Returns a dense `RBM` equivalent to `convrbm` with replicated weights.
"""
function hankel(rbm::ConvRBM, in_size::NTuple{N,Int}) where {N}
    @assert kernel_ndims(rbm) == N
    @assert rbm.groups == 1 # groups > 1 not supported
    w = hankel_weight(rbm.w, channel_size(rbm), in_size; rbm.stride, rbm.pad, rbm.dilation)
    out_size = output_size(rbm, in_size)
    vis = replicate(visible(rbm), in_size...)
    hid = replicate(hidden(rbm), out_size...)
    return RBM(vis, hid, w)
end

"""
    replicate(layer, n...)

Returns a new layer of size `(size(layer)..., n...)` by repeating the original `layer`
along the new dimensions.
"""
replicate(l::AbstractLayer, n::Int...) = repeat(l, map(Returns(1), size(l))..., n...)


text:
using Test

@testset "NDTensors.jl" begin
  @testset "$filename" for filename in (
    "linearalgebra.jl",
    "dense.jl",
    "blocksparse.jl",
    "diag.jl",
    "itensors.jl",
  )
    println("Running $filename")
    include(filename)
  end
end

nothing


text:
export PolarForm, bounds, powerflow
export State, Control, Parameters, NumberOfState, NumberOfControl

"""
    AbstractStructure

The user may specify a mapping to the single input vector `x` for AD.

"""
abstract type AbstractStructure end

"""
    AbstractFormulation

Second layer of the package, implementing the interface between
the first layer (the topology of the network) and the
third layer (implementing the callbacks for the optimization solver).

"""
abstract type AbstractFormulation end

"""
    AbstractFormAttribute

Attributes attached to an `AbstractFormulation`.
"""
abstract type AbstractFormAttribute end

"Number of states attached to a particular formulation."
struct NumberOfState <: AbstractFormAttribute end

"Number of controls attached to a particular formulation."
struct NumberOfControl <: AbstractFormAttribute end

"""
    AbstractVariable

Variables corresponding to a particular formulation.
"""
abstract type AbstractVariable end

"""
    State <: AbstractVariable

All variables ``x`` depending on the variables `Control` ``u`` through
a non-linear equation ``g(x, u) = 0``.

"""
struct State <: AbstractVariable end

"""
    Control <: AbstractVariable

Independent variables ``u`` used in the reduced-space
formulation.

"""
struct Control <: AbstractVariable end

"""
    PhysicalState <: AbstractVariable

All physical variables describing the current physical state
of the underlying network.

`PhysicalState` variables are encoded in a `AbstractNetworkBuffer`,
storing all the physical values needed to describe the current
state of the network.

"""
struct PhysicalState <: AbstractVariable end

# Templates
"""
    get(form::AbstractFormulation, attr::AbstractFormAttribute)

Return value of attribute `attr` attached to the particular
formulation `form`.

## Examples

```julia
get(form, NumberOfState())
get(form, NumberOfControl())

```
"""
function get end

"""
    setvalues!(form::AbstractFormulation, attr::PS.AbstractNetworkAttribute, values)

Update inplace the attribute's values specified by `attr`.

## Examples

```julia
setvalues!(form, ActiveLoad(), new_ploads)
setvalues!(form, ReactiveLoad(), new_qloads)

```
"""
function setvalues! end

"""
    bounds(form::AbstractFormulation, var::AbstractVariable)

Return the bounds attached to the variable `var`.

    bounds(form::AbstractFormulation, func::Function)

Return a tuple of vectors `(lb, ub)` specifying the admissible range
of the constraints specified by the function `cons_func`.

## Examples

```julia
u_min, u_max = bounds(form, Control())
h_min, h_max = bounds(form, reactive_power_constraints)

```
"""
function bounds end

"""
    initial(form::AbstractFormulation, var::AbstractVariable)

Return an initial position for the variable `var`.

## Examples

```julia
u₀ = initial(form, Control())
x₀ = initial(form, State())

```
"""
function initial end

"""
    powerflow(form::AbstractFormulation,
              algo::AbstractNonLinearSolver;
              kwargs...)

    powerflow(form::AbstractFormulation,
              jacobian::AutoDiff.Jacobian,
              buffer::AbstractNetworkBuffer,
              algo::AbstractNonLinearSolver;
              kwargs...) where VT <: AbstractVector

Solve the power flow equations ``g(x, u) = 0`` w.r.t. the state ``x``,
using the algorithm specified in `algo` ([`NewtonRaphson`](@ref) by default).
The initial state ``x`` is specified inside
`buffer`. The object `buffer` is modified inplace in the function.

The algorithm stops when a tolerance `tol` or a maximum number of
iterations `maxiter` is reached (these parameters being specified
in the object `algo`).

## Notes
If only the arguments `form` and `algo` are specified to the function,
then the Jacobian `jacobian` and the cache `buffer` are inferred
from the object `form`.

## Arguments

* `form::AbstractFormulation`: formulation of the power flow equation
* `jacobian::AutoDiff.Jacobian`: Jacobian
* `buffer::AbstractNetworkBuffer`: buffer storing current state `x` and control `u`
* `algo::AbstractNonLinearSolver`: non-linear solver. Currently only `NewtonRaphson` is being implemented.

## Optional arguments

* `linear_solver::AbstractLinearSolver` (default `DirectSolver()`): solver to solve the linear systems ``J x = y`` arising at each iteration of the Newton-Raphson algorithm.

"""
function powerflow end

# Cost function
"""
    cost_production(form::AbstractFormulation, buffer::AbstractNetworkBuffer)::Float64

Get operational cost.
"""
function cost_production end

# Generic constraints

"""
    voltage_magnitude_constraints(form::AbstractFormulation, cons::AbstractVector, buffer::AbstractNetworkBuffer)

Bounds the voltage magnitudes at PQ nodes:
```math
v_{pq}^♭ ≤ v_{pq} ≤ v_{pq}^♯ .
```
The result is stored inplace, inside `cons`.

## Note
The constraints on the voltage magnitudes at PV nodes ``v_{pv}``
are taken into account when bounding the control ``u``.
"""
function voltage_magnitude_constraints end

"""
    active_power_constraints(form::AbstractFormulation, cons::AbstractVector, buffer::AbstractNetworkBuffer)

Evaluate the constraints on the **active power production** at the generators
that are not already taken into account in the bound constraints.
```math
p_g^♭ ≤ p_g ≤ p_g^♯  .
```

The result is stored inplace, inside the vector `cons`.
"""
function active_power_constraints end

"""
    active_power_generation(form::AbstractFormulation, pg::AbstractVector, buffer::AbstractNetworkBuffer)

Evaluate the **active power production** of all generators.
Used to evaluate the operational cost.
The result is stored inplace, inside the vector `pg`.
"""
function active_power_generation end

"""
    reactive_power_constraints(form::AbstractFormulation, cons::AbstractVector, buffer::AbstractNetworkBuffer)

Evaluate the constraints on the **reactive power production** at the generators:
```math
q_g^♭ ≤ q_g ≤ q_g^♯  .
```
The result is stored inplace, inside the vector `cons`.
"""
function reactive_power_constraints end

"""
    flow_constraints(form::AbstractFormulation, cons::AbstractVector, buffer::AbstractNetworkBuffer)

Evaluate the thermal limit constraints porting on the lines of the network.
The result is stored inplace, inside the vector `cons`.
"""
function flow_constraints end

@doc raw"""
    power_balance(form::AbstractFormulation, cons::AbstractVector, buffer::AbstractNetworkBuffer)

Evaluate the power balance in the network:
```math
g(x, u) = 0 ,
```
corresponding to the balance equations
```math
\begin{aligned}
    p_i &= v_i \sum_{j}^{n} v_j (g_{ij}\cos{(\theta_i - \theta_j)} + b_{ij}\sin{(\theta_i - \theta_j})) \,, &
    ∀ i ∈ \{PV, PQ\} \\
    q_i &= v_i \sum_{j}^{n} v_j (g_{ij}\sin{(\theta_i - \theta_j)} - b_{ij}\cos{(\theta_i - \theta_j})) \,. &
    ∀ i ∈ \{PQ\} \\
\end{aligned}
```

The result is stored inplace, inside the vector `cons`.
"""
function power_balance end

# Interface for the constraints
"""
    size_constraint(cons_func::Function)::Bool
Return whether the function `cons_func` is a supported constraint
in the powerflow model.
"""
function is_constraint end

"""
    size_constraint(form::AbstractFormulation, cons_func::Function)::Int

Get number of constraints specified by the function `cons_func`
in the formulation `form`.
"""
function size_constraint end

"""
    adjoint!(form::AbstractFormulation, pbm::AutoDiff.TapeMemory, adj_h, h, buffer)

Return the adjoint w.r.t. the variables of the network (voltage magnitudes
and angles, power injection) for the constraint stored inside the [`AutoDiff.TapeMemory`](@ref)
object `pbm`. The results are stored directly inside the stack stored
inside `pbm`.
"""
function adjoint! end

"""
    jacobian_transpose_product!(form::AbstractFormulation, pbm::AutoDiff.TapeMemory, buffer, v)

Return the two transpose-Jacobian vector product ``(Jᵤ^⊤ v, Jₓ^⊤ v)``  w.r.t. the
control ``u`` and the state ``x``. Store the two resulting vectors directly inside
the [`AutoDiff.TapeMemory`](@ref) `pbm`.

"""
function jacobian_transpose_product! end

"""
    matpower_jacobian(form::AbstractFormulation, X::Union{State,Control}, cons_func::Function, V::Vector{Complex})
    matpower_jacobian(form::AbstractFormulation, X::Union{State,Control}, cons_func::Function, buffer::AbstractNetworkBuffer)

For the constraint `cons_func`, return the expression of the Jacobian ``J``
w.r.t. the state or the control (depending on the argument `X`),
as given by MATPOWER.
"""
function matpower_jacobian end

@doc raw"""
    matpower_hessian(form::AbstractFormulation, cons_func::Function, buffer::AbstractNetworkBuffer, λ::AbstractVector)

For constraint `cons_func`, return the three matrices ``(λ^⊤ H_{xx},
λ^⊤ H_{xu},λ^⊤ H_{uu})`` storing the product of the Hessian tensor ``H`` with the vector ``\lambda``.
The expressions of the Hessian matrices are given by MATPOWER.

"""
function matpower_hessian end

"""
    jacobian_sparsity(form::AbstractFormulation, cons_func::Function, X::Union{State,Control})

For the constraint `cons_func`, return the sparsity pattern of the Jacobian ``J``
w.r.t. the state or the control (depending on the argument `X`).

"""
function jacobian_sparsity end


text:
using Test
import GLPK

include("tests.jl")

Tests.runtests(GLPK.Optimizer)


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestMutableArithmetics

using Test

import MutableArithmetics
const MA = MutableArithmetics

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

_zero(::Type{MOI.VariableIndex}, T::Type) = zero(MOI.ScalarAffineFunction{T})
_zero(x, ::Type) = zero(x)

function _promote_operation_test(op::Function, T, x::Type, y::Type)
    f() = MA.promote_operation(op, x, y)
    @test typeof(op(_zero(x, T), _zero(y, T))) == f()
    @test 0 == @allocated f()
end

function _test_promote_operation_allocation(T)
    AffType = MOI.ScalarAffineFunction{T}
    QuadType = MOI.ScalarQuadraticFunction{T}
    for op in [+, -, *]
        _promote_operation_test(op, T, T, MOI.VariableIndex)
        _promote_operation_test(op, T, MOI.VariableIndex, T)
        _promote_operation_test(op, T, T, AffType)
        _promote_operation_test(op, T, AffType, T)
        _promote_operation_test(op, T, T, QuadType)
        _promote_operation_test(op, T, QuadType, T)
        _promote_operation_test(op, T, MOI.VariableIndex, AffType)
        _promote_operation_test(op, T, AffType, MOI.VariableIndex)
        if op != *
            _promote_operation_test(op, T, MOI.VariableIndex, QuadType)
            _promote_operation_test(op, T, QuadType, MOI.VariableIndex)
            _promote_operation_test(op, T, AffType, QuadType)
            _promote_operation_test(op, T, QuadType, AffType)
        end
    end
    return
end

function test_promote_operation_allocation_Int()
    return _test_promote_operation_allocation(Int)
end
function test_promote_operation_allocation_Float64()
    return _test_promote_operation_allocation(Float64)
end
function test_promote_operation_allocation_Float32()
    return _test_promote_operation_allocation(Float32)
end

function _test_promote_operation(T)
    @test MA.promote_operation(*, MOI.VariableIndex, T) ==
          MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, T, MOI.VariableIndex) ==
          MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, MOI.ScalarAffineFunction{T}, T) ==
          MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, T, MOI.ScalarAffineFunction{T}) ==
          MOI.ScalarAffineFunction{T}
    @test MA.promote_operation(*, MOI.ScalarQuadraticFunction{T}, T) ==
          MOI.ScalarQuadraticFunction{T}
    @test MA.promote_operation(*, T, MOI.ScalarQuadraticFunction{T}) ==
          MOI.ScalarQuadraticFunction{T}
    return
end

test_promote_operation_Int() = _test_promote_operation(Int)
test_promote_operation_Float64() = _test_promote_operation(Float64)
test_promote_operation_Float32() = _test_promote_operation(Float32)

function _test_scaling(T)
    x = MOI.VariableIndex(1)
    @test T(3) == MA.scaling(T(0)x + T(3))
    f = T(2)x + T(3)
    err = InexactError(:convert, T, f)
    @test_throws err MA.scaling(f)
end

test_scaling_Float64() = _test_scaling(Float64)
test_scaling_Float32() = _test_scaling(Float32)

function test_unary_minus()
    for T in [Float64, Float32]
        x = MOI.VariableIndex(1)
        for f in [T(2)x + T(3), T(4) * x * x + T(2)x + T(3)]
            g = -f
            @test g ≈ MA.operate!(-, f)
            @test g ≈ f
            @test -g ≈ MA.operate!(-, f)
            @test -g ≈ f
        end
    end
    return
end

function _run_all_tests(T::Type, a, b, c, d, e, f, g)
    exclude = ["scalar_uniform_scaling", "cube"]
    @testset "Scalar" begin
        MA.Test.scalar_test(a, exclude = exclude)
        MA.Test.scalar_test(b, exclude = exclude)
        MA.Test.scalar_test(c, exclude = exclude)
        MA.Test.scalar_test(d, exclude = exclude)
        MA.Test.scalar_test(e, exclude = exclude)
        MA.Test.scalar_test(f, exclude = exclude)
        MA.Test.scalar_test(g, exclude = exclude)
    end
    if !(a isa MOI.ScalarQuadraticFunction)
        exclude = String[]
        if T <: Integer
            push!(exclude, "quadratic_division")
        end
        @testset "Quadratic" begin
            MA.Test.quadratic_test(a, b, c, d, exclude = exclude)
            MA.Test.quadratic_test(b, c, d, e, exclude = exclude)
            MA.Test.quadratic_test(c, d, e, f, exclude = exclude)
            MA.Test.quadratic_test(d, e, f, g, exclude = exclude)
            MA.Test.quadratic_test(e, f, g, a, exclude = exclude)
            MA.Test.quadratic_test(f, g, a, b, exclude = exclude)
            MA.Test.quadratic_test(g, a, b, c, exclude = exclude)
        end
    end
    exclude = String[]
    if a isa MOI.ScalarQuadraticFunction
        push!(exclude, "sparse_quadratic")
    end
    @testset "Sparse" begin
        MA.Test.sparse_test(a, b, [a b c; b c a; a b a], exclude = exclude)
    end
    exclude = String[]
    if a isa MOI.ScalarQuadraticFunction
        push!(exclude, "matrix_vector")
        push!(exclude, "dot")
    end
    @testset "Vector" begin
        MA.Test.array_test([a, b, c], exclude = exclude)
        MA.Test.array_test([b, c, d], exclude = exclude)
        MA.Test.array_test([c, d, e], exclude = exclude)
        MA.Test.array_test([d, e, a], exclude = exclude)
        MA.Test.array_test([e, a, b], exclude = exclude)
    end
    exclude = [
        "broadcast_division",
        "matrix_uniform_scaling",
        "symmetric_matrix_uniform_scaling",
        "matrix_vector_division",
    ]
    if a isa MOI.ScalarQuadraticFunction
        push!(exclude, "dot")
        push!(exclude, "broadcast_multiplication")
        push!(exclude, "sum_multiplication")
    end
    @testset "Matrix" begin
        MA.Test.array_test([a b; c d], exclude = exclude)
        MA.Test.array_test([c e; e d], exclude = exclude)
        MA.Test.array_test([a b c; b c a; a b a], exclude = exclude)
        MA.Test.array_test([d b c; d c e; e b a], exclude = exclude)
    end
    return
end

function test_VariableIndex()
    x = MOI.VariableIndex(1)
    a = 2x
    MA.Test.@test_rewrite(a + x)
    MA.Test.@test_rewrite(a + 3 * x)
    MA.Test.@test_rewrite(a + x * 4)
    MA.Test.@test_rewrite(a + 3 * x * 4)
    MA.Test.@test_rewrite(x + a)
    MA.Test.@test_rewrite(a - x)
    MA.Test.@test_rewrite(a - 3 * x)
    MA.Test.@test_rewrite(a - x * 4)
    MA.Test.@test_rewrite(a - 3 * x * 4)
    MA.Test.@test_rewrite(x - a)
end

function test_ScalarAffineFunction()
    T = Int
    MA.Test.int_test(
        MOI.ScalarAffineFunction{T},
        exclude = ["int_mul", "int_add_mul"],
    )
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    a = T(2) * x + T(1)
    b = T(4) * y + T(2)
    c = T(3) * x - T(2) * y - T(3)
    d = T(3) * x - T(2) * y + T(4) * x
    e = T(5) * x - T(5)
    f = T(1) * y - T(2) * x + T(2)
    g = T(2) * x + T(3) * x
    _run_all_tests(T, a, b, c, d, e, f, g)
    return
end

function test_ScalarQuadraticFunction()
    T = Int
    MA.Test.int_test(
        MOI.ScalarQuadraticFunction{T},
        exclude = ["int_mul", "int_add_mul"],
    )
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    a = T(2) * x + T(1) + T(4) * x * y
    b = T(4) * y + T(3) * y * y - T(3) * y * x + T(2)
    c = T(2) * x * x + T(3) * x - T(2) * y - T(3)
    d = T(2) * x * x + T(3) * x - T(2) * y + T(4) * x - T(1) * y * y
    e = T(5) * x - T(5) - T(4) * x * y
    f = T(1) * y + T(2) * y * y - T(2) * x + T(2)
    g = T(2) * x + T(3) * x + T(3) * x * y
    _run_all_tests(T, a, b, c, d, e, f, g)
    return
end

end  # module

TestMutableArithmetics.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# This script is not an active part of the package.
# It uses Clang.jl package to parse Sundials C headers and generate
# Julia wrapper for Sundials API.
#
# To run the script from Julia console:
# include(joinpath(Pkg.dir("Sundials"), "src", "wrap_sundials.jl"));

using Clang.wrap_c

# `outpath` specifies, where the julian wrappers would be generated.
# If the generated .jl files are ok, they have to be copied to the "src" folder
# overwriting the old ones
const outpath = normpath(joinpath(dirname(@__FILE__), "..", "new"))
mkpath(outpath)

# Find all relevant Sundials headers
const incpath = normpath(joinpath(dirname(@__FILE__), "..", "deps", "usr", "include"))
if !isdir(incpath)
    error("Sundials C headers not found. Run Pkg.build(\"Sundials\") before trying to wrap C headers.")
end

info("Scanning Sundials headers in $incpath...")
const sundials_folders = filter!(isdir, map!(folder -> joinpath(incpath, folder),
                                 ["nvector", "sundials", "cvode", "cvodes",
                                  "ida", "idas", "kinsol", "arkode",
                                  "sunlinsol","sunmatrix"]))
const sundials_headers = similar(sundials_folders, 0)
for folder in sundials_folders
    info("Processing $folder...")
    append!(sundials_headers,
            map(x->joinpath(folder, x),
                sort!(convert(typeof(sundials_headers), readdir(folder)))))
end
# @show sundials_headers

const clang_path = "/usr/lib/clang/3.8.0" # change to your clang location
const clang_includes = [
    joinpath(clang_path, "include"),
]

# check_use_header(path) = true
# Callback to test if a header should actually be wrapped (for exclusion)
function wrap_header(top_hdr::AbstractString, cursor_header::AbstractString)
    !ismatch(r"(_parallel|_impl)\.h$", cursor_header) &&  # don't wrap parallel and implementation definitions
    (top_hdr == cursor_header) # don't wrap if header is included from the other header (e.g. nvector in cvode or cvode_direct from cvode_band)
end

function wrap_cursor(name::AbstractString, cursor)
    if typeof(cursor) == Clang.cindex.FunctionDecl
        # only wrap API functions
        return ismatch(r"^(CV|KIN|IDA|ARK|N_V|SUN)", name)
    else
        # skip problematic definitions
        return !ismatch(r"^(ABS|SQRT|EXP)$", name)
    end
end

function julia_file(header::AbstractString)
    src_name = basename(dirname(header))
    if src_name == "sundials"
        src_name = "libsundials" # avoid having both Sundials.jl and sundials.jl
    end
    return joinpath(outpath, string(src_name, ".jl"))
end
function library_file(header::AbstractString)
    header_name = basename(header)
    if startswith(header_name, "nvector")
        return "libsundials_nvecserial"
    else
        return string("libsundials_", basename(dirname(header)))
    end
end

const context = wrap_c.init(
    common_file = joinpath(outpath, "types_and_consts.jl"),
    clang_args = [
        "-D", "__STDC_LIMIT_MACROS",
        "-D", "__STDC_CONSTANT_MACROS",
        "-v"
    ],
    clang_diagnostics = true,
    clang_includes = [clang_includes; incpath],
    header_outputfile = julia_file,
    header_library = library_file,
    header_wrapped=wrap_header,
    cursor_wrapped=wrap_cursor
)
context.headers = sundials_headers

# 1st arg name to wrapped arg type map
const arg1_name2type = Dict(
    :arkode_mem => :(ARKODEMemPtr),
    :cvode_mem => :(CVODEMemPtr),
    :cv_mem => :(CVODEMemPtr),
    :kinmem => :(KINMemPtr),
    :kinmemm => :(KINMemPtr), # Sundials typo?
    :ida_mem => :(IDAMemPtr),
    :idaa_mem => :(IDAMemPtr), # Sundials typo?
    :idaadj_mem => :(IDAMemPtr), # Sundials typo?
)

const linear_solvers_and_matrices = [
    "dense",
    "band",
    "pcg",
    "spbcgs",
    "spfgmr",
    "spgmr",
    "sptfqmr",
    "sparse"
    ]

# substitute Ptr{Void} with the typed pointer
const ctor_return_type = Dict(
    "ARKodeCreate" => :(ARKODEMemPtr),
    "CVodeCreate" => :(CVODEMemPtr),
    "IDACreate" => :(IDAMemPtr),
    "KINCreate" => :(KINMemPtr)
)

# signatures for C function pointer types
const FnTypeSignatures = Dict(
    :ARKRhsFn => (:Cint, :((realtype, N_Vector, N_Vector, Ptr{Void}))),
    :CVRhsFn => (:Cint, :((realtype, N_Vector, N_Vector, Ptr{Void}))),
    :CVRootFn => (:Cint, :((realtype, N_Vector, Ptr{realtype}, Ptr{Void}))),
    :IDAResFn => (:Cint, :((realtype, N_Vector, N_Vector, N_Vector, Ptr{Void}))),
    :IDARootFn => (:Cint, :((realtype, N_Vector, N_Vector, Ptr{realtype}, Ptr{Void}))),
    :KINSysFn => (:Cint, :((N_Vector, N_Vector, Ptr{Void}))),
)

wrap_sundials_api(notexpr) = Any[notexpr]

function wrap_sundials_api(expr::Expr)
    if expr.head == :function &&
        expr.args[1].head == :call
        func_name = string(expr.args[1].args[1])
        convert_required = false
        if ismatch(r"^(ARK|CV|KIN|IDA|SUN|N_V)", func_name)
	        @show func_name
            if ismatch(r"Create$", func_name)
                # create functions return typed pointers
		        @assert expr.args[2].args[1].args[3] == :(Ptr{Void})
                expr.args[2].args[1].args[3] = ctor_return_type[func_name]
            elseif length(expr.args[1].args) > 1
		        arg1_type = expr.args[2].args[1].args[4].args[1]
                if arg1_type == :(Ptr{Void}) || arg1_type == :(Ptr{Ptr{Void}})
		            arg1_name = expr.args[1].args[2]
                    arg1_newtype = arg1_name2type[arg1_name]
                    if arg1_type == :(Ptr{Ptr{Void}})
                        # adjust for void** type (for CVodeFree()-like funcs)
                        arg1_newtype = Expr(:curly, :Ref, arg1_newtype)
                    end
		             expr.args[2].args[1].args[4].args[1] = arg1_newtype
                    convert_required = true
                end
            end
            if ismatch(r"UserDataB?$", func_name)
                # replace Ptr{Void} with Any to allow passing Julia objects through user data
                for (i, arg_expr) in enumerate(expr.args[2].args[1].args)
         	        if !(typeof(arg_expr)<:Symbol) && arg_expr.args[1] in values(ctor_return_type)
			            if arg_expr.args[2] == :(Ptr{Void})
		                    arg_expr.args[2] = Any
		                elseif arg_expr.args[3] == :(Ptr{Void})
		    	            arg_expr.args[3] = Any
		                end
			        end
                end
            end
            if !(typeof(expr)<:Symbol) && length(expr.args) > 1 && (expr.args[2].args[1].args[2].args[2] == :libsundials_sunlinsol || expr.args[2].args[1].args[2].args[2] == :libsundials_sunmatrix)
                if func_name[1:6] == "SUNMAT"
                    expr.args[2].args[1].args[2].args[2] =
                    Symbol(string(expr.args[2].args[1].args[2].args[2])*
                    lowercase(split(func_name,"_")[end]))
                else
                    name_i = findfirst(lsmn -> contains(lowercase(func_name),lsmn),linear_solvers_and_matrices)
                    @assert name_i > 0
                    name = linear_solvers_and_matrices[name_i]
                    expr.args[2].args[1].args[2].args[2] = Symbol(string(expr.args[2].args[1].args[2].args[2])*name)
                end
            end
            # generate a higher-level wrapper that converts 1st arg to XXXMemPtr
            # and all other args from Julia objects to low-level C wrappers (e.g. NVector to N_Vector)
            # create wrapers for all arguments that need wrapping
            args_wrap_exprs = map(Base.Iterators.drop(expr.args[1].args, 1)) do expr
                # process each argument
                # return a tuple:
                #   1) high-level wrapper arg name expr
                #   2) expr for local var definition, nothing if not required
                #   3) expr for low-level wrapper call
                # if 1)==3), then no wrapping is required
                if typeof(expr) <: Symbol
			arg_name_expr = expr
			arg_type_expr = Any
		else
			arg_name_expr = expr.args[1]
                	arg_type_expr = expr.args[2]
		end
                if arg_type_expr == :N_Vector
                    # first convert argument to NVector() and store in local var,
                    # this guarantees that the wrapper and associated Sundials object (e.g. N_Vector)
                    # is not removed by GC
                    return (arg_name_expr,
                            Expr(:call, :convert, :NVector, arg_name_expr), # convert arg to NVector to store in a local var
                            Expr(:call, :convert, arg_type_expr, Symbol(string("__", arg_name_expr)))) # convert NVector to N_Vector
                elseif arg_type_expr == :Clong || arg_type_expr == :Cint ||
                    ismatch(r"MemPtr$", string(arg_type_expr))
                    # convert(XXXMemPtr, mem), no local var required
                    return (arg_name_expr, nothing,
                            Expr(:call, :convert, arg_type_expr, arg_name_expr))
                elseif isa(arg_type_expr, Expr) && arg_type_expr.head == :curly &&
                    arg_type_expr.args[1] == :Ptr && arg_type_expr.args[2] != :FILE
                    # convert julia arrays to pointer, no local var required
                    # FIXME sometimes these arguments are not really arrays, but just a pointer to a var to be assigned
                    # by the function call. Does that make sense to detect such cases and assume that input arg is a reference to Julia var?
                    return (arg_name_expr, nothing,
                            Expr(:call, :pointer, arg_name_expr))
                elseif haskey(FnTypeSignatures, arg_type_expr) # wrap Julia function to C function using a defined signature
                    return (arg_name_expr, nothing,
                            Expr(:call, Symbol(string(arg_type_expr, "_wrapper")), arg_name_expr))
                else # any other case, no argument wrapping
                    return (arg_name_expr, nothing, arg_name_expr)
                end
            end
            # check which arguments are wrapped
            if any(map(arg_exprs -> arg_exprs[1] != arg_exprs[3], args_wrap_exprs))
                # mangle the name of the low-level wrapper to avoid recursion
                lowlevel_func_name = string("__", func_name)
                expr.args[1].args[1] = Symbol(lowlevel_func_name)

                # higher-level wrapper function
                wrapper_func_expr = Expr(:function,
                    # function declaration with argument types stripped, so it would accept any type
                    Expr(:call, Symbol(func_name), map(arg_exprs -> arg_exprs[1], args_wrap_exprs)...),
                    Expr(:block,
                    # local var defs
                    map(filter(arg_exprs -> arg_exprs[2] !== nothing, args_wrap_exprs)) do arg_exprs
                        Expr(:(=), Symbol("__", arg_exprs[1]), arg_exprs[2])
                    end...,
                    # low-level function call with Julia types converted to low-level arguments
                    Expr(:call, Symbol(lowlevel_func_name),
                         map(arg_exprs -> arg_exprs[3], args_wrap_exprs)...)
                    )
                )
                # write down both low-level and higher level wrappers
                return Any[expr, wrapper_func_expr]
            else
                return Any[expr]
            end
        end
    elseif expr.head == :typealias && haskey(FnTypeSignatures, expr.args[1])
        # C functional type, generate the wrapper from Julia to C function
        # dummy wrapper that accepts pointer
        fn_typename = expr.args[1]
        fn_rettype, fn_argtypes = FnTypeSignatures[fn_typename]
        wrapper_name = Symbol(string(string(fn_typename), "_wrapper"))
        c_wrapper_def = Expr(:(=),
                Expr(:call, wrapper_name, Expr(:(::), :fp, fn_typename)), :fp)
        jl_wrapper_def = Expr(:(=),
            Expr(:call, wrapper_name, :f),
            # function declaration with argument types stripped
            Expr(:call, :cfunction, :f, fn_rettype, fn_argtypes))
        return Any[expr, c_wrapper_def, jl_wrapper_def]
    elseif expr.head == :const && expr.args[1].head == :(=) &&
           isa(expr.args[1].args[2], Int)
        # fix integer constants should be Cint
        expr.args[1].args[2] = Expr(:call, :Cint, expr.args[1].args[2])
    end
    return Any[expr]
end

context.rewriter = function(exprs)
    mod_exprs = sizehint!(Vector{Any}(), length(exprs))
    for expr in exprs
        append!(mod_exprs, wrap_sundials_api(expr))
    end
    mod_exprs
end

info("Generating .jl wrappers for Sundials in $outpath...")
run(context)
info("Done generating .jl wrappers for Sundials in $outpath")


text:
export AbstractInvariantDistribution

abstract type AbstractInvariantDistribution end

####################
# Pretty printing
####################
function summarise(iv::AbstractInvariantDistribution)
    invdist_type = typeof(iv)
    n_bins = length(iv.dist)
    n_nonzero = length(iv.nonzero_inds)
    invdist_str = "$invdist_type where $n_nonzero out of $n_bins bins have positive measure\n"
    return invdist_str
end

Base.show(io::IO, iv::AbstractInvariantDistribution) = println(io, summarise(iv))


text:
"""
    RangeMatrix(rs::AbstractVector{T}) where T <: AbstractRange

A RangeMatrix is a simple matrix representation of a vector of ranges, with
each range representing one column. Construct a RangeMatrix with a vector of
ranges; the ranges must all have the same length.

Note that it is only efficient when all the ranges are of the same type and in
a concretely typed Vector.
"""
struct RangeMatrix{T,A} <: AbstractArray{T,2}
    rs::A # A <: AbstractVector{_<:Range{T}}
    dims::Tuple{Int,Int}
end
RangeMatrix(rs::AbstractRange...) = RangeMatrix(collect(rs)) # TODO: use tuple storage?
function RangeMatrix(rs::AbstractVector{T}) where T <: AbstractRange
    n = length(rs)
    n == 0 && return RangeMatrix{T}(rs, (0, 0))
    m = length(rs[1])
    for j=2:n
        m == length(rs[j]) || throw(ArgumentError("all ranges must have the same length; expected $m, got $(length(rs[j]))"))
    end
    RangeMatrix{eltype(T), typeof(rs)}(rs, (m, n))
end

Base.size(R::RangeMatrix) = R.dims
Base.IndexStyle(::Type{<:RangeMatrix}) = IndexCartesian()

# Scalar indexing
@inline function Base.getindex(R::RangeMatrix, i::Int, j::Int)
    @boundscheck checkbounds(R, i, j);
    @inbounds return R.rs[j][i]
end

# For non-scalar indexing, only specialize with inner Ranges and Colons to
# return Ranges or RangeMatrixes. For everything else, we can use the fallbacks.
@inline function Base.getindex(R::RangeMatrix, I::Union{AbstractRange, Colon}, J)
    @boundscheck checkbounds(R, I, J)
    unsafe_getindex(R, I, J)
end
@inline unsafe_getindex(R::RangeMatrix, I::Union{AbstractRange, Colon}, j::Real) =
    @inbounds return R.rs[j][I]
@inline unsafe_getindex(R::RangeMatrix, I::Union{AbstractRange, Colon}, ::Colon) =
    @inbounds return RangeMatrix([R.rs[j][I] for j=1:length(R.rs)])
@inline unsafe_getindex(R::RangeMatrix, I::Union{AbstractRange, Colon}, J) =
    @inbounds return RangeMatrix([R.rs[j][I] for j in J])

# We can also optimize bounds checks to only look at each range's endpoints
function Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, R::RangeMatrix)
    b = true
    @inbounds for r in R.rs
        b &= checkindex(Bool, inds, r[1])
        b &= checkindex(Bool, inds, r[end])
    end
    b
end


text:
gauss_nodes_weights(p::Type{P}, n) where {P<:Chebyshev} =
    FastGaussQuadrature.gausschebyshev(n)

gauss_nodes_weights(p::Type{P}, n) where {α,β,P<:Jacobi{α,β}} =
    FastGaussQuadrature.gaussjacobi(n, α, β)

include("Orthogonal/FastLegendre.jl")

eval_basis(::Type{P}, n, x::Float64) where {P<:Legendre} = FastLegendre.fastlegendre(n, x)
eval_basis(::Type{P}, n, x::Union{Int,Int32,Float16,Float32}) where {P<:Legendre} =
    FastLegendre.fastlegendre(n, float(x))
gauss_nodes_weights(p::Type{P}, n) where {P<:Legendre} =
    FastGaussQuadrature.gausslegendre(n)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

broadcast_trafo(f::Any, v_src::Any) = f.(v_src)

function broadcast_trafo(::typeof(identity), v_src::Any)
    v_src
end

function broadcast_trafo(
    trafo::Any,
    v_src::Union{ArrayOfSimilarVectors{<:Real},ShapedAsNTArray}
)
    vs_trg = trafo(elshape(v_src))
    R = eltype(unshaped(trafo(first(v_src)), vs_trg))
    v_src_us = unshaped.(v_src)
    trafo_us = unshaped(trafo)

    n = length(eachindex(v_src_us))
    v_trg_unshaped = nestedview(similar(flatview(v_src_us), R, totalndof(vs_trg), n))
    @assert axes(v_trg_unshaped) == axes(v_src)
    @assert v_trg_unshaped isa ArrayOfSimilarArrays
    @threads for i in eachindex(v_trg_unshaped, v_src)
        v = trafo_us(v_src_us[i])
        v_trg_unshaped[i] = v
    end
    vs_trg.(v_trg_unshaped)
end

function broadcast_trafo(
    ::typeof(identity),
    s_src::DensitySampleVector
)
    s_src
end

function broadcast_trafo(
    trafo::Any,
    s_src::DensitySampleVector
)
    vs_trg = trafo(elshape(s_src.v))
    R = eltype(unshaped(trafo(first(s_src.v)), vs_trg))
    s_src_us = unshaped.(s_src)
    trafo_us = unshaped(trafo)

    n = length(eachindex(s_src_us))
    s_trg_unshaped = DensitySampleVector((
        nestedview(similar(flatview(s_src_us.v), R, totalndof(vs_trg), n)),
        zero(s_src_us.logd),
        deepcopy(s_src_us.weight),
        deepcopy(s_src_us.info),
        deepcopy(s_src_us.aux),
    ))
    @assert axes(s_trg_unshaped) == axes(s_src)
    @assert s_trg_unshaped.v isa ArrayOfSimilarArrays
    @threads for i in eachindex(s_trg_unshaped, s_src)
        v, ladj = with_logabsdet_jacobian(trafo_us, s_src_us.v[i])
        s_trg_unshaped.v[i] = v
        s_trg_unshaped.logd[i] = s_src_us.logd[i] - ladj
    end
    vs_trg.(s_trg_unshaped)
end


text:
using Images, TestImages
using Deconvolution
using FFTW
using ZernikePolynomials #https://github.com/rdoelman/ZernikePolynomials.jl

img = channelview(testimage("cameraman"))

blurring = evaluateZernike(LinRange(-16,16,512), [12, 4, 0], [1.0, -1.0, 2.0], index=:OSA)
blurring = fftshift(blurring)
blurring = blurring ./ sum(blurring)

blurred_img = fft(img) .* fft(blurring) |> ifft |> real

@time result_200_iterations = lucy(blurred_img, blurring, iterations=200)
@time result_2000_iterations = lucy(blurred_img, blurring, iterations=2000)

(N, M) = size(blurred_img)

save("original.jpg", Images.clamp01nan.(img))
save("blurred.jpg", Images.clamp01nan.(blurred_img))
save("restored_200.jpg", Images.clamp01nan.(result_200_iterations))
save("restored_2000.jpg", Images.clamp01nan.(result_2000_iterations))


text:
function centering!(mechanism::Mechanism, αaff::T) where T
    system = mechanism.system

    n = 0
    ν = 0.0
    νaff = 0.0
    for contact in mechanism.contacts
        ν, νaff, n = centering!(ν, νaff, n, mechanism, contact, get_entry(system, contact.id), αaff)
    end

    for joint in mechanism.joints
        ν, νaff, n = centering!(ν, νaff, n, mechanism, joint, get_entry(system, joint.id), αaff)
    end

    ν /= n
    νaff /= n
    return ν, νaff
end

function centering!(ν, νaff, n, mechanism, contact::ContactConstraint{T,N,Nc,Cs,N½}, vector_entry::Entry, αaff) where {T,N,Nc,Cs,N½}
    s = contact.impulses_dual[2]
    γ = contact.impulses[2]
    Δs = vector_entry.value[1:N½]
    Δγ = vector_entry.value[N½ .+ (1:N½)]
    ν += dot(s, γ)
    νaff += dot(s + αaff * Δs, γ + αaff * Δγ) # plus or minus
    n += cone_degree(contact)
    return ν, νaff, n
end

function centering!(ν, νaff, n, mechanism, joint::JointConstraint{T,N,Nc}, vector_entry::Entry, αaff) where {T,N,Nc}
    for (i, element) in enumerate([joint.translational, joint.rotational])
        s, γ = split_impulses(element, joint.impulses[2][joint_impulse_index(joint,i)])
        Δs, Δγ = split_impulses(element, vector_entry.value[joint_impulse_index(joint,i)])
        ν += dot(s, γ)
        νaff += dot(s + αaff * Δs, γ + αaff * Δγ) # plus or minus
        n += length(s)
    end
    return ν, νaff, n
end


text:
"""
$(TYPEDEF)

Epigraph of real or complex matrix spectral norm (i.e. maximum singular value)
for a matrix (stacked column-wise) of `nrows` rows and `ncols` columns with
`nrows ≤ ncols`.

    $(FUNCTIONNAME){T, R}(nrows::Int, ncols::Int, use_dual::Bool = false)
"""
mutable struct EpiNormSpectral{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d1::Int
    d2::Int
    is_complex::Bool

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_aux_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}

    W_svd
    s::Vector{T}
    U::Matrix{R}
    Vt::Matrix{R}
    z::Vector{T}
    uzi::Vector{T}
    Urzi::Matrix{R}
    Vtsrzi::Matrix{R}
    cu::T
    zti1::T
    tzi::Vector{T}
    usti::Vector{T}
    zszidd::Matrix{T}
    zstidd::Matrix{T}

    w1::Matrix{R}
    w2::Matrix{R}
    w3::Matrix{R}
    s1::Vector{T}
    s2::Vector{T}
    U1::Matrix{R}
    U2::Matrix{R}
    U3::Matrix{R}

    function EpiNormSpectral{T, R}(
        d1::Int,
        d2::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert 1 <= d1 <= d2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.is_complex = (R <: Complex)
        cone.dim = 1 + vec_length(R, d1 * d2)
        cone.d1 = d1
        cone.d2 = d2
        return cone
    end
end

reset_data(cone::EpiNormSpectral) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_aux_updated =
    cone.hess_fact_updated = false)

function setup_extra_data!(
    cone::EpiNormSpectral{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    (d1, d2) = (cone.d1, cone.d2)
    cone.z = zeros(T, d1)
    cone.uzi = zeros(T, d1)
    cone.Urzi = zeros(R, d1, d1)
    cone.Vtsrzi = zeros(R, d1, d2)
    cone.tzi = zeros(T, d1)
    cone.usti = zeros(T, d1)
    cone.zszidd = zeros(T, d1, d1)
    cone.zstidd = zeros(T, d1, d1)
    cone.w1 = zeros(R, d1, d2)
    cone.w2 = zeros(R, d1, d2)
    cone.w3 = zeros(R, d1, d2)
    cone.s1 = zeros(T, d1)
    cone.s2 = zeros(T, d1)
    cone.U1 = zeros(R, d1, d1)
    cone.U2 = zeros(R, d1, d1)
    cone.U3 = zeros(R, d1, d1)
    return cone
end

get_nu(cone::EpiNormSpectral) = 1 + cone.d1

function set_initial_point!(
    arr::AbstractVector{T},
    cone::EpiNormSpectral{T},
    ) where {T <: Real}
    arr .= 0
    arr[1] = sqrt(T(get_nu(cone)))
    return arr
end

function update_feas(cone::EpiNormSpectral{T}) where {T <: Real}
    # TODO speed up using norm bounds, cholesky of Z?
    @assert !cone.feas_updated
    u = cone.point[1]

    if u > eps(T)
        @views vec_copyto!(cone.w1, cone.point[2:end])
        cone.W_svd = svd(cone.w1, full = false) # TODO in place
        cone.is_feas = (u - maximum(cone.W_svd.S) > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::EpiNormSpectral{T}) where {T <: Real}
    # TODO speed up using norm bound, sum sqrt eigvals of WW'?
    u = cone.dual_point[1]
    if u > eps(T)
        W = @views vec_copyto!(cone.w1, cone.dual_point[2:end])
        return (u - sum(svdvals!(W)) > eps(T))
    end
    return false
end

function update_grad(cone::EpiNormSpectral{T}) where T
    @assert cone.is_feas
    u = cone.point[1]
    U = cone.U = cone.W_svd.U
    Vt = cone.Vt = cone.W_svd.Vt
    s = cone.s = cone.W_svd.S
    z = cone.z
    uzi = cone.uzi
    s1 = cone.s1
    rz = cone.s2
    w1 = cone.w1
    g = cone.grad

    cone.cu = (cone.d1 - 1) / u
    @. z = (u - s) * (u + s)
    @. uzi = 2 * u / z
    @. rz = sqrt(z)
    @. s1 = inv(rz)
    mul!(cone.Urzi, U, Diagonal(s1))
    @. s1 = s / rz
    mul!(cone.Vtsrzi, Diagonal(s1), Vt)

    g[1] = cone.cu - sum(uzi)
    mul!(w1, cone.Urzi, cone.Vtsrzi, 2, false)
    @views vec_copyto!(g[2:end], w1)

    cone.grad_updated = true
    return cone.grad
end

function update_hess_aux(cone::EpiNormSpectral{T}) where T
    @assert !cone.hess_aux_updated
    @assert cone.grad_updated
    d1 = cone.d1
    u = cone.point[1]
    s = cone.s
    z = cone.z
    tzi = cone.tzi
    usti = cone.usti
    zszidd = cone.zszidd
    zstidd = cone.zstidd

    u2 = abs2(u)
    zti1 = one(u)
    @inbounds for j in 1:d1
        s_j = s[j]
        z_j = z[j]
        for i in 1:(j - 1)
            s_i = s[i]
            z_i = z[i]
            s_ij = s_i * s_j
            z_ij = u2 - s_ij
            t_ij = u2 + s_ij
            # zszidd and zstidd are nonsymmetric
            zszidd[i, j] = z_i / z_ij * s_j
            zszidd[j, i] = z_j / z_ij * s_i
            zstidd[i, j] = z_i / t_ij * s_j
            zstidd[j, i] = z_j / t_ij * s_i
        end
        t_j = u2 + abs2(s_j)
        zt_ij = z_j / t_j
        zti1 += zt_ij
        tzi[j] = t_j / z_j
        usti[j] = 2 * u * s_j / t_j
        zszidd[j, j] = s_j
        zstidd[j, j] = zt_ij * s_j
    end
    cone.zti1 = zti1

    cone.hess_aux_updated = true
end

function update_hess(cone::EpiNormSpectral)
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :hess) || alloc_hess!(cone)
    d1 = cone.d1
    d2 = cone.d2
    u = cone.point[1]
    z = cone.z
    uzi = cone.uzi
    Urzi = cone.Urzi
    Vtsrzi = cone.Vtsrzi
    tzi = cone.tzi
    w1 = cone.w1
    w2 = cone.w2
    w3 = cone.w3
    U1 = cone.U1
    U2 = cone.U2
    H = cone.hess.data

    # u, u
    @inbounds H[1, 1] = -cone.cu / u + 2 * sum(tzi[i] / z[i] for i in 1:d1)

    # u, w
    mul!(U1, Urzi, Diagonal(uzi), -2, false)
    mul!(w1, U1, Vtsrzi)
    @views vec_copyto!(H[1, 2:end], w1)

    # w, w
    Urzit = copyto!(w1, Urzi') # accessing columns below
    c_idx = 2
    reim1s = (cone.is_complex ? [1, im] : [1,])
    @inbounds for j in 1:d2, i in 1:d1, reim1 in reim1s
        @views Urzi_i = Urzit[:, i]
        @views mul!(U1, Urzi_i, Vtsrzi[:, j]', reim1, false)
        @. U2 = U1 + U1'
        mul!(w2, Hermitian(U2, :U), Vtsrzi)
        @. @views w2[:, j] += reim1 * Urzi_i
        mul!(w3, Urzi, w2, 2, false)
        @views vec_copyto!(H[2:end, c_idx], w3)
        c_idx += 1
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormSpectral,
    )
    cone.hess_aux_updated || update_hess_aux(cone)
    d1 = cone.d1
    u = cone.point[1]
    z = cone.z
    Urzi = cone.Urzi
    Vtsrzi = cone.Vtsrzi
    cu = cone.cu
    tzi = cone.tzi
    w1 = cone.w1
    w2 = cone.w2
    U1 = cone.U1
    U2 = cone.U2
    Duzi = Diagonal(cone.uzi)

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views vec_copyto!(w1, arr[2:end, j])
        mul!(w2, Urzi', w1)
        mul!(U1, w2, Vtsrzi')
        @. U2 = U1 + U1'

        prod[1, j] = -cu * p / u + 2 * sum((p * tzi[i] -
            u * real(U2[i, i])) / z[i] for i in 1:d1)

        @. U2 -= p * Duzi
        mul!(w2, Hermitian(U2, :U), Vtsrzi, true, true)
        mul!(w1, Urzi, w2, 2, false)
        @views vec_copyto!(prod[2:end, j], w1)
    end

    return prod
end

function update_inv_hess(cone::EpiNormSpectral{T}) where {T <: Real}
    cone.hess_aux_updated || update_hess_aux(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    d1 = cone.d1
    d2 = cone.d2
    u = cone.point[1]
    U = cone.U
    Vt = cone.Vt
    z = cone.z
    zti1 = cone.zti1
    usti = cone.usti
    zszidd = cone.zszidd
    zstidd = cone.zstidd
    s1 = cone.s1
    w1 = cone.w1
    w2 = cone.w2
    w3 = cone.w3
    U1 = cone.U1
    U2 = cone.U2
    Hi = cone.inv_hess.data

    # u, u
    huu = u / zti1 * u
    Hi[1, 1] = huu

    # u, w
    @. s1 = huu * usti
    mul!(U1, U, Diagonal(s1))
    mul!(w1, U1, Vt)
    @views vec_copyto!(Hi[1, 2:end], w1)

    # w, w
    Ut = copyto!(w1, U') # accessing columns below
    c_idx = 2
    reim1s = (cone.is_complex ? [1, im] : [1,])
    @inbounds for j in 1:d2, i in 1:d1, reim1 in reim1s
        @views U_i = Ut[:, i]
        @views mul!(U1, U_i, Vt[:, j]', reim1, false)
        U1 .*= zszidd
        @. U2 = (U1 + U1') * zstidd
        mul!(w2, U2, Vt, -1, false)
        @. @views w2[:, j] += reim1 * z * U_i
        mul!(w3, U, w2, T(0.5), false)
        @views vec_copyto!(Hi[2:end, c_idx], w3)
        c_idx += 1
    end

    rthuu = sqrt(huu)
    @. s1 = rthuu * usti
    mul!(U1, U, Diagonal(s1))
    mul!(w1, U1, Vt)
    @views Hiuw2vec = Hi[2:end, 1]
    vec_copyto!(Hiuw2vec, w1)
    @views mul!(Hi[2:end, 2:end], Hiuw2vec, Hiuw2vec', true, true)

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiNormSpectral{T},
    ) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    d1 = cone.d1
    u = cone.point[1]
    U = cone.U
    Vt = cone.Vt
    z = cone.z
    zti1 = cone.zti1
    usti = cone.usti
    zszidd = cone.zszidd
    zstidd = cone.zstidd
    w1 = cone.w1
    w2 = cone.w2
    U1 = cone.U1
    U2 = cone.U2
    Dusti = Diagonal(usti)

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views vec_copyto!(w1, arr[2:end, j])
        mul!(w2, U', w1)
        mul!(U1, w2, Vt')

        c1 = u * (p + sum(usti[i] * real(U1[i, i]) for i in 1:d1)) / zti1 * u
        prod[1, j] = c1

        U1 .*= zszidd
        @. U2 = (U1 + U1') * zstidd - 2 * c1 * Dusti

        lmul!(Diagonal(z), w2)
        mul!(w2, U2, Vt, -1, true)
        mul!(w1, U, w2, T(0.5), false)
        @views vec_copyto!(prod[2:end, j], w1)
    end

    return prod
end

function dder3(cone::EpiNormSpectral{T}, dir::AbstractVector{T}) where T
    cone.hess_aux_updated || update_hess_aux(cone)
    d1 = cone.d1
    u = cone.point[1]
    z = cone.z
    uzi = cone.uzi
    Urzi = cone.Urzi
    Vtsrzi = cone.Vtsrzi
    w1 = cone.w1
    w2 = cone.w2
    s1 = cone.s1
    s2 = cone.s2
    U1 = cone.U1
    U2 = cone.U2
    U3 = cone.U3
    dder3 = cone.dder3
    Ds1 = Diagonal(s1)

    p = dir[1]
    @views vec_copyto!(w1, dir[2:end])

    @. s1 = p * uzi
    @. s2 = 2 * u * s1 - p

    mul!(w2, Urzi', w1)
    mul!(U3, w2, Vtsrzi')
    @. U1 = U3 + U3'

    mul!(U3, Diagonal(uzi), U1)
    @. U2 = U3 + U3'

    mul!(U3, w2, w2')
    mul!(U3, Hermitian(U1), U1, true, true)
    @inbounds tr1 = sum(real(U3[i, i]) * uzi[i] for i in 1:d1)

    @. U1 -= Ds1
    @inbounds tr2 = sum((s1[i] * (s2[i] + 2 * p) +
        2 * s2[i] * real(U1[i, i])) / z[i] for i in 1:d1)

    dder3[1] = -cone.cu * abs2(p / u) + tr1 - tr2

    @. s1 = s2 / z
    @. U3 += p * (Ds1 - U2)

    mul!(w1, Hermitian(U1, :U), w2)
    mul!(w1, Hermitian(U3, :U), Vtsrzi, true, true)
    mul!(w2, Urzi, w1, -2, false)
    @views vec_copyto!(dder3[2:end], w2)

    return dder3
end


text:
module TestScikitLearn

import MLJModels, ScikitLearn
using MLJBase
using Tables
using Test
using LinearAlgebra
using Random
using MLJModels.ScikitLearn_
using CategoricalArrays

include("../testutils.jl")
include("clustering.jl")
include("svm.jl")
include("linear-regressors.jl")
include("linear-classifiers.jl")
include("gaussian-process.jl")
include("ensemble.jl")
include("discriminant-analysis.jl")
include("misc.jl")

end
true


text:
# Matrix type
@cenum(MatrixType,
    REAL_SYM            = 1,
    REAL_SYM_POSDEF     = 2,
    REAL_SYM_INDEF      = -2,
    COMPLEX_STRUCT_SYM  = 3,
    COMPLEX_HERM_POSDEF = 4,
    COMPLEX_HERM_INDEF  = -4,
    COMPLEX_SYM         = 6,
    REAL_NONSYM         = 11,
    COMPLEX_NONSYM      = 13,
)

Base.isreal(v::MatrixType) = v in (REAL_SYM, REAL_SYM_POSDEF, REAL_SYM_INDEF, REAL_NONSYM)
issymmetric(v::MatrixType) = v in (REAL_SYM, REAL_SYM_POSDEF, REAL_SYM_INDEF, COMPLEX_STRUCT_SYM,
                                        COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF, COMPLEX_SYM)
Base.ishermitian(v::MatrixType) = v in (REAL_SYM_POSDEF, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF)
isposornegdef(v::MatrixType) = v in (REAL_SYM_POSDEF, REAL_SYM_INDEF, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF)

const MATRIX_STRING = Dict{MatrixType, String}(
    REAL_SYM => "Real structurally symmetric",
    REAL_SYM_POSDEF => "Real symmetric positive definite",
    REAL_SYM_INDEF => "Real symmetric indefinite",
    COMPLEX_STRUCT_SYM => "Complex structurally symmetric",
    COMPLEX_HERM_POSDEF => "Complex Hermitian postive definite",
    COMPLEX_HERM_INDEF => "Complex Hermitian indefinite",
    COMPLEX_SYM  => "Complex symmetric",
    REAL_NONSYM  => "Real nonsymmetric",
    COMPLEX_NONSYM  => "Complex nonsymmetric"
)

const REAL_MATRIX_TYPES = [REAL_SYM, REAL_SYM_POSDEF, REAL_SYM_INDEF, REAL_NONSYM]
const COMPLEX_MATRIX_TYPES = [COMPLEX_STRUCT_SYM, COMPLEX_HERM_POSDEF, COMPLEX_HERM_INDEF, COMPLEX_NONSYM]

# Messages
@cenum(MessageLevel,
    MESSAGE_LEVEL_OFF = 0,
    MESSAGE_LEVEL_ON = 1
)

# Solver
@cenum(Solver,
    DIRECT_SOLVER = 0,
    ITERATIVE_SOLVER = 1
)

const SOLVER_STRING = Dict{Solver, String}(
    DIRECT_SOLVER => "Direct solver",
    ITERATIVE_SOLVER => "Iterative solver"
)

# Phase
@cenum(Phase,
    ANALYSIS                             = 11,
    ANALYSIS_NUM_FACT                    = 12,
    ANALYSIS_NUM_FACT_SOLVE_REFINE       = 13,
    NUM_FACT                             = 22,
    SELECTED_INVERSION                   = -22,
    NUM_FACT_SOLVE_REFINE                = 23,
    SOLVE_ITERATIVE_REFINE               = 33,
    SOLVE_ITERATIVE_REFINE_ONLY_FORWARD  = 331,
    SOLVE_ITERATIVE_REFINE_ONLY_DIAG     = 332,
    SOLVE_ITERATIVE_REFINE_ONLY_BACKWARD = 333,
    RELEASE_LU_MNUM                      = 0,
    RELEASE_ALL                           = -1
)

const PHASE_STRING = Dict{Phase, String}(
    ANALYSIS                             => "Analysis",
    ANALYSIS_NUM_FACT                    => "Analysis, numerical factorization",
    ANALYSIS_NUM_FACT_SOLVE_REFINE       => "Analysis, numerical factorization, solve, iterative refinement",
    NUM_FACT                             => "Numerical factorization",
    SELECTED_INVERSION                   => "Selected Inversion",
    NUM_FACT_SOLVE_REFINE                => "Numerical factorization, solve, iterative refinement",
    SOLVE_ITERATIVE_REFINE               => "Solve, iterative refinement",
    RELEASE_LU_MNUM                      => "Release internal memory for L and U matrix number MNUM",
    RELEASE_ALL                           => "Release all internal memory for all matrices",
    SOLVE_ITERATIVE_REFINE_ONLY_FORWARD  => "like phase=33, but only forward substitution",
    SOLVE_ITERATIVE_REFINE_ONLY_DIAG     => "like phase=33, but only diagonal substitution (if available)",
    SOLVE_ITERATIVE_REFINE_ONLY_BACKWARD => "like phase=33, but only backward substitution",
)