text:
"""
    Euler(; h = 0.0) :: ExplicitRungeKuttaSolver
    ExplicitEuler(args...; kwargs...) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 1st-order Euler method.
```
"""
function Euler(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        0 0;
        1 1;
    ])
    return ERK(tableau, h)
end
@doc (@doc Euler) ExplicitEuler(args...; kwargs...) = Euler(args...; kwargs...)

"""
    Heun2(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 2nd-order Heun method.
"""
function Heun2(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        0   0   0 ;
        1   1   0 ;
        2  1/2 1/2;
    ])
    return ERK(tableau, h)
end

"""
    Ralston2(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 2nd-order Ralston method.
"""
function Ralston2(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0 ;
        2/3 2/3  0 ;
         2  1/4 3/4;
    ])
    return ERK(tableau, h)
end

"""
    Midpoint(; h = 0.0) :: ExplicitRungeKuttaSolver
    ExplicitMidpoint(args...; kwargs...) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 2nd-order midpoint method.
"""
function Midpoint(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0;
        1/2 1/2  0;
         2   0   1;
    ])
    return ERK(tableau, h)
end
@doc (@doc Midpoint) ExplicitMidpoint(args...; kwargs...) = Midpoint(args...; kwargs...)

"""
    Heun3(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 3rd-order Heun method.
"""
function Heun3(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0   0 ;
        1/3 1/3  0   0 ;
        2/3  0  2/3  0 ;
         3  1/4  0  3/4;
    ])
    return ERK(tableau, h)
end

"""
    Kutta3(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 3rd-order Kutta method.
"""
function Kutta3(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0   0 ;
        1/2 1/2  0   0 ;
         1  -1   2   0 ;
         3  1/6 2/3 1/6;
    ])
    return ERK(tableau, h)
end

"""
    Ralston3(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 3rd-order Ralston method.
"""
function Ralston3(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0   0 ;
        1/2 1/2  0   0 ;
        3/4  0  3/4  0 ;
         3  2/9 1/3 4/9;
    ])
    return ERK(tableau, h)
end

"""
    SSPRK3(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 3rd-order Strong-Stability-Preserving Runge-Kutta method.
"""
function SSPRK3(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0   0 ;
         1   1   0   0 ;
        1/2 1/4 1/4  0 ;
         3  1/6 1/6 2/3;
    ])
    return ERK(tableau, h)
end

"""
    RK4(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 4th-order Runge-Kutta method.
"""
function RK4(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0   0   0   0 ;
        1/2 1/2  0   0   0 ;
        1/2  0  1/2  0   0 ;
         1   0   0   1   0 ;
         4  1/6 1/3 1/3 1/6;
    ])
    return ERK(tableau, h)
end

"""
    Rule38(; h = 0.0) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 4th-order 3/8-rule method.
"""
function Rule38(; h = 0.0)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0    0    0    0    0 ;
        1/3  1/3   0    0    0 ;
        2/3 -1/3   1    0    0 ;
         1    1   -1    1    0 ;
         4   1/8  3/8  3/8  1/8;
    ])
    return ERK(tableau, h)
end

"""
    HeunEuler(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 2nd-order Heun-Euler method with 1st-order error estimate.
"""
function HeunEuler(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        0   0   0 ;
        1   1   0 ;
        2  1/2 1/2;
        1   1   0 ;
    ])
    adaptive = AdaptiveParameters(δ=δ, ϵ=ϵ, K=K)
    return ERK(tableau, h, adaptive)
end

"""
    Fehlberg45(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100) :: ExplicitRungeKuttaSolver
    F45(args...; kwargs...) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 4th-order Fehlberg method with 5th-order error estimate.
"""
function Fehlberg45(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
          0       0          0          0           0         0    0  ;
         1/4     1/4         0          0           0         0    0  ;
         3/8     3/32       9/32        0           0         0    0  ;
        12/13 1932/2197 -7200/2197  7296/2197       0         0    0  ;
          1    439/216      -8      3680/513    -845/4104     0    0  ;
         1/2    -8/27        2     -3544/2565   1859/4104  -11/40  0  ;
          4     25/216       0      1408/2565   2197/4104   -1/5   0  ;
          5     16/135       0      6656/12825 28561/56430  -9/50 2/55;
    ])
    adaptive = AdaptiveParameters(δ=δ, ϵ=ϵ, K=K)
    return ERK(tableau, h, adaptive)
end
@doc (@doc Fehlberg45) F45(args...; kwargs...) = Fehlberg45(args...; kwargs...)

"""
    DormandPrince54(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100) :: ExplicitRungeKuttaSolver
    DP54(args...; kwargs...) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 5th-order Dormand-Prince method with 4th-order error estimate.
"""
function DormandPrince54(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0        0            0          0          0          0          0      0  ;
        1/5      1/5           0          0          0          0          0      0  ;
        3/10     3/40         9/40        0          0          0          0      0  ;
        4/5     44/45       -56/15      32/9         0          0          0      0  ;
        8/9  19372/6561  -25360/2187 64448/6561  -212/729       0          0      0  ;
         1    9017/3168    -355/33   46732/5247    49/176  -5103/18656     0      0  ;
         1      35/384         0       500/1113   125/192  -2187/6784    11/84    0  ;
         5      35/384         0       500/1113   125/192  -2187/6784    11/84    0  ;
         4    5179/57600       0      7571/16695  393/640 -92097/339200 187/2100 1/40;
    ])
    adaptive = AdaptiveParameters(δ=δ, ϵ=ϵ, K=K)
    return ERK(tableau, h, adaptive)
end
@doc (@doc DormandPrince54) DP54(args...; kwargs...) = DormandPrince54(args...; kwargs...)

"""
    Verner65(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100) :: ExplicitRungeKuttaSolver
    V65(args...; kwargs...) :: ExplicitRungeKuttaSolver

returns an [`ExplicitRungeKuttaSolver`](@ref) for the 6th-order Verner method with 5th-order error estimate.
"""
function Verner65(; h = 0.0, δ = 0.0, ϵ = 1e-5, K = 100)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0        0          0         0          0          0       0       0        0   ;
        1/6      1/6         0         0          0          0       0       0        0   ;
        4/15     4/75      16/75       0          0          0       0       0        0   ;
        2/3      5/6       -8/3       5/2         0          0       0       0        0   ;
        5/6   -165/64      55/6    -425/64      85/96        0       0       0        0   ;
         1      12/5        -8     4015/612    -11/36      88/255    0       0        0   ;
        1/15 -8263/15000  124/75   -643/680    -81/250   2484/10625  0       0        0   ;
         1    3501/1720  -300/43 297275/52632 -319/2322 24068/84065  0   3850/26703   0   ;
         6       3/40        0      875/2244    23/72     264/1955   0    125/11592 43/616;
         5      13/160       0     2375/5984     5/16      12/85    3/44     0        0   ;
    ])
    adaptive = AdaptiveParameters(δ=δ, ϵ=ϵ, K=K)
    return ERK(tableau, h, adaptive)
end
@doc (@doc Verner65) V65(args...; kwargs...) = Verner65(args...; kwargs...)

"""
    BackwardEuler(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver
    ImplicitEuler(args...; kwargs...) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 1st-order backward Euler method.
"""
function BackwardEuler(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        1 1;
        1 1;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end
@doc (@doc BackwardEuler) ImplicitEuler(args...; kwargs...) = BackwardEuler(args...; kwargs...)

"""
    ImplicitMidpoint(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 2nd-order implicit midpoint method.
"""
function ImplicitMidpoint(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        1/2 1/2;
         2   1 ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    CrankNicolson(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 2nd-order Crank-Nicolson method.
"""
function CrankNicolson(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        0   0   0 ;
        1  1/2 1/2;
        2  1/2 1/2;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    SDIRK3(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 3rd-order SDIRK method.
"""
function SDIRK3(; h = 0.0, ϵ = 1e-3, K = 10)
    γ = 1/2 + √3/6
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         γ    γ   0 ;
        1-γ 1-2γ  γ ;
         3   1/2 1/2;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    GaussLegendre4(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver
    GL4(args...; kwargs...) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 4th-order Gauss-Legendre method.
"""
function GaussLegendre4(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        1/2-√3/6   1/4    1/4-√3/6;
        1/2+√3/6 1/4+√3/6   1/4   ;
           4       1/2      1/2   ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end
@doc (@doc GaussLegendre4) GL4(args...; kwargs...) = GaussLegendre4(args...; kwargs...)

"""
    GaussLegendre6(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver
    GL6(args...; kwargs...) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 6th-order Gauss–Legendre method.
"""
function GaussLegendre6(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        1/2-√15/10     5/36    2/9-√15/15 5/36-√15/30;
           1/2     5/36+√15/24    2/9     5/36-√15/24;
        1/2+√15/10 5/36+√15/30 2/9+√15/15     5/36   ;
            6          5/18       4/9         5/18   ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end
@doc (@doc GaussLegendre6) GL6(args...; kwargs...) = GaussLegendre6(args...; kwargs...)

"""
    LobattoIIIA4(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 4th-order Lobatto IIIA method.
"""
function LobattoIIIA4(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   0    0    0  ;
        1/2 5/24 1/3 -1/24;
         1  1/6  2/3  1/6 ;
         4  1/6  2/3  1/6 ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    LobattoIIIB2(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 2nd-order Lobatto IIIB method.
"""
function LobattoIIIB2(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        1/2 1/2  0 ;
        1/2 1/2  0 ;
         2  1/2 1/2;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    LobattoIIIB4(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 4th-order Lobatto IIIB method.
"""
function LobattoIIIB4(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0  1/6 -1/6  0 ;
        1/2 1/6  1/3  0 ;
         1  1/6  5/6  0 ;
         4  1/6  2/3 1/6;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    LobattoIIIC2(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 2nd-order Lobatto IIIC method.
"""
function LobattoIIIC2(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        0  1/2 -1/2;
        1  1/2  1/2;
        2  1/2  1/2;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    LobattoIIIC4(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 4th-order Lobatto IIIC method.
"""
function LobattoIIIC4(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0   1/6 -1/3   1/6 ;
        1/2  1/6  5/12 -1/12;
         1   1/6  2/3   1/6 ;
         4   1/6  2/3   1/6 ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    RadauIA3(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 3rd-order Radau IA method.
"""
function RadauIA3(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
         0  1/4 -1/4 ;
        2/3 1/4  5/12;
         3  1/4  3/4 ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    RadauIA5(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 5th-order Radau IA method.
"""
function RadauIA5(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
            0     1/9   -1/18-√6/18     -1/18+√6/18  ;
        3/5-√6/10 1/9  11/45+7*√6/360 11/45-43*√6/360;
        3/5+√6/10 1/9 11/45+43*√6/360  11/45-7*√6/360;
            5     1/9    4/9+√6/36       4/9-√6/36   ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    RadauIIA3(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 3rd-order Radau IIA method.
"""
function RadauIIA3(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        1/3  5/12 -1/12;
         1   3/4   1/4 ;
         3   3/4   1/4 ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    RadauIIA5(; h = 0.0, ϵ = 1e-3, K = 10) :: ImplicitRungeKuttaSolver

returns an [`ImplicitRungeKuttaSolver`](@ref) for the 5th-order Radau IIA method.
"""
function RadauIIA5(; h = 0.0, ϵ = 1e-3, K = 10)
    h_T = typeof(h)
    tableau = ButcherTableau(h_T[
        2/5-√6/10   11/45-7*√6/360   37/225-169*√6/1800 -2/225+√6/75;
        2/5+√6/10 37/225+169*√6/1800   11/45+7*√6/360   -2/225-√6/75;
           1           4/9-√6/36          4/9+√6/36          1/9    ;
           5           4/9-√6/36          4/9+√6/36          1/9    ;
    ])
    newton = NewtonParameters(ϵ=ϵ, K=K)
    return IRK(tableau, h, newton)
end

"""
    ExponentialRK4(; h = 0.0) :: ExplicitExponentialRungeKuttaSolver

returns an [`ExplicitExponentialRungeKuttaSolver`](@ref) for the 4th-order EERK method.
"""
function ExponentialRK4(; h = 0.0)
    return EERK(h)
end
@doc (@doc ExponentialRK4) ERK4(args...; kwargs...) = ExponentialRK4(args...; kwargs...)


text: println("Goodbye World")

ext: jl

max_stars_repo_name: KorfLab/julie

max_issues_repo_name: KorfLab/julie

max_forks_repo_name: KorfLab/julie


text:
# Generic implementation: arbitrary list of elements on the intervals between the nodes
struct Piecewise{T,F<:FunctionSpace{T},S<:AbstractVector{<:Real}} <: AbstractPiecewise{T,F}
    nodes::S
    elements::Vector{F}
    function Piecewise(nodes::S, elements::Vector{F}) where {T, F<:FunctionSpace{T}, S<:AbstractVector}
        @assert length(nodes) == length(elements) + 1
        @assert eltype(nodes) <: Real
        if nodes isa AbstractRange
            @assert step(nodes) > 0
        else
            @assert issorted(nodes)
        end
        return new{T, F, S}(nodes, elements)
    end
end

# Basic properties
nodes(p::Piecewise) = p.nodes
elements(p::Piecewise) = p.elements
nodevalues(p::Piecewise) = Base.Generator(p, nodes(p))

Base.length(p::Piecewise) = length(p.elements)

Base.:(==)(p1::Piecewise, p2::Piecewise) =
    nodes(p1) == nodes(p2) && elements(p1) == elements(p2)

# Indexing, getting and setting elements
@inline Base.getindex(p::Piecewise, i) = getindex(p.elements, i)
@inline Base.setindex!(p::Piecewise, v, i) = setindex!(p.elements, v, i)

# Use as function
function (p::Piecewise)(x)
    x₀ = first(nodes(p))
    x₁ = last(nodes(p))
    (x < x₀ || x > x₁) && throw(DomainError())
    x == x₀ && return p[1](x)
    x == x₁ && return p[length(p)](x)
    if nodes isa AbstractRange
        ir = (x-x₀)/step(nodes(p))
        i = floor(Int, ir)
        return p[i](x)
    else
        j = findlast(<=(x), nodes(p))
        @assert j !== nothing && j !== length(nodes(p))
        return p[j](x)
    end
end

# Change number of coefficients
function truncate!(p::Piecewise; kwargs...)
    for i = 1:length(p)
        truncate!(p[i]; kwargs..., dx = (p.nodes[i+1]-p.nodes[i])/2)
    end
    return p
end

# Special purpose constructor
Base.similar(p::Piecewise, args...) =
    Piecewise(nodes(p), map(x->similar(x, args...), elements(p)))
Base.zero(p::Piecewise) = Piecewise(nodes(p), map(zero, elements(p)))
Base.one(p::Piecewise) = Piecewise(nodes(p), map(one, elements(p)))

# Arithmetic (out of place)
Base.copy(p::AbstractPiecewise) = Piecewise(copy(nodes(p)), map(copy, elements(p)))

Base.:+(p::AbstractPiecewise) = Piecewise(nodes(p), map(+, elements(p)))
Base.:-(p::AbstractPiecewise) = Piecewise(nodes(p), map(-, elements(p)))

function Base.:+(p1::AbstractPiecewise, p2::AbstractPiecewise)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), map(+, elements(p1), elements(p2)))
end

function Base.:-(p1::AbstractPiecewise, p2::AbstractPiecewise)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), map(-, elements(p1), elements(p2)))
end

Base.:*(p::AbstractPiecewise, a::Const) =
    Piecewise(nodes(p), map(Base.Fix2(*, a), elements(p)))
Base.:*(a::Const, p::AbstractPiecewise) =
    Piecewise(nodes(p), map(Base.Fix1(*, a), elements(p)))
Base.:/(p::AbstractPiecewise, a) =
    Piecewise(nodes(p), map(Base.Fix2(/, a), elements(p)))
Base.:\(a, p::AbstractPiecewise) =
    Piecewise(nodes(p), map(Base.Fix1(\, a), elements(p)))

Base.:*(p1::AbstractPiecewise, p2::AbstractPiecewise) = truncmul(p1, p2)
function truncmul(p1::AbstractPiecewise, p2::AbstractPiecewise; kwargs...)
    n = nodes(p1)
    @assert n == nodes(p2)
    return Piecewise(n, [truncmul(p1[i], p2[i]; kwargs..., dx = n[i+1]-n[i]) for i = 1:length(p1)])
end

for f in (:conj, :adjoint, :transpose, :real, :imag)
    @eval Base.$f(p::AbstractPiecewise) = Piecewise(nodes(p), map($f, elements(p)))
end

LinearAlgebra.tr(p::AbstractPiecewise) = Piecewise(nodes(p), map(tr, elements(p)))

# Arithmetic (in place / mutating methods)
function LinearAlgebra.rmul!(p::Piecewise, α)
    for i = 1:length(p)
        rmul!(p[i], α)
    end
    return p
end

function LinearAlgebra.lmul!(α, p::Piecewise)
    for i = 1:length(p)
        lmul!(α, p[i])
    end
    return p
end

function LinearAlgebra.mul!(pdst::Piecewise, α::Const, psrc::Piecewise)
    @assert nodes(pdst) == nodes(psrc)
    for i = 1:length(psrc)
        mul!(pdst[i], α, psrc[i])
    end
    return pdst
end

function LinearAlgebra.mul!(pdst::Piecewise, psrc::Piecewise, α::Const)
    @assert nodes(pdst) == nodes(psrc)
    for i = 1:length(psrc)
        mul!(pdst[i], psrc[i], α)
    end
    return pdst
end

function LinearAlgebra.axpy!(α, px::AbstractPiecewise, py::Piecewise)
    @assert nodes(px) == nodes(py)
    for i = 1:length(px)
        axpy!(α, px[i], py[i])
    end
    return py
end
function LinearAlgebra.axpby!(α, px::AbstractPiecewise, β, py::Piecewise)
    @assert nodes(px) == nodes(py)
    for i = 1:length(px)
        axpby!(α, px[i], β, py[i])
    end
    return py
end

function truncmul!(p::Piecewise, p1::AbstractPiecewise, p2::AbstractPiecewise,
                                α = true, β = false; kwargs...)
    @assert nodes(p) == nodes(p1) == nodes(p2)
    n = nodes(p)
    @inbounds for i = 1:length(p)
        truncmul!(p[i], p1[i], p2[i], α, β; kwargs..., dx = n[i+1]-n[i])
    end
    return p
end

function LinearAlgebra.mul!(p::Piecewise, p1::AbstractPiecewise, p2::AbstractPiecewise,
                                α = true, β = false)
    @assert nodes(p) == nodes(p1) == nodes(p2)
    @inbounds for i = 1:length(p)
        mul!(p[i], p1[i], p2[i], α, β)
    end
    return p
end

# Inner product and norm
function localdot(p1::AbstractPiecewise, p2::AbstractPiecewise)
    @assert nodes(p1) == nodes(p2)
    return Piecewise(nodes(p1), map(localdot, elements(p1), elements(p2)))
end

LinearAlgebra.dot(p1::AbstractPiecewise, p2::AbstractPiecewise) = integrate(localdot(p1, p2), domain(p1))

LinearAlgebra.norm(p::AbstractPiecewise) = sqrt(dot(p, p))

# Differentiate and integrate
differentiate(p::AbstractPiecewise) = Piecewise(nodes(p), map(differentiate, elements(p)))
function integrate(p::AbstractPiecewise, interval = domain(p))
    @assert interval == domain(p)
    s = integrate(p[1], (nodes(p)[1], nodes(p)[2]))
    for i = 2:length(p)
        s += integrate(p[i], (nodes(p)[i], nodes(p)[i+1]))
    end
    return s
end


text:
using TrajGWAS
using Test
using CSV, DataFrames, WiSER

filepath = normpath(joinpath(dirname(pathof(WiSER)), "../data/"))
df = DataFrame(CSV.File(filepath * "sbp.csv"), copycols = false)
vlmm = WSVarLmmModel(
    @formula(sbp ~ 1 + agegroup + bmi_std + meds), # removed obswt, gender
    @formula(sbp ~ 1 + bmi_std),
    @formula(sbp ~ 1 + agegroup + bmi_std), # removed obswt, gender, meds
    :id, df
);
WiSER.fit!(vlmm)
df_cut = df[!, [:obswt, :gender, :meds]]
df_cut[!, :gender_bool] = map(x -> x == "Male", df_cut[!, :gender])
df_cut[!, :meds_bool] = map(x -> x == "OnMeds", df_cut[!, :meds])
df_cut = df_cut[!, setdiff(names(df_cut), ["meds", "gender"])];
X1vec = Matrix{Float64}[]
W1vec = Matrix{Float64}[]
for i in unique(df[!, :id])
    push!(X1vec, Matrix{Float64}(df_cut[df[!, :id] .== i, :]))
    push!(W1vec, Matrix{Float64}(df_cut[df[!, :id] .== i, :]))
end

@testset "generic" begin
st = TrajGWAS.WSVarScoreTest(vlmm, 3, 3)
_, v1, _, _ = TrajGWAS.test!(st, X1vec, W1vec)
@test isapprox(v1, 6.236443105947029e-11)

st = TrajGWAS.WSVarScoreTest(vlmm, 1, 1);

# testing for obswt
X1vec_obswt = map(x -> x[:, 1], X1vec)
v1, v2, _, _ = TrajGWAS.test!(st, X1vec_obswt, X1vec_obswt)
@test all(isapprox.((v1, v2),
    (0.4409730648525017, 0.7156750808059059)))

# testing for gender
X1vec_gender = map(x -> x[:, 2], X1vec)
v1, v2, _, _ = TrajGWAS.test!(st, X1vec_gender, X1vec_gender)
@test isapprox(v2, 0.8535798925691855)

# testing for meds. meds is already included in the null model.
st = TrajGWAS.WSVarScoreTest(vlmm, 0, 1)
X1vec_meds = map(x -> x[:, 3], X1vec)
v1, v2, _, _ = TrajGWAS.test!(st, nothing, X1vec_meds)
@test all(isapprox.((v1, v2),
    (-1.0, 1.3935485985991808e-12)))

end

@testset "generic_invariant" begin
for i in 1:length(X1vec)
    for j in 1:size(X1vec[i], 1)
        @assert all(X1vec[i][j, :] .== X1vec[i][1, :]) # in fact, we are dealing with time-invariant features
    end
end
X1 = Matrix{Float64}(undef, length(X1vec), 3)
for i in 1:length(X1vec)
    X1[i, :] .= X1vec[i][1, :]
end
W1 = copy(X1);

st_i = TrajGWAS.WSVarScoreTestInvariant(vlmm, 3, 3)
_, v1, _, _ = TrajGWAS.test!(st_i, X1, W1)
@test isapprox(v1, 6.236443105947029e-11)

st_i = TrajGWAS.WSVarScoreTestInvariant(vlmm, 1, 1);
# obswt
v1, v2, _, _ = TrajGWAS.test!(st_i, X1[:, 1], W1[:, 1])
@test all(isapprox.((v1, v2),
    (0.4409730648525017, 0.7156750808059059)))

# testing for gender
v1, v2, _, _ = TrajGWAS.test!(st_i, X1[:, 2], W1[:, 2])
@test isapprox(v2, 0.8535798925691855)

# testing for meds. meds is already included in the null model.
# meds
st_i = TrajGWAS.WSVarScoreTestInvariant(vlmm, 0, 1);
v1, v2, _, _ = TrajGWAS.test!(st_i, nothing, W1[:, 3])
@test all(isapprox.((v1, v2),
    (-1.0, 1.3935485985991808e-12)))
end

filepath = joinpath(dirname(pathof(TrajGWAS)), "../data/")

pvalpath = "testpval"
@testset "trajgwas_singlesnp_plink" begin
trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_plinkex.csv",
        filepath * "hapmap3",
        pvalfile = pvalpath,
        usespa = false)
results = CSV.read(pvalpath, DataFrame)
@test all(isapprox.((mean(results.betapval),
    mean(results.taupval)),
    (0.20892290777343045, 0.1905578803635017)))

#spa
trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_plinkex.csv",
        filepath * "hapmap3",
        pvalfile = pvalpath,
        usespa = true)
results = CSV.read(pvalpath, DataFrame)

@test all(isapprox.((mean(results.betapval),
    mean(results.taupval)),
    (0.20841900621550116, 0.18913803849088423); rtol=1e-5))
end

@testset "trajgwas_singlesnp_vcf" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_vcfex.csv",
        filepath * "test_vcf",
        pvalfile = pvalpath,
        geneticformat = "VCF",
        vcftype = :DS,
        usespa = false)
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.47580382024875084, 0.4752535666464088); rtol=1e-5))

    #spa version
    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_vcfex.csv",
        filepath * "test_vcf",
        pvalfile = pvalpath,
        geneticformat = "VCF",
        vcftype = :DS,
        usespa = true)
    results = CSV.read(pvalpath, DataFrame)

    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.4709957279235137, 0.4691934843663542); rtol=1e-5))
end

@testset "trajgwas_singlesnp_bgen" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_bgen_ex.csv",
        filepath * "example.8bits",
        usespa = false,
        geneticformat = "BGEN",
        pvalfile = pvalpath,
        ref_dosage=false)
    results = CSV.read(pvalpath, DataFrame)

    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.46429181963000665, 0.4931759050587741); rtol=1e-5))

    #spa version
    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_bgen_ex.csv",
        filepath * "example.8bits",
        usespa = true,
        geneticformat = "BGEN",
        pvalfile = pvalpath,
        reportchisq=true,
        ref_dosage=false)#Adjustor(ones(500, 1)))
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.46428062017147764, 0.49230409589144597); rtol=1e-3))
end

@testset "trajgwas_snpset_plink" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
            @formula(y ~ 1),
            @formula(y ~ 1 + sex),
            :id,
            filepath * "trajgwas_plinkex.csv",
            filepath * "hapmap3",
            pvalfile = pvalpath,
            analysistype = "snpset",
            snpset = 2)
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.10426027515117509, 0.08863103713602266); rtol=1e-5))

    trajgwas(@formula(y ~ 1 + sex + onMeds),
            @formula(y ~ 1),
            @formula(y ~ 1 + sex),
            :id,
            filepath * "trajgwas_plinkex.csv",
            filepath * "hapmap3",
            pvalfile = pvalpath,
            analysistype = "snpset",
            snpset = filepath * "hapmap_snpsetfile.txt")
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.08547206177854592, 0.07865194590141143); rtol=1e-5))

    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex),
        :id,
        filepath * "trajgwas_plinkex.csv",
        filepath * "hapmap3",
        pvalfile = pvalpath,
        analysistype = "snpset",
        snpset = 1:5)
    pvalfile = open(pvalpath)
    pvals = split(readline(pvalfile), r"[, ]")[[end-3, end]]
    close(pvalfile)
    @test all(isapprox.(parse.(Float64, pvals), (9.948443168811026e-6,
    1.2760805146705195e-13), rtol=1e-4))
end

@testset "trajgwas_snpset_vcf" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
            @formula(y ~ 1),
            @formula(y ~ 1 + sex + onMeds),
            :id,
            filepath * "trajgwas_vcfex.csv",
            filepath * "test_vcf",
            pvalfile = pvalpath,
            geneticformat = "VCF",
            vcftype = :DS,
            analysistype = "snpset",
            snpset = 20)
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.25347423101832717, 0.19024751040999272); rtol=1e-5))

    trajgwas(@formula(y ~ 1 + sex + onMeds),
            @formula(y ~ 1),
            @formula(y ~ 1 + sex),
            :id,
            filepath * "trajgwas_vcfex.csv",
            filepath * "test_vcf",
            pvalfile = pvalpath,
            geneticformat = "VCF",
            vcftype = :DS,
            analysistype = "snpset",
            snpset = filepath * "snpsetfile_vcf.txt")
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.24940609907211, 0.22513112148296); rtol=1e-5))

    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex),
        :id,
        filepath * "trajgwas_vcfex.csv",
        filepath * "test_vcf",
        pvalfile = pvalpath,
        geneticformat = "VCF",
        vcftype = :DS,
        analysistype = "snpset",
        snpset = 1:5)
    pvalfile = open(pvalpath)
    pvals = split(readline(pvalfile), r"[, ]")[[end-3, end]]
    close(pvalfile)
    @test all(isapprox.(parse.(Float64, pvals), (0.079501633468,
    0.01825642314521), rtol=1e-4))
end

@testset "trajgwas_snpset_bgen" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_bgen_ex.csv",
        filepath * "example.8bits",
        geneticformat = "BGEN",
        pvalfile = pvalpath,
        analysistype = "snpset",
        snpset = 20,
        ref_dosage=false)
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.4444189887865223, 0.43837358376815266); rtol=1e-5))

    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex),
        :id,
        filepath * "trajgwas_bgen_ex.csv",
        filepath * "example.8bits",
        geneticformat = "BGEN",
        pvalfile = pvalpath,
        analysistype = "snpset",
        snpset = filepath * "bgen_snpsetfile.txt",
        ref_dosage=false,
        r=2.0,
        adjustor=Adjustor(ones(500, 1)))
        #adjustor=Adjustor(ones(500, 1)))
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.43077865459266, 0.38136082362710); rtol=1e-5))

    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex),
        :id,
        filepath * "trajgwas_bgen_ex.csv",
        filepath * "example.8bits",
        geneticformat = "BGEN",
        pvalfile = pvalpath,
        analysistype = "snpset",
        snpset = 1:5,
        ref_dosage=false)
    pvalfile = open(pvalpath)
    pvals = split(readline(pvalfile), r"[, ]")[[end-3, end]]
    close(pvalfile)
    @test all(isapprox.(parse.(Float64, pvals), (7.975567153147366e-7,
    7.722715354932478e-20), rtol=1e-4))
end

@testset "trajgwas_gxe_plink" begin
trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_plinkex.csv",
        filepath * "hapmap3",
        pvalfile = pvalpath,
        analysistype = "gxe",
        e = :sex,
        snpinds = 1:10)
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.5917252588216351, 0.4047112266291169); rtol=1e-5))
end

@testset "trajgwas_gxe_vcf" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
            @formula(y ~ 1),
            @formula(y ~ 1 + sex + onMeds),
            :id,
            filepath * "trajgwas_vcfex.csv",
            filepath * "test_vcf",
            pvalfile = pvalpath,
            geneticformat = "VCF",
            vcftype = :DS,
            analysistype = "gxe",
            e = :sex,
            snpinds = 1:10)
    results = CSV.read(pvalpath, DataFrame)
    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.6296884957250368, 0.711909935890595); rtol=1e-5))
    end

@testset "trajgwas_gxe_bgen" begin
    trajgwas(@formula(y ~ 1 + sex + onMeds),
        @formula(y ~ 1),
        @formula(y ~ 1 + sex + onMeds),
        :id,
        filepath * "trajgwas_bgen_ex.csv",
        filepath * "example.8bits",
        geneticformat = "BGEN",
        pvalfile = pvalpath,
        analysistype = "gxe",
        e = :sex,
        snpinds = 1:5)
    results = CSV.read(pvalpath, DataFrame)

    @test all(isapprox.((mean(results.betapval),
        mean(results.taupval)),
        (0.6339056748172345, 0.34364220119286115); rtol=1e-3))
    @test all(isapprox.((mean(results.snpeffectnullbeta),
        mean(results.snpeffectnulltau)),
        (-0.0144768768325708, 0.005777868071713646); rtol=1e-5))
end

rm(pvalpath)
rm("trajgwas.null.txt")


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/method/lohners_qr.jl
# Defines the lohner's method for contracting in discretize and relax.
#############################################################################

"""
LohnersFunctor

A functor used in computing bounds and relaxations via Lohner's method. The
implementation of the parametric Lohner's method described in the paper in (1)
based on the non-parametric version given in (2).

1. [Sahlodin, Ali M., and Benoit Chachuat. "Discretize-then-relax approach for
convex/concave relaxations of the solutions of parametric ODEs." Applied Numerical
Mathematics 61.7 (2011): 803-820.](https://www.sciencedirect.com/science/article/abs/pii/S0168927411000316)
2. [R.J. Lohner, Computation of guaranteed enclosures for the solutions of
ordinary initial and boundary value problems, in: J.R. Cash, I. Gladwell (Eds.),
Computational Ordinary Differential Equations, vol. 1, Clarendon Press, 1992,
pp. 425–436.](http://www.goldsztejn.com/old-papers/Lohner-1992.pdf)
"""
mutable struct LohnersFunctor{F <: Function, K, T <: Real, S <: Real, NY} <: AbstractStateContractor
    set_tf!::TaylorFunctor!{F, K, T, S}
    real_tf!::TaylorFunctor!{F, K, T, T}
    jac_tf!::JacTaylorFunctor!{F, K, T, S, NY}
    η::Interval{Float64}
    μX::Vector{S}
    ρP::Vector{S}
    f̃::Vector{Vector{S}}
    f̃val::Vector{Vector{Float64}}
    Vⱼ₊₁::Vector{Float64}
    Rⱼ₊₁::Vector{S}
    mRⱼ₊₁::Vector{Float64}
    Δⱼ₊₁::Vector{S}
    Jxmat::Matrix{S}
    Jxvec::Vector{S}
    Jpvec::Vector{S}
    rRⱼ₊₁::Vector{S}
    YdRⱼ₊₁::Vector{S}
    YJxmat::Matrix{S}
    YJxvec::Vector{S}
    YJpmat::Matrix{S}
    YJpvec::Vector{S}
end
function LohnersFunctor(f!::F, nx::Int, np::Int, k::Val{K}, s::S, t::T) where {F, K, S <: Number, T <: Number}
    set_tf! = TaylorFunctor!(f!, nx, np, k, zero(S), zero(T))
    real_tf! = TaylorFunctor!(f!, nx, np, k, zero(T), zero(T))
    jac_tf! = JacTaylorFunctor!(f!, nx, np, k, zero(S), zero(T))
    μX = zeros(S, nx)
    ρP = zeros(S, np)
    f̃ = Vector{S}[]
    f̃val = Vector{Float64}[]
    for i = 1:(K+1)
        push!(f̃, zeros(S, nx))
        push!(f̃val, zeros(Float64, nx))
    end
    Vⱼ₊₁ = zeros(Float64, nx)
    Rⱼ₊₁ = zeros(S, nx)
    mRⱼ₊₁ = zeros(Float64, nx)
    Δⱼ₊₁ = zeros(S, nx)
    Jxmat = zeros(S, nx, nx)
    Jxvec = zeros(S, nx)
    Jpvec = zeros(S, nx)
    rRⱼ₊₁ = zeros(S, nx)
    YdRⱼ₊₁ = zeros(S, nx)
    YJxmat = zeros(S, nx, nx)
    YJxvec = zeros(S, nx)
    YJpmat = zeros(S, nx, np)
    YJpvec = zeros(S, nx)
    LohnersFunctor{F, K+1, T, S, nx+np}(set_tf!, real_tf!, jac_tf!, Interval{Float64}(0.0,1.0), μX, ρP,
                                        f̃, f̃val, Vⱼ₊₁, Rⱼ₊₁, mRⱼ₊₁, Δⱼ₊₁, Jxmat, Jxvec, Jpvec,
                                        rRⱼ₊₁, YdRⱼ₊₁, YJxmat, YJxvec, YJpmat, YJpvec)
end

"""
LohnerContractor{K}

An `AbstractStateContractorName` used to specify a parametric Lohners method contractor
of order K.
"""
struct LohnerContractor{K} <: AbstractStateContractorName end
LohnerContractor(k::Int) = LohnerContractor{k}()
function state_contractor(m::LohnerContractor{K}, f, Jx!, Jp!, nx, np, style, s, h) where K
    LohnersFunctor(f, nx, np, Val{K}(), style, s)
end

state_contractor_k(m::LohnerContractor{K}) where K = K
state_contractor_γ(m::LohnerContractor{K}) where K = 1.0
state_contractor_steps(m::LohnerContractor{K}) where K = 2
state_contractor_integrator(m::LohnerContractor{K}) where K = CVODE_Adams()

function (d::LohnersFunctor{F,K,S,T,NY})(contract::ContractorStorage{T},
                                         result::StepResult{T},
                                         count::Int) where {F <: Function, K, S <: Real, T <: Real, NY}

    set_tf! = d.set_tf!
    real_tf! = d.real_tf!
    Jf! = d.jac_tf!

    k = set_tf!.k
    nx = set_tf!.nx
    np = length(contract.pval)

    hⱼ = contract.hj_computed
    hjk = hⱼ^k
    t = contract.times[1]

    set_tf!(d.f̃, contract.Xj_apriori, contract.P, t)

    # computes Rj and it's midpoint
    for i = 1:nx
        @inbounds d.Rⱼ₊₁[i] = hjk*d.f̃[k + 1][i]
        @inbounds d.mRⱼ₊₁[i] = mid(d.Rⱼ₊₁[i])
    end

    # defunes new x point... k corresponds to k - 1 since taylor
    # coefficients are zero indexed
    real_tf!(d.f̃val, contract.xval, contract.pval, t)
    hji1 = 0.0
    fill!(d.Vⱼ₊₁, 0.0)
    for i = 1:k
        hji1 = hⱼ^i
        @__dot__ d.Vⱼ₊₁ += hji1*d.f̃val[i + 1]
    end
    @__dot__ d.Vⱼ₊₁ += contract.xval

    # compute extensions of taylor cofficients for rhs
    μ!(d.μX, contract.Xj_0, contract.xval, d.η)
    ρ!(d.ρP, contract.P, contract.pval, d.η)
    set_JxJp!(Jf!, d.μX, d.ρP, t)
    for i = 1:k
        hji1 = hⱼ^(i - 1)
        if i == 1
            fill!(Jf!.Jxsto, zero(S))
            for j = 1:nx
                Jf!.Jxsto[j,j] = one(S)
            end
        else
            @__dot__ Jf!.Jxsto += hji1*Jf!.Jx[i]
        end
        @__dot__ Jf!.Jpsto += hji1*Jf!.Jp[i]
    end

    # update x floating point value
    @__dot__ contract.xval_computed = d.Vⱼ₊₁ + d.mRⱼ₊₁

    # update bounds on X at new time
    mul_split!(d.Jxmat, Jf!.Jxsto, contract.A_Q[2], nx)
    mul_split!(d.Jxvec, d.Jxmat, contract.Δ[2], nx)
    mul_split!(d.Jpvec, Jf!.Jpsto, contract.rP, nx)

    @__dot__ contract.X_computed = d.Vⱼ₊₁ + d.Rⱼ₊₁ + d.Jxvec + d.Jpvec
    affine_contract!(contract.X_computed, contract.P, contract.pval, nx, np)

    # calculation block for computing Aⱼ₊₁ and inv(Aⱼ₊₁)
    @__dot__ contract.B = mid(d.Jxmat)
    calculateQ!(contract.A_Q[1], contract.B, nx)
    calculateQinv!(contract.A_inv[1], contract.A_Q[1], nx)

    # update Delta
    @__dot__ d.rRⱼ₊₁ = d.Rⱼ₊₁ - d.mRⱼ₊₁
    mul_split!(d.YdRⱼ₊₁, contract.A_inv[1], d.rRⱼ₊₁, nx)
    mul_split!(d.YJxmat, contract.A_inv[1], d.Jxmat, nx)
    mul_split!(d.YJxvec, d.YJxmat, contract.Δ[2], nx)
    mul_split!(d.YJpmat, contract.A_inv[1], Jf!.Jpsto, nx)
    mul_split!(d.YJpvec, d.YJpmat, contract.rP, nx)
    @__dot__ d.Δⱼ₊₁ = d.YdRⱼ₊₁ + d.YJxvec + d.YJpvec

    contract.Δ[1] = copy(d.Δⱼ₊₁)

    return RELAXATION_NOT_CALLED
end

get_Δ(lf::LohnersFunctor{F,K,S,T,NY}) where {F <: Function, K, S <: Real, T <: Real, NY} = copy(lf.Δⱼ₊₁)
advance_contractor_buffer!(lf::LohnersFunctor) = nothing


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

nearestpsd_insts(use_nat::Bool) = [
    [(side, use_completable, false, use_nat, true)
    for side in vcat(30, 50:50:200, 300:100:800)] # includes compile run
    for use_completable in (false, true)
    ]

insts = OrderedDict()
insts["nat"] = (nothing, nearestpsd_insts(true))
insts["ext"] = (:SOCExpPSD, nearestpsd_insts(false))
return (NearestPSDJuMP, insts)


text:
# # File I/O
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/fileio.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/fileio.ipynb)

# ## Import Data

using RHEOS

# RHEOS has a convenience function for importing data from CSV files: [`importcsv`](@ref). The default column delimiter is `,` but an alternative can be specified as a keyword argument. The row delimiter is a newline character (`\n`). For standard time-domain viscoelastic testing data RHEOS expects either stress, strain and time data, just stress and time, or just strain and time. Arguments must be identified by providing the number of the column in which they are contained. The function returns a [`RheoTimeData`](@ref) object.

## Import file
data_1 = importcsv("assets/data_time.csv", t_col = 1, ϵ_col = 2, σ_col = 3)
## Check data type
rheotimedatatype(data_1)

## Import incomplete data
data_2 = importcsv("assets/data_time.csv", t_col = 1, ϵ_col = 2)
## Check data type
rheotimedatatype(data_2)

# The function [`importcsv`](@ref) can also be used to import frequency, storage modulus, and loss modulus data (as a complete set). In this case, the function returns a [`RheoFreqData`](@ref) object.

## Import file
data_f = importcsv("assets/data_freq.csv", ω_col = 1, Gp_col = 2, Gpp_col = 3)
## Check data type
rheofreqdatatype(data_f)

# ## Export Data

# If you want to analyse or plot your data in software other than Julia you will likely want to export it to a CSV file. To export [`RheoTimeData`](@ref) and [`RheoFreqData`](@ref) objects to CSV files we can use the [`exportcsv`](@ref) function. For the two complete data-sets we imported above, we can export them into new files in the following way. As with [`importcsv`](@ref), the order of the columns can be specified by the user.

## Export file
exportcsv(data_1, "assets/export_timedata.csv")
exportcsv(data_f, "assets/export_frequdata.csv")


text:
mutable struct SoysterMaximinModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::SoysterMaximinObjective
    output::Output
    SoysterMaximinModel(settings, constraints, obj, output) =
        new(settings, constraints, obj, output)
    SoysterMaximinModel() =
        new(Settings(), Constraint[], SoysterMaximinObjective(), Output())
end


text:
file_dir = Pkg.dir("D4M")*"/examples/2Apps/1EntityAnalysis/Entity.csv";
E = ReadCSV(file_dir);
print(E[1:5,:]);

~,~,doc_val      = find(E[:,"doc,"]);
~,~,entity_val   = find(E[:,"entity,"]);
~,~,position_val = find(E[:,"position,"]);
~,~,type_val     = find(E[:,"type,"]);

typeEntity_val = CatStr(type_val, "/" , entity_val);

E = Assoc(doc_val,typeEntity_val,position_val);

using JLD
save(Pkg.dir("D4M")*"/examples/2Apps/1EntityAnalysis/Entity.jld" , "E",E );

print(E[1,:]);

using PyPlot
spy(E');


text:
using Test
using QuantumOptics
using LinearAlgebra, SparseArrays, Random

mutable struct SpectralanalysisTestOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR} end

@testset "spectralanalysis" begin

Random.seed!(0)

sprandop(b) = sparse(DenseOperator(b, rand(ComplexF64, length(b), length(b))))

# Test diagonalization
@test_throws ArgumentError eigenstates(SpectralanalysisTestOperator{Basis,Basis}())
@test_throws QuantumOpticsBase.IncompatibleBases eigenstates(DenseOperator(GenericBasis(3), GenericBasis(4)))
@test_throws ArgumentError eigenenergies(SpectralanalysisTestOperator{Basis,Basis}())
@test_throws QuantumOpticsBase.IncompatibleBases eigenenergies(DenseOperator(GenericBasis(3), GenericBasis(4)))

# Test hermitian diagonalization
b = GenericBasis(5)
a = randoperator(b)
H = (a+dagger(a))/2
U = exp(1im*H)
d = [-3, -2.6, -0.1, 0.0, 0.6]
D = DenseOperator(b, diagm(0 => d))
Dsp = sparse(D)
@test eigenenergies(D) ≈ d
@test eigenenergies(Dsp, 3, info=false) ≈ d[1:3]

R = U*D*dagger(U)
Rsp = sparse(R)
@test eigenenergies((R+dagger(R))/2) ≈ d
@test eigenenergies((Rsp+dagger(Rsp))/2, 3; info=false) ≈ d[1:3]

E, states = eigenstates((R+dagger(R))/2, 3)
Esp, states_sp = eigenstates((Rsp+dagger(Rsp))/2, 3, info=false)
for i=1:3
    @test E[i] ≈ d[i]
    @test Esp[i] ≈ d[i]
    v = U.data[1,i]/states[i].data[1]
    @test states[i].data*v ≈ U.data[:,i]
    v = U.data[1,i]/states_sp[i].data[1]
    @test states_sp[i].data*v ≈ U.data[:,i]
end

# Test nonhermitian diagonalization
b = GenericBasis(5)
a = randoperator(b)
H = (a+dagger(a))/2
U = exp(1im*H)
d = [-3+0.2im, -2.6-0.1im, -0.1+0.5im, 0.0, 0.6+0.3im]
D = DenseOperator(b, diagm(0 => d))
Dsp = sparse(D)
@test eigenenergies(D; warning=false) ≈ d
@test eigenenergies(Dsp, 3; warning=false, info=false) ≈ d[1:3]

R = U*D*dagger(U)
Rsp = sparse(R)
@test eigenenergies(R; warning=false) ≈ d
@test eigenenergies(Rsp, 3; warning=false, info=false) ≈ d[1:3]

E, states = eigenstates(R, 3; warning=false)
Esp, states_sp = eigenstates(Rsp, 3; warning=false, info=false)
for i=1:3
    @test E[i] ≈ d[i]
    @test Esp[i] ≈ d[i]
    v = U.data[1,i]/states[i].data[1]
    @test states[i].data*v ≈ U.data[:,i]
    v = U.data[1,i]/states_sp[i].data[1]
    @test states_sp[i].data*v ≈ U.data[:,i]
end

# Test simdiag
spinbasis = SpinBasis(1//2)
sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
twospinbasis = spinbasis ⊗ spinbasis
Sx = dense(sum([embed(twospinbasis, i, sx) for i=1:2]))/2.
Sy = dense(sum([embed(twospinbasis, i, sy) for i=1:2]))/2.
Sz = dense(sum([embed(twospinbasis, i, sz) for i=1:2]))/2.
Ssq = Sx^2 + Sy^2 + Sz^2
d, v = simdiag([Sz, Ssq])
@test d[1] == [-1.0, 0, 0, 1.0]
@test d[2] ≈ [2, 0.0, 2, 2]
@test_throws ErrorException simdiag([Sx, Sy])

threespinbasis = spinbasis ⊗ spinbasis ⊗ spinbasis
Sx3 = dense(sum([embed(threespinbasis, i, sx) for i=1:3])/2.)
Sy3 = dense(sum([embed(threespinbasis, i, sy) for i=1:3])/2.)
Sz3 = dense(sum([embed(threespinbasis, i, sz) for i=1:3])/2.)
Ssq3 = Sx3^2 + Sy3^2 + Sz3^2
d3, v3 = simdiag([Ssq3, Sz3])
dsq3_std = eigenenergies(dense(Ssq3))
@test diagm(0 => dsq3_std) ≈ v3'*Ssq3.data*v3

fockbasis = FockBasis(4)
@test_throws ErrorException simdiag([Sy3, Sz3])
@test_throws ErrorException simdiag([dense(destroy(fockbasis)), dense(create(fockbasis))])

end # testset


text:
using SymPyCall
using LinearAlgebra
using Test

## SymPyCall only is adding methods for those found in SpecialFunctions
## so not the special polynomials, etc.
using SpecialFunctions

@testset "SpecialFuns" begin
    @syms a b x y
    @syms n m θ ϕ
    @syms ν::integer

    @test sympy.fresnels(Sym(0)) == 0
    @test sympy.fresnels(Sym(oo)) == Sym(1)/2
    @test diff(sympy.fresnels(x), x) == sin(PI*x^2/2)
    #@test evalf(fresnels(Sym(2)), 30) == Sym(parse(BigFloat, "0.343415678363698242195300815958"))

    @test  diff(sympy.fresnelc(x), x) == cos(PI*x^2/2)

    @test sympy.Ei(Sym(-1)) == sympy.Ei(exp(1im*PI))
    @test diff(sympy.Ei(x), x) == exp(x)/x
    @test diff(sympy.Si(x), x) == sin(x)/x
    @test diff(sympy.Ci(x), x) == cos(x)/x

    @test airyai(Sym(0)) == 3^(Sym(1)/3)/(3*gamma(Sym(2)/3))
    @test airybi(Sym(0)) == 3^(Sym(5)/6)/(3*gamma(Sym(2)/3))

    @test sympy.jacobi(0, a, b, x) == 1
    @test sympy.jacobi(1, a, b, x) == a/2 - b/2 + x*(a/2 + b/2 + 1)
    @test sympy.jacobi(n, 0, 0, x) == sympy.legendre(n, x)
    @test sympy.jacobi(n, a, b, -x) == (-1)^n*sympy.jacobi(n, b, a, x)
    @test conj(sympy.jacobi(n, a, b, x)) == sympy.jacobi(n, conj(a), conj(b), conj(x))

    @test sympy.gegenbauer(0, a, x) == 1
    @test sympy.gegenbauer(1, a, x) == 2a*x
    @test sympy.gegenbauer(2, a, x) == -a + x^2*(2a^2 + 2a)
    @test conj(sympy.gegenbauer(n, a, x)) == sympy.gegenbauer(n, conj(a), conj(x))
    @test diff(sympy.gegenbauer(n, a, x), x) == 2a*sympy.gegenbauer(n-1, a+1, x)

    @test sympy.chebyshevt(0, x) == 1
    @test sympy.chebyshevt(1, x) == x
    @test sympy.chebyshevt(2, x) == 2x^2-1

    @test sympy.chebyshevu(0, x) == 1
    @test sympy.chebyshevu(1, x) == 2x
    @test sympy.chebyshevu(2, x) == 4x^2 - 1
    @test sympy.chebyshevu(n, 0) == cos(π*n/2)
    @test sympy.chebyshevu(n, 1) == n + 1

    @test sympy.legendre(0, x) == 1
    @test sympy.legendre(1, x) == x
    @test sympy.legendre(2, x) == 3x^2/2 - 1//2
    @test sympy.legendre(n, x) == sympy.legendre(n, x)
    @test diff(sympy.legendre(n,x), x) == n*(x*sympy.legendre(n, x) - sympy.legendre(n - 1, x))/(x^2 - 1)

    @test sympy.assoc_legendre(0,0, x) == 1
    @test sympy.assoc_legendre(1,0, x) == x
    @test sympy.assoc_legendre(1,1, x) == -sqrt(-x^2 + 1)
    #@test assoc_legendre(n,m,x) == assoc_legendre(n, m, x)

    @test sympy.hermite(0, x) == 1
    @test sympy.hermite(1, x) == 2x
    @test sympy.hermite(2, x) == 4x^2 - 2
    @test sympy.hermite(n, x) == sympy.hermite(n, x)
    @test diff(sympy.hermite(n,x), x) == 2n*sympy.hermite(n - 1, x)
    @test sympy.hermite(n, -x) == (-1)^n*sympy.hermite(n, x)

    @test sympy.laguerre(0, x) == 1
    @test sympy.laguerre(1, x) == -x + 1
    @test sympy.laguerre(2, x) == x^2/2 - 2x + 1
    @test sympy.laguerre(3, x) == -x^3/6 + 3x^2/2 - 3x + 1
    @test diff(sympy.laguerre(n, x), x) == -sympy.assoc_laguerre(n - 1, 1, x)

    @test sympy.assoc_laguerre(0, a, x) == 1
    @test sympy.assoc_laguerre(1, a, x) == a - x + 1
    @test sympy.assoc_laguerre(2, a, x) == a^2/2 + 3a/2 + x^2/2 + x*(-a - 2) + 1
    @test sympy.assoc_laguerre(n, 0, x) == sympy.laguerre(n, x)
    @test diff(sympy.assoc_laguerre(n, a, x), x) == -sympy.assoc_laguerre(n - 1, a + 1, x)

    @test sympy.Ynm(n, m, θ, ϕ) == sympy.Ynm(n, m, θ, ϕ)
    @test sympy.Ynm(n, -m, θ, ϕ) == (-1)^m*exp(-2im*m*ϕ)*sympy.Ynm(n, m, θ, ϕ)
    @test sympy.Ynm(n, m, -θ, ϕ) == sympy.Ynm(n, m, θ, ϕ)
    @test sympy.Ynm(n, m, θ, -ϕ) == exp(-2im*m*ϕ)*sympy.Ynm(n, m, θ, ϕ)
    @test expand(simplify(sympy.Ynm(0, 0, θ, ϕ)), func=true) == 1/(2*sqrt(PI))

    @test sympy.Ynm_c(n, m, θ, ϕ) == conj(sympy.Ynm(n, m, θ, ϕ))

    @test diff(hankelh1(n, x), x) == hankelh1(n - 1, x)/2 - hankelh1(n + 1, x)/2
    @test diff(hankelh2(n, x), x) == hankelh2(n - 1, x)/2 - hankelh2(n + 1, x)/2

    @test sympy.expand_func(sympy.jn(0, x)) == sin(x)/x
    @test sympy.expand_func(sympy.jn(1, x)) == sin(x)/x^2 - cos(x)/x
    @test rewrite(sympy.jn(ν, x), "besselj") == sqrt(2PI/x)*besselj(ν + Sym(1)/2, x)/2
    VERSION < v"0.6.0-dev" && @test rewrite(sympy.jn(ν, x), "bessely") == (-1)^ν*sqrt(2PI/x)*bessely(-ν - Sym(1)/2, x)/2
    u = N(sympy.jn(2, 5.2+0.3im).evalf(20))
    @test norm(real(u) - 0.099419756723640344491) <= 1e-15 && norm(imag(u) + 0.054525080242173562897) <= 1e-15

    @test sympy.expand_func(sympy.yn(0, x)) == -cos(x)/x
    @test sympy.expand_func(sympy.yn(1, x)) == -cos(x)/x^2-sin(x)/x
    VERSION < v"0.6.0-dev" && @test rewrite(sympy.yn(ν, x), "besselj") == (-1)^(ν + 1)*sqrt(2PI/x)*besselj(-ν - Sym(1)/2, x)/2
    @test rewrite(sympy.yn(ν, x), "bessely") == sqrt(2PI/x)*bessely(ν + Sym(1)/2, x)/2
    @test N(sympy.yn(2, 5.2+0.3im)) ≈ 0.18525034196069722536 + 0.014895573969924817587im

    # gamma, beta and related functions
    @test gamma(Sym(1)) == 1
    @test gamma(Sym(3)/2) == sqrt(PI)/2

    @test diff(polygamma(Sym(0), x), x) == polygamma(Sym(1), x)
    @test diff(polygamma(Sym(0), x), x, 2) == polygamma(Sym(2), x)

    @test diff(beta(x, y), x) == (polygamma(Sym(0), x) - polygamma(Sym(0), x + y)) * beta(x, y)
    @test diff(beta(x, y), y) == (polygamma(Sym(0), y) - polygamma(Sym(0), x + y)) * beta(x, y)

    # test numerical consistency with Julia functions
    @test N(gamma(Sym(4.1))) ≈ gamma(4.1)
    @test N(polygamma(Sym(2), Sym(3.2))) ≈ polygamma(2, 3.2)
    VERSION >= v"0.5.0" && @test N(beta(Sym(1)+1im, Sym(1)+1im)) ≈ beta(1.0+1im, 1.0+1im)

    # Elliptic-type functions

    @test sympy.elliptic_k(Sym(0)) == PI/2
    @test N(sympy.elliptic_k(Sym(1.0 + im))) ≈ 1.50923695405127 + 0.625146415202697im

    @test N(sympy.elliptic_f(Sym(3.0 + im/2), Sym(1.0 + im))) ≈ 2.909449841483 + 1.74720545502474im

    @test sympy.elliptic_e(Sym(0)) == PI/2
    @test N(sympy.elliptic_e(Sym(2.0 - im))) ≈ 0.991052601328069 + 0.81879421395609im

    @test sympy.elliptic_pi(Sym(0), Sym(0)) == PI/2
    @test N(sympy.elliptic_pi(Sym(1.0 - im/3), Sym(2.0 + im))) ≈ 3.29136443417283 + 0.32555634906645im

    # Bessel-type functions
    @test diff(besselj(n, x), x) == (besselj(n - 1, x) - besselj(n + 1, x))/2
    @test rewrite(besselj(n, x), "jn") == sqrt(2x/PI)*sympy.jn(n - 1//2, x)

    @test diff(bessely(n, x), x) == (bessely(n - 1, x) - bessely(n + 1, x))/2
    @test rewrite(bessely(n, x), "yn") == sqrt(2x/PI)*sympy.yn(n - 1//2, x)

    @test diff(besseli(n, x), x) == (besseli(n - 1, x) + besseli(n + 1, x))/2

    @test diff(besselk(n, x), x) == -(besselk(n - 1, x) + besselk(n + 1, x))/2

    # test numerical consistency with Julia functions
    @test N(besselj(Sym(3.2), Sym(1.5))) ≈ besselj(3.2, 1.5)
    @test N(bessely(Sym(3.2), Sym(1.5))) ≈ bessely(3.2, 1.5)
    @test N(besseli(Sym(3.2), Sym(1.5))) ≈ besseli(3.2, 1.5)
    @test N(besselk(Sym(3.2), Sym(1.5))) ≈ besselk(3.2, 1.5)

    @test sympy.expand_func(x*sympy.hyper([1, 1], [2], -x)) == log(x + 1)
end

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
using CUDA
using PaddedViews

function k_index_kernel_fft_3d(data, out, offx, offy, offz)
    dx = (blockIdx().x-1)*blockDim().x + threadIdx().x
    dy = (blockIdx().y-1)*blockDim().y + threadIdx().y
    dz = (blockIdx().z-1)*blockDim().z + threadIdx().z

    sx,sy,sz = size(data)
    sx = sx - offx
    sy = sy - offy
    sz = sz - offz
    if dx > sx || dy > sy || dz > sz
        return nothing
    end

    @inbounds out[dx,dy,dz] = data[dx+offx, dy+offy, dz + offz]

    return nothing
end

function k_zeroing_3d(data, offx, offy, offz)
    sx,sy,sz = size(data)
    dx = (blockIdx().x-1)*blockDim().x + threadIdx().x
    dy = (blockIdx().y-1)*blockDim().y + threadIdx().y
    dz = 1
    for dz = 1:offx
        if dx <= sx && dy <= sy && dz <= offz
            @inbounds data[dx, dy, end-dz+1] = 0
        end
        if dx <= sx && dy <= sz && dz <= offy
            @inbounds data[dx, end-dz+1, dy] = 0
        end
        if dx <= sz && dy <= sy && dz <= offx
            @inbounds data[end-dz+1, dy, dx] = 0
        end
    end
    return nothing
end

function cu_offset_result_3d(data, out, offx, offy, offz)
    dims = size(data)
    tx = 16; ty = 8; tz = 8
    bx = 1 + dims[1]÷tx
    by = 1 + dims[2]÷ty
    bz = 1 + dims[3]÷tz
    @cuda blocks=(bx,by, bz) threads=(tx,ty, tz) k_index_kernel_fft_3d(data,out,offx,offy,offz)
    tx = 32
    ty = 32
    #tz = 2max(offy,offx)
    bx = 1 + max(dims...)÷tx
    by = 1 + max(dims...)÷ty
    bz = 1
    #println("BLOCKS THREADS ",(bx,by,bz) ,(tx,ty,tz) , size(out), typeof(out))
    @cuda blocks=(bx,by) threads=(tx,ty) k_zeroing_3d(out,2offx,2offy,2offz)
end

function k_elem_mul_3d!(A, B)
    tx = (blockIdx().x-1)*blockDim().x + threadIdx().x
    ty = (blockIdx().y-1)*blockDim().y + threadIdx().y
    tz = (blockIdx().z-1)*blockDim().z + threadIdx().z
    sx,sy,sz = size(A)
    if tx > sx || ty > sy || tz > sz
        return nothing
    end
    @inbounds a = A[tx,ty,tz]
    @inbounds b = B[tx,ty,tz]
    @inbounds A[tx,ty,tz] = a*b
    return nothing
end

function elem_mul_3d!(A,B)
    sx,sy,sz = size(A)
    tx = 16; ty = 8; tz = 8

    bx = 1 + sx÷tx
    by = 1 + sy÷ty
    bz = 1 + sz÷tz
    @cuda blocks=(bx,by,bz) threads=(tx,ty,tz) k_elem_mul_3d!(A,B)
    nothing
end

function k_cmplx_to_real_3d!(A)
    dx = (blockIdx().x-1)*blockDim().x + threadIdx().x
    dy = (blockIdx().y-1)*blockDim().y + threadIdx().y
    dz = (blockIdx().z-1)*blockDim().z + threadIdx().z
    sx,sy,sz = size(A)
    if dx > sx || dy > sy || dz > sz
        return nothing
    end
    @inbounds a = A[dx,dy,dz]
    @inbounds A[dx,dy,dz] = a.re
    return nothing
end

function cmplx_to_real_3d!(A)
    sx,sy,sz = size(A)
    tx = 16; ty = 8; tz = 8
    bx = 1 + sx÷tx
    by = 1 + sy÷ty
    bz = 1 + sz÷tz
    @cuda blocks=(bx,by,bz) threads=(tx,ty,tz) k_cmplx_to_real_3d!(A)
    nothing
end

#
# function fft_stencil_3d(data, template, t_steps=1)
#     # template is a 3d square array
#     r = (size(template,1)-1)÷2
#     dx,dy,dz = size(data)
#     padded_data = PaddedView(0, data, (1:dx+2r, 1:dy+2r, 1:dz+2r), (1:dx,1:dy,1:dz))
#     padded_template = PaddedView(0, template, (1:dx+2r, 1:dy+2r, 1:dz+2r), (1:2r+1,1:2r+1,1:2r+1))
#     ## Initial Upload
#     fd = CuArray{CUDA.cuFloatComplex}(padded_data)
#     ft = CuArray{CUDA.cuFloatComplex}(padded_template)
#     temp = similar(ft)
#     CUFFT.fft!(fd)
#     CUFFT.fft!(ft)
#     elem_mul_3d!(fd,ft)
#     ## Time Loop
#     CUFFT.ifft!(fd)
#     cu_offset_result_3d(fd, temp, r,r,r)
#     cmplx_to_real_3d!(temp)
#     fd,temp = temp,fd
#     for t = 2:t_steps
#         CUFFT.fft!(fd)
#         elem_mul_3d!(fd,ft)
#         CUFFT.ifft!(fd)
#         cu_offset_result_3d(fd, temp, r,r,r)
#         cmplx_to_real_3d!(temp)
#         fd,temp = temp,fd
#     end
#     return view(real(Array(fd)), 1:dx, 1:dy, 1:dz)
# end

function fft_stencil_3d(data, template, t_steps=1)
    # template is a 3d square array
    #r = (size(template,1)-1)÷2

    r2 = size(template,1)
    r = 32
    while r <= r2
        r += 32
    end
    r = r + 1
    dx,dy,dz = size(data)
    padded_data = PaddedView(0, data, (1:dx+r-1, 1:dy+r-1, 1:dz+r-1), (1:dx,1:dy,1:dz))
    padded_template = PaddedView(0, template, (1:dx+r-1, 1:dy+r-1, 1:dz+r-1), (1:r2,1:r2,1:r2))
    ## Initial Upload
    fd = CuArray{CUDA.cuFloatComplex}(padded_data)
    ft = CuArray{CUDA.cuFloatComplex}(padded_template)
    temp = similar(ft)
    CUFFT.fft!(fd)
    CUFFT.fft!(ft)
    #@time m = fd.*ft
    elem_mul_3d!(fd,ft)
    ## Time Loop
    CUFFT.ifft!(fd)
    r2 = r2÷2 + 1
    rn = r2:(r2+dx-1)
    return real(Array(fd))[rn,rn,rn]
end


text:
struct InterfaceNotImplemented <: Exception
    family::Symbol
    method::String
end

Base.showerror(io::IO, err::InterfaceNotImplemented) =
    print(io, "Missing method from $(err.family) interface: `$(err.method)`")

struct InfiniteOrder{T} <: Exception
    x::T
    msg
    InfiniteOrder(g::Union{GroupElement, Group}) = new{typeof(g)}(g)
    InfiniteOrder(g::Union{GroupElement, Group}, msg) = new{typeof(g)}(g, msg)
end

function Base.showerror(io::IO, err::InfiniteOrder{T}) where T
    println(io, "Infinite order exception with ", err.x)
    if isdefined(err, :msg)
        print(io, err.msg)
    else
        print(io, "order will only return a value when it is finite. ")
        f = if T <: Group
            "isfinite(G)"
        elseif T <: GroupElement
            "isfiniteorder(g)"
        end
        print(io, "You should check with `$f` first.")
    end
end


text:
println("Eb) Simulation of the four-step cascade after 1s-3p photo-excitation of Si^+: Configuration model.")

@warn("\n\n !!! This example does not work properly at present !!! \n\n")
#
using JLD
wr = load("zzz-Cascade-2019-04-26T08.jld");   data = wr["cascade data:"];
wa = Cascade.Simulation([Cascade.IonDist(), Cascade.FinalDist()], Cascade.ProbPropagation(), Cascade.SimulationSettings() )
wb = perform(wa, data; output=false)

nothing


text:
#############################################
# Get the braking points
#############################################
function get_edges(X::Matrix{T}, nbins=250) where {T}
    edges = Vector{Vector{T}}(undef, size(X,2))
    @threads for i in 1:size(X, 2)
        edges[i] = unique(quantile(view(X, :,i), (0:nbins)/nbins))[2:end]
        if length(edges[i]) == 0
            edges[i] = [minimum(view(X, :,i))]
        end
    end
    return edges
end

####################################################
# Transform X matrix into a UInt8 binarized matrix
####################################################
function binarize(X, edges)
    X_bin = zeros(UInt8, size(X))
    @threads for i in 1:size(X, 2)
        X_bin[:,i] = searchsortedlast.(Ref(edges[i][1:end-1]), view(X,:,i)) .+ 1
    end
    X_bin
end

function find_bags(x_bin::Vector{T}) where T <: Real
    𝑖 = 1:length(x_bin) |> collect
    bags = [BitSet() for _ in 1:maximum(x_bin)]
    for bag in 1:length(bags)
        bags[bag] = BitSet(𝑖[x_bin .== bag])
    end
    return bags
end

function update_bags!(bins, set)
    for bin in bins
        intersect!(bin, set)
    end
end

function find_split_static!(hist_δ::Vector{SVector{L,T}}, hist_δ²::Vector{SVector{L,T}}, hist_𝑤::Vector{SVector{1,T}}, bins::Vector{BitSet}, X_bin, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, params::EvoTreeRegressor, info::SplitInfo{L,T,S}, edges::Vector{T}, set::Vector{S}) where {L,T,S}

    # initialize histogram
    hist_δ .*= 0.0
    hist_δ² .*= 0.0
    hist_𝑤 .*= 0.0

    # initialize tracking
    ∑δL = ∑δ * 0
    ∑δ²L = ∑δ² * 0
    ∑𝑤L = ∑𝑤 * 0
    ∑δR = ∑δ
    ∑δ²R = ∑δ²
    ∑𝑤R = ∑𝑤

    # build histogram
    @inbounds for i in set
        hist_δ[X_bin[i]] += δ[i]
        hist_δ²[X_bin[i]] += δ²[i]
        hist_𝑤[X_bin[i]] += 𝑤[i]
    end

    @inbounds for bin in 1:(length(bins)-1)
        ∑δL += hist_δ[bin]
        ∑δ²L += hist_δ²[bin]
        ∑𝑤L += hist_𝑤[bin]
        ∑δR -= hist_δ[bin]
        ∑δ²R -= hist_δ²[bin]
        ∑𝑤R -= hist_𝑤[bin]

        gainL, gainR = get_gain(params.loss, ∑δL, ∑δ²L, ∑𝑤L, params.λ), get_gain(params.loss, ∑δR, ∑δ²R, ∑𝑤R, params.λ)
        gain = gainL + gainR

        if gain > info.gain && ∑𝑤L[1] >= params.min_weight && ∑𝑤R[1] >= params.min_weight
            info.gain = gain
            info.gainL = gainL
            info.gainR = gainR
            info.∑δL = ∑δL
            info.∑δ²L = ∑δ²L
            info.∑𝑤L = ∑𝑤L
            info.∑δR = ∑δR
            info.∑δ²R = ∑δ²R
            info.∑𝑤R = ∑𝑤R
            info.cond = edges[bin]
            info.𝑖 = bin
        end
    end
    return
end

function update_hist!(hist_δ::Vector{Vector{SVector{L,T}}}, hist_δ²::Vector{Vector{SVector{L,T}}}, hist_𝑤::Vector{Vector{SVector{1,T}}}, X_bin, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}, set::Vector{S}, j::Int) where {L,T,S}
    # build histogram
    hist_δ[j] .*= 0.0
    hist_δ²[j] .*= 0.0
    hist_𝑤[j] .*= 0.0
    @inbounds @simd for i in set
        hist_δ[j][X_bin[i,j]] += δ[i]
        hist_δ²[j][X_bin[i,j]] += δ²[i]
        hist_𝑤[j][X_bin[i,j]] += 𝑤[i]
    end
end

function find_split!(hist_δ::Vector{SVector{L,T}}, hist_δ²::Vector{SVector{L,T}}, hist_𝑤::Vector{SVector{1,T}}, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, params::EvoTreeRegressor, info::SplitInfo{L,T,S}, edges::Vector{T}, j::Int) where {L,T,S}

    # initialize tracking
    ∑δL = ∑δ * 0
    ∑δ²L = ∑δ² * 0
    ∑𝑤L = ∑𝑤 * 0
    ∑δR = ∑δ
    ∑δ²R = ∑δ²
    ∑𝑤R = ∑𝑤

    @inbounds for bin in 1:(length(hist_δ)-1)
        ∑δL += hist_δ[bin]
        ∑δ²L += hist_δ²[bin]
        ∑𝑤L += hist_𝑤[bin]
        ∑δR -= hist_δ[bin]
        ∑δ²R -= hist_δ²[bin]
        ∑𝑤R -= hist_𝑤[bin]

        gainL, gainR = get_gain(params.loss, ∑δL, ∑δ²L, ∑𝑤L, params.λ), get_gain(params.loss, ∑δR, ∑δ²R, ∑𝑤R, params.λ)
        gain = gainL + gainR

        if gain > info.gain && ∑𝑤L[1] >= params.min_weight && ∑𝑤R[1] >= params.min_weight
            info.gain = gain
            info.gainL = gainL
            info.gainR = gainR
            info.∑δL = ∑δL
            info.∑δ²L = ∑δ²L
            info.∑𝑤L = ∑𝑤L
            info.∑δR = ∑δR
            info.∑δ²R = ∑δ²R
            info.∑𝑤R = ∑𝑤R
            info.cond = edges[bin]
            info.𝑖 = bin
        end
    end
end

function find_split_narrow!(hist_δ::Vector{SVector{L,T}}, hist_δ²::Vector{SVector{L,T}}, hist_𝑤::Vector{SVector{1,T}}, bins::Vector{BitSet}, X_bin, δ::Vector{SVector{L,T}}, δ²::Vector{SVector{L,T}}, 𝑤::Vector{SVector{1,T}}, ∑δ::SVector{L,T}, ∑δ²::SVector{L,T}, ∑𝑤::SVector{1,T}, params::EvoTreeRegressor, info::SplitInfo{L,T,S}, edges::Vector{T}, set::Vector{S}) where {L,T,S}

    # initialize histogram
    hist_δ .*= 0.0
    hist_δ² .*= 0.0
    hist_𝑤 .*= 0.0

    # build histogram
    @inbounds for i in set
        hist_δ[X_bin[i]] += δ[i]
        hist_δ²[X_bin[i]] += δ²[i]
        hist_𝑤[X_bin[i]] += 𝑤[i]
    end
    return
end


text:
@test callback(fsm, "onbeforeevent", "hop", "first", "second") == 1
@test before_any_event(fsm, "hop", "first", "second") == 1
@test before_event(fsm, "hop", "first", "second") == 1

@test callback(fsm, "onbeforejump", "jump", "third", "fourth") == 3
@test before_any_event(fsm, "jump", "third", "fourth") == 1
@test before_this_event(fsm, "jump", "third", "fourth") == 3
@test before_event(fsm, "jump", "third", "fourth") == 1

@test callback(fsm, "onafterevent", "skip", "second", "third") == 2
@test after_any_event(fsm, "skip", "second", "third") == 2
@test after_this_event(fsm, "skip", "second", "third") == 5
@test after_event(fsm, "skip", "second", "third") == 2

@test change_state(fsm, "hop", "first", "second") == 8
@test enter_any_state(fsm, "hop", "first", "second") == 13
@test enter_this_state(fsm, "hop", "first", "second") == 21
@test leave_any_state(fsm, "hop", "first", "second") == 34
@test leave_this_state(fsm, "hop", "first", "second") == 55

@test after_any_event(mfsm, "hop", "first", "second") == 1
@test enter_any_state(mfsm, "hop", "first", "second") == -1


text:
const STATIC_URL = "https://api-fxpractice.oanda.com/v3/accounts/" * ACCOUNT * "/pricing/?instruments="
const STREAM_URL = "https://stream-fxpractice.oanda.com/v3/accounts/" * ACCOUNT * "/pricing/stream?instruments="
const URL_DELIMITER = "%2C"

function price(instrument::String)
    Requests.get(STATIC_URL * instrument; headers = HEADERS)
end

function price(instruments::Array{String,1})
    price(join(instruments, URL_DELIMITER))
end

function price_stream(instrument::String)
    Requests.get_streaming(STREAM_URL * instrument; headers = HEADERS)
end

function price_stream(instruments::Array{String,1})
    price_stream(join(instruments, URL_DELIMITER))
end

function Base.next(stream::Requests.ResponseStream)
    while nb_available(stream.buffer) == 0
        sleep(.001) #TODO do something useful like yield() instead of just sleep
    end
    tick = JSON.parse(readline(stream.buffer))
    get(tick, "type", "") == "PRICE" ? Tick(tick) : next(stream)
end

#function flatten(tick::Dict{String,Any}) #NOTE ["asks"][1] and ["bids"][1] assume one liquidity provider per practice API, TODO test on live API
#    Dict(
#        "arrival" => Dates.format(now(Dates.UTC), "yyyy-mm-ddTHH:MM:SS.sss"),
#        "instrument" => tick["instrument"],
#        "time" => tick["time"][1:end-7],
#        "liquidityAsk" => tick["asks"][1]["liquidity"],
#        "liquidityBid" => tick["bids"][1]["liquidity"],
#        "closeoutAsk" => parse(Float64,tick["closeoutAsk"]),
#        "closeoutBid" => parse(Float64,tick["closeoutBid"]),
#        "bid" => parse(Float64,tick["bids"][1]["price"]),
#        "ask" => parse(Float64,tick["asks"][1]["price"]),
#        "tradeable" => tick["tradeable"]
#    )
#end

max_stars_repo_name: cadwallader/Oanda

max_issues_repo_name: cadwallader/Oanda

max_forks_repo_name: cadwallader/Oanda


text:
import Base.Broadcast: broadcasted

# ZonedDateTime arithmetic
Base.:(+)(x::ZonedDateTime) = x
Base.:(-)(x::ZonedDateTime, y::ZonedDateTime) = x.utc_datetime - y.utc_datetime

function Base.:(+)(zdt::ZonedDateTime, p::DatePeriod)
    return ZonedDateTime(DateTime(zdt, Local) + p, timezone(zdt))
end
function Base.:(+)(zdt::ZonedDateTime, p::TimePeriod)
    return ZonedDateTime(DateTime(zdt, UTC) + p, timezone(zdt); from_utc=true)
end
function Base.:(-)(zdt::ZonedDateTime, p::DatePeriod)
    return ZonedDateTime(DateTime(zdt, Local) - p, timezone(zdt))
end
function Base.:(-)(zdt::ZonedDateTime, p::TimePeriod)
    return ZonedDateTime(DateTime(zdt, UTC) - p, timezone(zdt); from_utc=true)
end

function broadcasted(::typeof(+), r::StepRange{ZonedDateTime}, p::DatePeriod)
    start, step, stop = first(r), Base.step(r), last(r)

    # Since the local time + period can result in an invalid local datetime when working with
    # VariableTimeZones we will use `first_valid` and `last_valid` which avoids issues with
    # non-existent and ambiguous dates.

    tz = timezone(start)
    if isa(tz, VariableTimeZone)
        start = first_valid(DateTime(start, Local) + p, tz, step)
    else
        start = start + p
    end

    tz = timezone(stop)
    if isa(tz, VariableTimeZone)
        stop = last_valid(DateTime(stop, Local) + p, tz, step)
    else
        stop = stop + p
    end

    return StepRange(start, step, stop)
end

broadcasted(::typeof(-), r::StepRange{ZonedDateTime}, p::DatePeriod) = broadcast(+, r, -p)


text:
using Test, SineWaves, FFTW

@testset "SineWave constructor tests" begin
    frequency = 440.0
    samplerate = 48000.0
    delta = 2*pi*frequency/samplerate
    sinewave = SineWave(frequency, samplerate)

    @test sinewave.current ≈ 0.0
    @test sinewave.previous ≈ -sin(delta)
    @test sinewave.cosine ≈ cos(delta)

    too_hi_frequency = 40000.0

    @test_throws ErrorException SineWave(too_hi_frequency, samplerate)
end

@testset "fill! function tests" begin
    frequency = 440.0
    samplerate = 48000.0
    delta = 2*pi*frequency/samplerate
    N = 512

    sinewave = SineWave(frequency, samplerate)
    buffer = zeros(Float64, N)
    fill!(buffer, sinewave)

    @test buffer ≈ sin.(range(delta; step=delta, length = N))
end

@testset "spectrum function tests" begin
    frequency = 440.0
    samplerate = 48000.0
    delta = 2*pi*frequency/samplerate
    N = 512

    sinewave = SineWave(frequency, samplerate)
    buffer = zeros(Float64, N)
    fill!(buffer, sinewave)

    buffer_ft = FFTW.rfft(buffer)
    expected = abs.(buffer_ft)

    @test spectrum(buffer) ≈ expected
end


text:
using Distances
export delete_insert!, vectorisable_metrics

function get_vectorisable_metrics(model::MultigraphSNF)
    return [Cityblock]
end

function is_vectorisable(model::MultigraphSNF)
    return typeof(model.d) ∈ get_vectorisable_metrics(model)
end

# Accept reject for vectorised scheme
function accept_reject!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    mcmc::SnfMcmcInsertDelete,
    model::MultigraphSNF,
    x_mode::Vector{Int} # Must also pass vectorised mode
    )
    M_curr = sum(x_curr)    # Total number of edges in current state
    M_prop = M_curr         # Total number of edges in proposal (to vary)
    δ = rand(1:mcmc.ν)
    d = rand(0:min(δ,M_curr)) # Number of edges to delete
    log_ratio = 0.0

    # Delete edges
    # This uses random sampling from urn. That is, imagine an urn containing all edges, with multiplicity.
    # We sample edges to delete by drawing from this urn.
    # See sample_urn_without_replacement() in InteractionNetworkModels\src\Graph Models\utils.jl
    for i in 1:d
        counter = 0
        z = rand(1:M_prop)
        for (edge, edge_weight) in enumerate(x_prop)
            counter += edge_weight
            if counter ≥ z
                # We are going to take away a single edge
                log_ratio -= log(x_prop[edge])
                x_prop[edge] -= 1
                M_prop -= 1
                break
            end
        end
    end
    # Add edges
    # println("Additons:", δ-d)
    N = length(x_curr)
    for i in 1:(δ-d)
        edge = rand(1:N)
        x_prop[edge] += 1
        log_ratio += log(x_prop[edge])
        M_prop += 1
    end
    # Add log(M_curr!/M_prop!)
    log_size_factorial = 0.0
    if M_curr > M_prop
        for i in (M_prop+1):M_curr
            log_size_factorial += log(i)
        end
    elseif M_prop > M_curr
        for i in (M_curr+1):M_prop
            log_size_factorial -= log(i)
        end
    end
    # Now add term for insertions
    log_insertion = (M_prop-M_curr)*log(N)     # log((1/V)^(d-a)) = (M_prop-M_curr)*log(V) (n.b. M_curr = M_prop - d + a)

    log_ratio += log_size_factorial + log_insertion

    log_ratio += log(min(M_curr, δ)+1) - log(min(M_prop, δ)+1)
    # With log_α denoting the log acceptance probability and
    #       f(x) = exp(- γ d(x, xᵐ))
    # denoting the un-normalised probability of x, we have
    # log_α = f(x_prop) - f(x_curr) + log_ratio
    # thus the following...
    γ, d = (model.γ, model.d)
    # println(x_curr==x_prop)
    log_model = γ * (d(x_curr, x_mode) - d(x_prop, x_mode))
    log_α =  log_model + log_ratio
    # @show log_model, log_ratio, M_curr, M_prop, log_size_factorial, log_insertion
    if log(rand()) < log_α
        copy!(x_curr, x_prop)
        return true
    else
        copy!(x_prop, x_curr)
        return false
    end
end

function accept_reject_test!(
    x_curr::Vector{Int},
    x_prop::Vector{Int},
    mcmc::SnfMcmcInsertDelete,
    model::MultigraphSNF,
    x_mode::Vector{Int} # Must also pass vectorised mode
    )
    M_curr = sum(x_curr)    # Total number of edges in current state
    δ = rand(1:mcmc.ν)
    d = rand(0:min(δ,M_curr)) # Number of edges to delete
    log_ratio = 0.0

    nz_edges = Set(findall(x->x>0, x_prop))
    # Delete edges
    for i in 1:d
        # @show edge_ind
        edge = rand(edge_ind)
        log_ratio += log(length(edge_ind))
        x_prop[edge] -= 1
        if x_prop[edge] == 0
            delete!(nz_edges, edge)
        end
    end
    # Add edges
    # println("Additons:", δ-d)
    N = length(x_curr)
    for i in 1:(δ-d)
        # @show edge_ind
        edge = rand(1:N)
        if x_prop[edge] == 0
            push!(nz_edges, edge)
        end
        x_prop[edge] += 1
        log_ratio -= log(length(edge_ind))
    end
    # Now add term for insertions

    log_ratio += (δ - 2d)*log(N)     # log((1/N)^(d-a)) = (M_prop-M_curr)*log(N) (n.b. M_curr = M_prop - d + a)
    M_prop = sum(x_prop)
    log_ratio += log(min(M_curr, δ)+1) - log(min(M_prop, δ)+1)
    # With log_α denoting the log acceptance probability and
    #       f(x) = exp(- γ d(x, xᵐ))
    # denoting the un-normalised probability of x, we have
    # log_α = f(x_prop) - f(x_curr) + log_ratio
    # thus the following...
    γ, d = (model.γ, model.d)
    # println(x_curr==x_prop)
    log_model = γ * (d(x_curr, x_mode) - d(x_prop, x_mode))
    log_α =  log_model + log_ratio
    # @show log_model, log_ratio, M_curr, M_prop
    if log(rand()) < log_α
        copy!(x_curr, x_prop)
        return true
    else
        copy!(x_prop, x_curr)
        return false
    end
end

function draw_sample!(
    out::Vector{Vector{Int}},
    mcmc::SnfMcmcInsertDelete,
    model::MultigraphSNF;
    burn_in::Int=1000, lag::Int=10,
    init::Vector{Int}=adj_mat_to_vec(model.mode, directed=model.directed, self_loops=model.self_loops)
    )

    x_curr = copy(init)
    x_prop = copy(x_curr)
    x_mode = adj_mat_to_vec(
        model.mode,
        directed=model.directed,
        self_loops=model.self_loops
    )
    sample_count = 1
    accept_count = 0
    iter_count = 0
    while sample_count ≤ length(out)
        iter_count += 1
        # Store value if beyond burn-in and coherent with lags
        if (iter_count > burn_in) & (((iter_count-1) % lag)==0)
            @inbounds copy!(out[sample_count], x_curr)
            sample_count += 1
        end
        was_accepted = accept_reject!(
            x_curr, x_prop,
            mcmc, model,
            x_mode
        )
        accept_count += was_accepted

    end
    return accept_count / iter_count
end

function draw_sample(
    mcmc::T,
    model::MultigraphSNF;
    desired_samples::Int=1000,
    burn_in::Int=1000, lag::Int=10,
    init::Union{Vector{Int},Matrix{Int}}=model.mode
    ) where {T<:SnfMcmcSampler}

    # We first determine whether for this metric we have a vectorised implementation
    vectorised = is_vectorisable(model)

    if vectorised
        if typeof(init)==Matrix{Int}
            init_vec = adj_mat_to_vec(
                init,
                directed=model.directed,
                self_loops=model.self_loops
            )
        else
            init_vec = init
        end
        out = [Int[] for i in 1:desired_samples]
        accept_prob = draw_sample!(
            out, mcmc, model,
            burn_in=burn_in, lag=lag,
            init=init_vec
        )
        # Now convert back to matrices
        directed, self_loops = (model.directed, model.self_loops)
        out = [vec_to_adj_mat(x, directed=directed, self_loops=self_loops) for x in out]
        return out, accept_prob

    else
        # To implement
        error("Non-vectorised scheme not yet implemented")
    end

end

function (mcmc::T where {T<:SnfMcmcSampler})(
    model::SNF;
    desired_samples::Int=mcmc.desired_samples,
    burn_in::Int=mcmc.burn_in,
    lag::Int=mcmc.lag,
    init::Union{Vector{Int},Matrix{Int}}=model.mode
    )

    out, a = draw_sample(
        mcmc,
        model,
        desired_samples=desired_samples,
        burn_in=burn_in,
        lag=lag,
        init=init
    )

    p_measures = Dict(
            "Acceptance probability" => a
    )

    output = SnfMcmcOutput(
            model,
            out,
            p_measures
    )
    return output

end

max_stars_repo_path: src/Graph Models/SNF/SNF_multigraph.jl

max_issues_repo_path: src/Graph Models/SNF/SNF_multigraph.jl

max_forks_repo_path: src/Graph Models/SNF/SNF_multigraph.jl


text:
using Revise
using BridgeLandmarks
const BL = BridgeLandmarks
using Distributions
using LinearAlgebra
using Random
using FileIO

Random.seed!(9)

workdir = @__DIR__
cd(workdir)
outdir = workdir

#-------- read data ----------------------------------------------------------
dat = load("data_cardiac.jld2")
xobs0 = dat["xobs0"]
xobsT = dat["xobsT"]

# set pars
p_ms = Pars_ms(δpos=0.00002, δmom=0.1, cinit=0.2, γinit=1.0, σobs = 0.01, ρinit=0.3)
p_ahs = Pars_ahs(δpos=0.01, δmom=0.1, cinit=0.05, γinit=0.5, db=[2.5,2.5],stdev=0.75, σobs = 0.01, ρinit=0.4)

ITER = 2000 # taken in paper
template_estimation(xobsT; ainit=0.1, xinitq=xobsT[1],
pars = p_ms, outdir=outdir, ITER=200, updatescheme = [:innov, :rmmala_pos, :parameter])  # deliberately initialise badly to show it works

plottemplate_estimation(outdir)


text:
### ConvAsymPadding Layer ###

struct ConvAsymPadding{T, N}
    c::T
    pad::NTuple{N,Int}
end
function ConvAsymPadding(k::NTuple{N,Integer}, ch::Pair{<:Integer,<:Integer}, σ = identity;
    init = Flux.glorot_uniform,  stride = 1, pad = 0, dilation = 1, bias=true) where N
    length(pad) < 2N || all(i -> length(unique(pad[i:i+1])) == 1, 1:2:2N) == 1&& return Conv(k, ch , σ, init=init, stride=stride, pad=pad, dilation=dilation, bias=bias)

    pad_manual = Tuple(map(i -> abs(pad[i] - pad[i+1]), 1:2:2N))
    pad_auto = Tuple(map(i -> minimum(pad[i:i+1]), 1:2:2N))
    return ConvAsymPadding(Conv(k, ch, σ, init=init, stride=stride, pad=pad_auto, dilation=dilation), pad_manual)
end
function (c::ConvAsymPadding)(x::AbstractArray)
    # Maybe there are faster ways to do this as well...
    padding = similar(x, c.pad..., size(x)[end-1:end]...)
    fill!(padding, 0)
    c.c(cat(x, padding, dims=1:length(c.pad)))
end

Flux.@functor ConvAsymPadding

### Cropping Layer ###

struct Crop{N}
    crop::NTuple{N,Int}
end

function (c::Crop{4})(x::AbstractArray)
  return x[(1+c.crop[1]):(end-c.crop[2]), (1+c.crop[3]):(end-c.crop[4]),:,:]
end

function (c::Crop{6})(x::AbstractArray)
  return x[(1+c.crop[1]):(end-c.crop[2]), (1+c.crop[3]):(end-c.crop[4]), (1+c.crop[5]):(end-c.crop[6]),:,:]
end

Flux.@functor Crop

### BatchNormWrap ###

expand_dims(x,n::Int) = reshape(x,ones(Int64,n)...,size(x)...)

function _random_normal(shape...)
  return Float32.(rand(Normal(0.f0,0.02f0),shape...)) |> device
end

function BatchNormWrap(out_ch)
    Chain(x->expand_dims(x,3),
	  BatchNorm(out_ch),
	  x->reshape(x, size(x)[4:end]))
end

### UNetConvBlock ###

UNetConvBlock(in_chs, out_chs; kernel = (3,3,3), pad = (1,1,1), stride=(1,1,1)) =
  Chain(Conv(kernel, in_chs=>out_chs, pad = pad, stride=stride; init=_random_normal),
	BatchNormWrap(out_chs),
	x->leakyrelu.(x,0.02f0))


text:
function readData(fname::AbstractString, shape, minFeatureInd::Integer)
	features = zeros(shape)
	labels = Float64[]
	fi = open(fname, "r")
	cnt = 1
	for line in eachline(fi)
		line = split(line, " ")
		push!(labels, parse(Float64,line[1]))
		line = line[2:end]
		for itm in line
			itm = split(itm, ":")
			features[cnt, parse(Int,itm[1]) + 1-minFeatureInd] = parse(Float64,chomp(itm[2]))
		end
		cnt += 1
	end
	close(fi)
	(features, labels)
end

function readDataSparse(fname::AbstractString, shape, minFeatureInd::Integer)
	A = Int64[]
	B = Int64[]
	C = Float64[]
	labels = Float64[]
	# featuresVec =  SparseMatrixCSC{Float64,Int64}[]
	# featuresVec2 =  SparseMatrixCSC{Float64,Int64}[]
	fi = open(fname, "r")
	cnt = 1

	for line in eachline(fi)
		# A1 = Int64[]
		# B1= Int64[]
		# C1 = Float64[]
		line = split(line, " ")
		push!(labels, float64(line[1]))
		line = line[2:end]
		for itm in line
			itm = split(itm, ":")
			push!(A, cnt)
			push!(B, int(itm[1]) + 1-minFeatureInd)
			push!(C, float64(chomp(itm[2])))

			# push!(A1, int(itm[1]) + 1-minFeatureInd)
			# push!(B1, 1)
			# push!(C1, float64(chomp(itm[2])))

		end

		# push!(featuresVec, sparse(A1, B1, C1,shape[2],1))
		# push!(featuresVec2, sparse( B1,A1, C1,1,shape[2]))
		cnt += 1
	end
	close(fi)

	features = sparse(A,B,C)

	# (features, featuresVec, featuresVec2,  labels)
	(features, labels)
end

function writeBin(fname::AbstractString, features, labels)
	fi = open(fname, "w")
	write(fi, features)
	write(fi, labels)
	close(fi)
end

function readBin(fname::AbstractString,numDatapoints::Integer,numVars::Integer)
	fi = open(fname, "r")
	features = read(fi, Float64, (numDatapoints,numVars))
	labels = read(fi, Float64, (numDatapoints,))
	close(fi)
	(features,labels)
end

function getStats(fname::AbstractString)
	fi = open(fname, "r")
	n = 0
	numTotal=0
	minfeatureInd=typemax(Int)
	maxfeatureInd=typemin(Int)
	minfeature=typemax(Int)
	maxfeature=typemin(Int)
	mySetOfClasses = Set()
	for line in eachline(fi)
		n += 1
		line = split(line, " ")
		push!(mySetOfClasses, line[1])
		line = line[2:end]
		for itm in line
			itm = split(itm, ":")
			minfeatureInd = min(minfeatureInd,parse(Int64,itm[1]))
			maxfeatureInd = max(maxfeatureInd,parse(Int64,itm[1]))
			minfeature = min(minfeature,parse(Float64,chomp(itm[2])))
			maxfeature = max(maxfeature,parse(Float64,chomp(itm[2])))
			numTotal+=1
		end
	end

	if numTotal!=n*(maxfeatureInd-minfeatureInd+1)
		minfeature = min(minfeature,0)
		maxfeature = max(maxfeature,0)
	end

	close(fi)

	(n,minfeatureInd,maxfeatureInd-minfeatureInd+1,numTotal,minfeature,maxfeature, length(mySetOfClasses))
end

function getStatsHIGGSSUSY(fname::AbstractString)
	fi = open(fname, "r")
	n = 0
	numTotal=0
	minfeature=typemax(Int)
	maxfeature=typemin(Int)
	mySetOfClasses = Set()
	for line in eachline(fi)
		n += 1
		line = split(line, ",")
		push!(mySetOfClasses, line[1])
		line = line[2:end]
		for itm in line
			minfeature = min(minfeature,float64(chomp(itm)))
			maxfeature = max(maxfeature,float64(chomp(itm)))
			numTotal+=1
		end
	end
	close(fi)
	(n,numTotal,minfeature,maxfeature, length(mySetOfClasses))
end

function readDataHIGGSSUSY(fname::AbstractString, shape)
	features = zeros(shape)
	labels = Float64[]
	fi = open(fname, "r")
	cnt = 1
	for line in eachline(fi)
		line = split(line, ",")
		push!(labels, float64(line[1]))
		line = line[2:end]
		for i in 1:length(line)
			features[cnt, i] = float64(chomp(line[i]))
		end
		cnt += 1
	end
	close(fi)
	(features, labels)
end

function getStatsIRIS(fname::AbstractString)
	fi = open(fname, "r")
	n = 0
	numTotal=0
	minfeature=typemax(Int)
	maxfeature=typemin(Int)
	mySetOfClasses = Set()
	for line in eachline(fi)
		n += 1
		line = split(line, ",")
		push!(mySetOfClasses, chomp(line[end]))
		line = line[1:end-1]
		for itm in line
			minfeature = min(minfeature,float64(chomp(itm)))
			maxfeature = max(maxfeature,float64(chomp(itm)))
			numTotal+=1
		end
	end
	close(fi)
	println(mySetOfClasses)
	(n,numTotal,minfeature,maxfeature, length(mySetOfClasses))
end

function readDataIRIS(fname::AbstractString, shape)
	features = zeros(shape)
	labels = String[]
	fi = open(fname, "r")
	cnt = 1
	for line in eachline(fi)
		line = split(line, ",")
		push!(labels, chomp(line[end]))
		line = line[1:end-1]
		for i in 1:length(line)
			features[cnt, i] = float64(chomp(line[i]))
		end
		cnt += 1
	end
	close(fi)
	(features, labels)
end


text:
"""
    AbstractLink

Abstract type defining maps from R^n -> X.
They can be applied by calling `link(x)`.

A series of definitions are given in http://web.pdx.edu/~newsomj/mvclass/ho_link.pdf
"""
abstract type AbstractLink end

struct ChainLink{Tls} <: AbstractLink
    links::Tls
end

(l::ChainLink)(x) = foldl((x, l) -> l(x), l.ls; init=x)

"""
    Link(f)

General construction for a link with a function `f`.
"""
struct Link{F} <: AbstractLink
    f::F
end

(l::Link)(x) = l.f(x)

"""
    LogLink()

`log` link, f:ℝ⁺->ℝ . Its inverse is the [`ExpLink`](@ref).
"""
struct LogLink <: AbstractLink end

(::LogLink)(x) = log(x)

Base.inv(::LogLink) = ExpLink()

"""
    ExpLink()

`exp` link, f:ℝ->ℝ⁺. Its inverse is the [`LogLink`](@ref).
"""
struct ExpLink <: AbstractLink end

(::ExpLink)(x) = exp(x)

Base.inv(::ExpLink) = LogLink()

"""
    InvLink()

`inv` link, f:ℝ/{0}->ℝ/{0}. It is its own inverse.
"""
struct InvLink <: AbstractLink end

(::InvLink)(x) = inv(x)

Base.inv(::InvLink) = InvLink()

"""
    SqrtLink()

`sqrt` link, f:ℝ⁺∪{0}->ℝ⁺∪{0}. Its inverse is the [`SquareLink`](@ref).
"""
struct SqrtLink <: AbstractLink end

(::SqrtLink)(x) = sqrt(x)

Base.inv(::SqrtLink) = SquareLink()

"""
    SquareLink()

`^2` link, f:ℝ->ℝ⁺∪{0}. Its inverse is the [`SqrtLink`](@ref).
"""
struct SquareLink <: AbstractLink end

(::SquareLink)(x) = x^2

Base.inv(::SquareLink) = SqrtLink()

"""
    LogitLink()

`log(x/(1-x))` link, f:[0,1]->ℝ. Its inverse is the [`LogisticLink`](@ref).
"""
struct LogitLink <: AbstractLink end

(::LogitLink)(x) = logit(x)

Base.inv(::LogitLink) = LogisticLink()

"""
    LogisticLink()

`exp(x)/(1+exp(-x))` link. f:ℝ->[0,1]. Its inverse is the [`Logit`](@ref).
"""
struct LogisticLink <: AbstractLink end

(::LogisticLink)(x) = logistic(x)

Base.inv(::LogisticLink) = LogitLink()

"""
    ProbitLink()

`ϕ⁻¹(y)` link, where `ϕ⁻¹` is the `invcdf` of a `Normal` distribution, f:[0,1]->ℝ.
Its inverse is the [`NormalCDFLink`](@ref).
"""
struct ProbitLink <: AbstractLink end

(::ProbitLink)(x) = norminvcdf(x)

Base.inv(::ProbitLink) = NormalCDFLink()

"""
    NormalCDFLink()

`ϕ(y)` link, where `ϕ` is the `cdf` of a `Normal` distribution, f:ℝ->[0,1].
Its inverse is the [`ProbitLink`](@ref).
"""
struct NormalCDFLink <: AbstractLink end

(::NormalCDFLink)(x) = normcdf(x)

Base.inv(::NormalCDFLink) = ProbitLink()

(::ChainLink{<:Tuple{LogLink,NormalCDFLink}})(x) = normlogcdf(x) # Specialisation for log + normal cdf

"""
    SoftMaxLink()

`softmax` link, i.e `f(x)ᵢ = exp(xᵢ)/∑ⱼexp(xⱼ)`.
f:ℝⁿ->Sⁿ⁻¹, where Sⁿ⁻¹ is an [(n-1)-simplex](https://en.wikipedia.org/wiki/Simplex)
It has no defined inverse
"""
struct SoftMaxLink <: AbstractLink end

(::SoftMaxLink)(x::AbstractVector{<:Real}) = softmax(x)


text:
# GFG3015 digital multimeter driver

mutable struct GFG3015
    handle
	initialized
    instr_dict
	wfm_dict
	rev_wfm_dict
	unit_dict
	fct
end
GFG3015(handle) = GFG3015(
                handle,
				false,
                Dict(
                    "C1"=>"1"
                    ),
                Dict(
					"sinusoid"=>"1",
					"triangle"=>"2",
					"square"=>"3",
                    ),
                Dict(
                    "1"=>"sinusoid",
					"2"=>"triangle",
					"3"=>"square",
                    ),
                Dict(
					"Vpp"=>"1",
					"Vrms"=>"2",
					"dBm"=>"3",
					"1"=>"Vpp",
					"2"=>"Vrms",
					"3"=>"dBm",
                    ),
					"sinusoid",
                    )

# MISC
function connect!(obj::GFG3015, handle)
	obj.handle = handle
    obj.connected = 1
	return nothing
end
function disconnect!(obj::GFG3015)
	obj.handle = 0
    obj.connected = 0
	return nothing
end
reset_instr(obj::GFG3015) = write(obj.handle, "*RST")
clear_instr(obj::GFG3015) = write(obj.handle, "*CLS")
get_idn(obj::GFG3015) = query(obj.handle, "*IDN?")

## WFM
function set_wfm(obj::GFG3015, ch, fct="sinusoid")
	# set_wfm(fgen, "C1", "sinusoid")
	ch = obj.instr_dict[ch]
	fct = obj.wfm_dict[fct]
	cmd = ":FUNCtion:WAVeform $fct"
	write(obj.handle, cmd)
end
function get_wfm(obj::GFG3015, ch)
	# get_wfm(fgen, "C1")
	ch = obj.instr_dict[ch]
	cmd = ":FUNCtion:WAVeform?"
	value = strip(query(obj.handle, cmd))
	value = obj.rev_wfm_dict[value]
end

# UNIT
function set_amplit_unit(obj::GFG3015, ch, unit="Vpp")
	# set_amplit_unit(fgen, "C1", "Vpp")
	ch = obj.instr_dict[ch]
	unit = obj.unit_dict[unit]
	cmd = ":AMPLitude:UNIT $unit"
	write(obj.handle, cmd)
end
function get_amplit_unit(obj::GFG3015, ch)
	# get_amplit_unit(fgen, "C1")
	ch = obj.instr_dict[ch]
	cmd = ":AMPLitude:UNIT?"
	value = strip(query(obj.handle, cmd))
	value = obj.unit_dict[value]
end

# AMPLITUDE
function set_amplit(obj::GFG3015, ch, value)
	# set_amplit(fgen, "C1", 2.2)
	ch = obj.instr_dict[ch]
	cmd = ":AMPLitude:VOLTage $value"
	write(obj.handle, cmd)
end
function get_amplit(obj::GFG3015, ch)
	# get_amplit(fgen, "C1")
	ch = obj.instr_dict[ch]
	cmd = ":AMPLitude:VOLTage?"
	value = parse(Float64, strip(query(obj.handle, cmd)))
end

# OFFSET
function set_offs(obj::GFG3015, ch, value)
	# set_offs(fgen, "C1", 1.1)
	ch = obj.instr_dict[ch]
	cmd = ":OFFSet $value"
	write(obj.handle, cmd)
end
function get_offs(obj::GFG3015, ch)
	# get_offs(fgen, "C1")
	ch = obj.instr_dict[ch]
	cmd = ":OFFSet?"
	value = parse(Float64, strip(query(obj.handle, cmd)))
end

# FREQUENCY
function set_freq(obj::GFG3015, ch, value)
	# set_freq(fgen, "C1", 1333)
	ch = obj.instr_dict[ch]
	cmd = ":FREQuency $value"
	write(obj.handle, cmd)
end
function get_freq(obj::GFG3015, ch)
	# get_freq(fgen, "C1")
	ch = obj.instr_dict[ch]
	cmd = ":FREQuency?"
	value = parse(Float64, strip(query(obj.handle, cmd)))
end

# DUTY CYCLE
function set_duty(obj::GFG3015, ch, value)
	# set_duty(fgen, "C1", 40)
	ch = obj.instr_dict[ch]
	cmd = ":DUTY $value"
	write(obj.handle, cmd)
end
function get_duty(obj::GFG3015, ch)
	# get_duty(fgen, "C1")
	ch = obj.instr_dict[ch]
	cmd = ":DUTY?"
	value = parse(Float64, strip(query(obj.handle, cmd)))
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
@testset "Spectral downsampling" begin
    # Test Fourier downsampling and signal truncation which is implemented to
    # be compatible with the Siemens way of doing this.
    t = (0:15)/16
    z = exp.(2*π*1im .* t)
    z = reshape(z, length(z), 1) # Add trivial channels dim

    # Know decimated z is numerically equal to Fourier downsampled z, as it's
    # got no frequency content outside the window.
    @test z[1:2:end] ≈ MagneticResonanceSignals.downsample_and_truncate(t, z, 0, 0, 2)[2]  rtol=0 atol=1e-14

    # Currently can't downsample to an odd number of samples
    @test_throws InexactError MagneticResonanceSignals.downsample_and_truncate(t, z, 1, 0, 2)
    @test_throws InexactError MagneticResonanceSignals.downsample_and_truncate(t, z, 0, 1, 2)
    @test_throws InexactError MagneticResonanceSignals.downsample_and_truncate(t, z, 0, 0, 3)

    # Test timing.
    @test t[3:2:end-2] == MagneticResonanceSignals.downsample_and_truncate(t, z, 2, 2, 2)[1]
    @test t[2:end-1] == MagneticResonanceSignals.downsample_and_truncate(t, z, 1, 1, 1)[1]
end

@testset "coil_combination" begin
    # TODO!!
end

@testset "lcosy" begin
    lcosy = @test_logs (:warn,) (:warn,) mr_load("twix/sub-SiemensBrainPhantom_seq-svslcosy_inc-1.twix")

    signal = simple_averaging(lcosy)
    @test size(signal) == (2048,1)
    @test axisnames(signal) == (:time2, :time1)

    # NB: The value of signal[1,1] here is very sensitive to the detail of how
    # channel combination is done.
    #
    # If this test breaks, check channel combination first.
    #
    # Note that we check this in float32 precision, because part of the signal
    # processing is done there. In particular, reduction order in sum()/mean()
    # is implementation defined, so there can be some 1ULP or so difference in
    # simple_averaging() which are system-dependent and cause the tests to be
    # unreliable.
    @test signal[1,1] ≈ 7.738675f-7 - 3.6804267f-6im

    signal = simple_averaging(lcosy, downsample=2)
    @test size(signal) == (1024,1)
    # See comment above regarding channel combination.
    @test signal[1,1] ≈ 5.918123f-6 - 3.929549f-6im
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
"""
    SourceVertex{T} <: AbstractVertex
    SourceVertex(data::T, name::String)

`AbstractVertex` which returns `data` when invoked with no arguments.
"""
struct SourceVertex{T} <: AbstractVertex
    data::T
    name::String
end
SourceVertex(data) = SourceVertex(data, "SourceVertex")

(v::SourceVertex)() = v.data

NaiveNASlib.name(v::SourceVertex) = v.name
NaiveNASlib.trait(v::SourceVertex) = Immutable()

NaiveNASflux.layertype(v::SourceVertex) = guess_layertype(ndims(v.data))

NaiveNASlib.nout(v::SourceVertex) = size(v.data, NaiveNASflux.outdim(layertype(v)))
NaiveNASlib.nin(::SourceVertex) = Integer[]

# TODO: Move to NaiveNASlib
NaiveNASlib.nin(v::OutputsVertex) = nin(base(v))
NaiveNASlib.nout(v::OutputsVertex) = nout(base(v))

NaiveNASlib.inputs(::SourceVertex) = AbstractVertex[]

function sourcevertex_with_outputs(data, name)
    sv = SourceVertex(data, name)
    ov = OutputsVertex(sv)
    NaiveNASlib.init!(ov, sv)
    return ov
end


text:
using GridArrays, Test
using GridArrays: PeriodicCartesianIndices

@testset "PeriodicCartesianIndices" begin
    modcart = PeriodicCartesianIndices((30,40,50),10:20,10:30,10:40);
    cart = modcart.iter
    t1 = @timed(for i in cart
        i ∈ modcart
    end)[3]
    t2 = @timed(for i in cart
        i ∈ cart
    end)[3]

    t1 = @timed(for i in cart
        i ∈ modcart
    end)[3]
    t2 = @timed(for i in cart
        i ∈ cart
    end)[3]

    @test t1 == t2

    for i in modcart
        @test i ∈ cart
    end

    modcart = PeriodicCartesianIndices((30,40,50),CartesianIndex(1,1,1),CartesianIndex(30,40,50));
    cart = modcart.iter

    @test modcart == cart
end

@testset "PeriodicUnitRange" begin
    modcart = GridArrays.PeriodicIndexing.PeriodicUnitRange(20,-3:3)
    cart = modcart.iter
    t1 = @timed(for i in cart
        i ∈ modcart
    end)[3]
    t2 = @timed(for i in cart
        i ∈ cart
    end)[3]

    t1 = @timed(for i in cart
        i ∈ modcart
    end)[3]
    t2 = @timed(for i in cart
        i ∈ cart
    end)[3]

    @test t1 == t2

    for i in cart
        @test i ∈ modcart
    end

    j = 0
    for i in modcart
        j+= 1
    end
    @test j==length(modcart)

    @test first(modcart) == 17
    @test last(modcart) == 3
end

@testset "boundary" begin
    g = PeriodicEquispacedGrid(11,-2,2)^2
    d = UnitDisk()
    b = GridArrays.boundary(g,d)
    for bi in b
        @test approx_in(bi, UnitCircle(),1e-10)
    end

    g = PeriodicEquispacedGrid(11,-2,2)^2
    d = UnitDisk()
    b = GridArrays.boundary_grid(g,d)
    @test subindices(b) == CartesianIndex{2}[CartesianIndex(6, 4), CartesianIndex(7, 4), CartesianIndex(5, 5),
    CartesianIndex(6, 5), CartesianIndex(7, 5), CartesianIndex(8, 5), CartesianIndex(4, 6),
    CartesianIndex(5, 6), CartesianIndex(8, 6), CartesianIndex(9, 6), CartesianIndex(4, 7),
    CartesianIndex(5, 7), CartesianIndex(8, 7), CartesianIndex(9, 7), CartesianIndex(5, 8),
    CartesianIndex(6, 8), CartesianIndex(7, 8), CartesianIndex(8, 8), CartesianIndex(6, 9),
    CartesianIndex(7, 9)]
end