text:
export DDPGPolicy

mutable struct DDPGPolicy{
    BA<:NeuralNetworkApproximator,
    BC<:NeuralNetworkApproximator,
    TA<:NeuralNetworkApproximator,
    TC<:NeuralNetworkApproximator,
    P,
    R<:AbstractRNG,
} <: AbstractPolicy

    behavior_actor::BA
    behavior_critic::BC
    target_actor::TA
    target_critic::TC
    γ::Float32
    ρ::Float32
    batch_size::Int
    start_steps::Int
    start_policy::P
    update_after::Int
    update_every::Int
    act_limit::Float64
    act_noise::Float64
    step::Int
    rng::R
    # for logging
    actor_loss::Float32
    critic_loss::Float32
end

Flux.functor(x::DDPGPolicy) = (
    ba = x.behavior_actor,
    bc = x.behavior_critic,
    ta = x.target_actor,
    tc = x.target_critic,
),
y -> begin
    x = @set x.behavior_actor = y.ba
    x = @set x.behavior_critic = y.bc
    x = @set x.target_actor = y.ta
    x = @set x.target_critic = y.tc
    x
end

"""
    DDPGPolicy(;kwargs...)

# Keyword arguments

- `behavior_actor`,
- `behavior_critic`,
- `target_actor`,
- `target_critic`,
- `start_policy`,
- `γ = 0.99f0`,
- `ρ = 0.995f0`,
- `batch_size = 32`,
- `start_steps = 10000`,
- `update_after = 1000`,
- `update_every = 50`,
- `act_limit = 1.0`,
- `act_noise = 0.1`,
- `step = 0`,
- `rng = Random.GLOBAL_RNG`,
"""
function DDPGPolicy(;
    behavior_actor,
    behavior_critic,
    target_actor,
    target_critic,
    start_policy,
    γ = 0.99f0,
    ρ = 0.995f0,
    batch_size = 32,
    start_steps = 10000,
    update_after = 1000,
    update_every = 50,
    act_limit = 1.0,
    act_noise = 0.1,
    step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(behavior_actor, target_actor)  # force sync
    copyto!(behavior_critic, target_critic)  # force sync
    DDPGPolicy(
        behavior_actor,
        behavior_critic,
        target_actor,
        target_critic,
        γ,
        ρ,
        batch_size,
        start_steps,
        start_policy,
        update_after,
        update_every,
        act_limit,
        act_noise,
        step,
        rng,
        0.0f0,
        0.0f0,
    )
end

# TODO: handle Training/Testing mode
function (p::DDPGPolicy)(env)
    p.step += 1

    if p.step <= p.start_steps
        p.start_policy(env)
    else
        D = device(p.behavior_actor)
        s = state(env)
        s = Flux.unsqueeze(s, ndims(s) + 1)
        action = p.behavior_actor(send_to_device(D, s)) |> vec |> send_to_host
        clamp(action[] + randn(p.rng) * p.act_noise, -p.act_limit, p.act_limit)
    end
end

function RLBase.update!(p::DDPGPolicy, traj::CircularArraySARTTrajectory)
    length(traj) > p.update_after || return
    p.step % p.update_every == 0 || return
    inds, batch = sample(p.rng, traj, BatchSampler{SARTS}(p.batch_size))
    update!(p, batch)
end

function RLBase.update!(p::DDPGPolicy, batch::NamedTuple{SARTS})
    s, a, r, t, s′ = send_to_device(device(p), batch)

    A = p.behavior_actor
    C = p.behavior_critic
    Aₜ = p.target_actor
    Cₜ = p.target_critic

    γ = p.γ
    ρ = p.ρ

    # !!! we have several assumptions here, need revisit when we have more complex environments
    # state is vector
    # action is scalar
    a′ = Aₜ(s′)
    qₜ = Cₜ(vcat(s′, a′)) |> vec
    y = r .+ γ .* (1 .- t) .* qₜ
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(C)) do
        q = C(vcat(s, a)) |> vec
        loss = mean((y .- q) .^ 2)
        ignore() do
            p.critic_loss = loss
        end
        loss
    end

    update!(C, gs1)

    gs2 = gradient(Flux.params(A)) do
        loss = -mean(C(vcat(s, A(s))))
        ignore() do
            p.actor_loss = loss
        end
        loss
    end

    update!(A, gs2)

    # polyak averaging
    for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
        dest .= ρ .* dest .+ (1 - ρ) .* src
    end
end


text:
import ChainRulesCore: rrule, no_rrule
import ChainRulesCore: rrule_via_ad, RuleConfig, NoForwardsMode, HasReverseMode
import Ghost: make_name, Input, to_expr

###############################################################################
#                              Primitives                                     #
###############################################################################

function function_signatures(fn)
    rrule_methods = methods(fn).ms
    sigs = [remove_first_parameter(rr.sig) for rr in rrule_methods]
    # add keyword version of these functions as well
    kw_sigs = [kwsig for kwsig in map(kwfunc_signature, sigs) if kwsig !== Tuple{}]
    return [sigs; kw_sigs]
end

const CHAINRULES_PRIMITIVES = Ref(FunctionResolver{Bool}())
const NUM_CHAINRULES_METHODS = Ref{Int}(0)

function update_chainrules_primitives!(;force=false)
    num_methods = length(methods(rrule)) + length(methods(no_rrule))
    if force || num_methods != NUM_CHAINRULES_METHODS[]
        sigs_flags = [
            [sig => true for sig in function_signatures(rrule)];
            [sig => false for sig in function_signatures(no_rrule)]  # override rrule(sig...)
            ]
        P = FunctionResolver{Bool}(sigs_flags)
        CHAINRULES_PRIMITIVES[] = P
        NUM_CHAINRULES_METHODS[] = num_methods
    end
end

is_chainrules_primitive(sig) = CHAINRULES_PRIMITIVES[][sig] == true

###############################################################################
#                              RuleConfig                                     #
###############################################################################

struct YotaRuleConfig <: RuleConfig{Union{NoForwardsMode, HasReverseMode}} end

function to_rrule_expr(tape::Tape)
    fn_name = gensym("rrule_$(tape[V(1)].val)")
    header = Expr(:call, fn_name)
    for v in inputs(tape)
        op = tape[v]
        push!(header.args, Expr(:(::), make_name(op), op.typ))
    end
    body = Expr(:block)
    # generate transformed forward pass
    seed_id = tape.meta[:seed].id
    for op in tape.ops[1:seed_id - 1]
        op isa Input && continue
        ex = to_expr(op)
        if ex isa Vector
            push!(body.args, ex...)
        else
            push!(body.args, ex)
        end
    end
    # generate pullback
    pb_name = gensym("pullback_$(tape[V(1)].val)")
    pb_ex = :(function $pb_name(dy) end)
    pb_body = pb_ex.args[2]
    empty!(pb_body.args)  # clean from useless linenumber nodes
    push!(pb_body.args, Expr(:(=), make_name(tape.meta[:seed].id), :dy))
    for op in tape.ops[seed_id + 1:length(tape) - 2]
        op isa Input && continue
        ex = to_expr(op)
        if ex isa Vector
            push!(pb_body.args, ex...)
        else
            push!(pb_body.args, ex)
        end
    end
    push!(body.args, pb_ex)
    # generate return
    result_name = make_name(tape[tape.result].args[1].id)
    push!(body.args, Expr(:tuple, result_name, pb_name))
    fn_ex = Expr(:function, header, body)
    return fn_ex
end

"""
    make_rrule(tape::Tape)
    make_rrule(f, args...)

Generate a function equivalent to (but not extending) ChainRulesCore.rrule(),
i.e. returning the primal value and the pullback.

### Examples:

```
foo(x) = 2x + 1
rr = make_rrule(foo, 2.0)
val, pb = rr(foo, 3.0)
pb(1.0)
```
"""
make_rrule(tape::Tape) = Base.eval(@__MODULE__, to_rrule_expr(tape))
make_rrule(f, args...) = make_rrule(gradtape(f, args...))

const GENERATED_RRULE_CACHE = Dict()

function ChainRulesCore.rrule_via_ad(::YotaRuleConfig, f, args...)
    sig = call_signature(f, args...)
    if haskey(GENERATED_RRULE_CACHE, sig)
        rr = GENERATED_RRULE_CACHE[sig]
        return Base.invokelatest(rr, f, args...)
    else
        rr = make_rrule(f, args...)
        GENERATED_RRULE_CACHE[sig] = rr
        return Base.invokelatest(rr, f, args...)
    end
end

###############################################################################
#                                 Rules                                       #
###############################################################################

function ChainRulesCore.rrule(::typeof(broadcasted), ::typeof(identity), x)
    identity_pullback(dy) = (NoTangent(), NoTangent(), dy)
    return x, identity_pullback
end

# test_rrule(broadcasted, identity, [1.0, 2.0]) -- fails at the moment

function ChainRulesCore.rrule(
    ::typeof(Core._apply_iterate), ::typeof(iterate), f::F, args...) where F
    flat = []
    for a in args
        push!(flat, a...)
    end
    return rrule(f, flat...)
end

function ChainRules.rrule(::typeof(tuple), args...)
    y = tuple(args...)
    return y, dy -> (NoTangent(), dy...)
end

# test_rrule(tuple, 1, 2, 3)

ext: jl


text:
# check the Bernoulli numbers and polynomials
#    this is mostly just a check that code is correct
#    most errors are arising because of errors in the functions being used
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot
using SpecialFunctions
import Base: parse
using Printf

Q = Polylogarithms.Q
c = Polylogarithms.c
# b = Polylogarithms.b
# f = Polylogarithms.f
# g = Polylogarithms.g

# from https://github.com/JuliaLang/julia/issues/18328
function parse(::Type{Complex{T}}, s::AbstractString) where {T<:Real}
    s = replace(s, "*^" => "e") # Mathematica seems to export in forms like "3.061616997868383*^-18"
    real_r = r"([+-]?\d*\.\d*(e[+-]\d+)?)([\s+-]|$)"
    imag_r = r"([+-]?\d*\.\d*(e[+-]\d+)?)\s*\*\s*I"

    z = zero(Complex{T})

    m = match( imag_r, s )
    if m != nothing
        z += im*parse(T, m[1])
    end

    m = match( real_r, s )
    if m != nothing
        z += parse(T, m[1])
    end
    return z
end

# input data from Mathematica and reparse into complex numbers
data1 = CSV.read("../data/Q_test_data_1.csv"; delim=",", type=String)
#    has trouble reading in numbers like "2." so read all into strings, and parse
data1[!,:n] = Int.(parse.(Float64, data1[!,:n] ))
data1[!,:s] = parse.(Complex{Float64}, data1[!,:s] )
data1[!,:r] = parse.(Float64, data1[!,:r] )
data1[!,:z] = parse.(Complex{Float64}, data1[!,:z] )
data1[!,:Q] = parse.(Complex{Float64}, data1[!,:Q] )
n = data1[!,:n]
m = size(data1,1)
s = data1[!,:s]
r = data1[!,:r]
μ = log.( data1[!,:z] )
ℒ = log.(complex.(-μ))  # '\u2112'

# figure(1)
# semilogx( data1[!,:r], real( data1[!,:Q] ), "b" )
# semilogx( data1[!,:r], imag( data1[!,:Q] ), "r--" )

# figure(2)
# loglog( data1[!,:r], abs.(real( data1[!,:Q] .- data1[1,:Q] )), "b" )
# loglog( data1[!,:r], abs.(imag( data1[!,:Q] .- data1[1,:Q] )), "r--" )

n_terms = 7

result = zeros(Complex{Float64}, m, n_terms)
error = zeros(Float64, m, n_terms)
rel_error = zeros(Float64, m, n_terms)

result2 = zeros(Complex{Float64}, m)
error2 = zeros(Float64, m)
rel_error2 = zeros(Float64, m)

for i=1:m
    for j=1:n_terms
        print("$i.")
        result[i,j] = Q(n[i], r[i], ℒ[i]; n_terms = j)
        error[i,j] = abs( result[i,j] - data1[i,:Q]  )
        rel_error[i,j] = error[i,j]/abs( data1[i,:Q] )

        result2[i] = (-1)^n[i] * factorial(n[i]) * gamma(-n[i] .- r[i]) * (-μ[i])^r[i]  + zeta(1 + r[i])
        error2[i] = abs( result2[i] - data1[i,:Q]  )
        rel_error2[i] = error2[i]/abs( data1[i,:Q] )
    end
end
println("")

for i=1:5
    @printf(" r1 = %.20f, r2 = %.20f, r3 = %.20f\n", abs(result[i,1]), abs(result[i,2]) , abs(result[i,3])  )
end

fig = figure("Q_bench", figsize=(5,5))
clf()
d = 0.00
for j=1:n_terms
    loglog( data1[!,:r] .+ d, rel_error[:,j], "-"; label=@sprintf("%d terms",j) )
end
loglog( data1[!,:r] .+ d, rel_error2[:], "r-"; label="direct")
legend()
xlabel("tau")
ylabel("relative absolute error")


text:
using Test
using Graphs
using PeriodicGraphs
import PeriodicGraphs: ofs, vertex_permutation, dimensionality, LoopException, extended_gcd
import StaticArrays: SVector

# using Aqua
# Aqua.test_all(PeriodicGraphs)

@static if VERSION < v"1.4-"
    function only(t)
        length(t) == 1 || throw(ArgumentError("Not exactly 1 element"))
        return t[1]
    end
end

@testset "PeriodicVertex" begin
    @test_throws MethodError PeriodicVertex2D() # Always require a vertex identifier
    @test iszero(ofs(PeriodicVertex2D(1)))
    @test ndims(PeriodicVertex2D(1)) == 2
    @test PeriodicVertex2D(1) != PeriodicVertex3D(1)
    @test PeriodicVertex{4}(5) == PeriodicVertex{4}(5)
    @test PeriodicVertex1D(1, SVector{1,Int}(3)) == PeriodicVertex1D(1, [3])
    @test PeriodicVertex(3, (1,2)) == PeriodicVertex2D(3, (1,2)) == PeriodicVertex(3, SVector{2,Int}(1,2))
    @test PeriodicVertex[(1, (-1,0))] == [PeriodicVertex2D(1, (-1,0))] == PeriodicVertex2D[(1, (-1,0))]
    @test_throws MethodError PeriodicVertex(2, [4]) # Performance trap if the dimension is unknown
    @test_throws DimensionMismatch PeriodicVertex2D(1, (2,))
    @test_throws DimensionMismatch PeriodicVertex3D(3, (3, 1, 4, 0))
    @test ndims(PeriodicVertex(0, SVector{2,Int}(1,4))) == 2 # Inference from SVector
    @test ndims(PeriodicVertex{0}(1)) == 0 # No need of a second argument here
    @test PeriodicVertex{0}(1) < PeriodicVertex{0}(2)
    @test PeriodicVertex(2, (3,)) >  PeriodicVertex(2, (1,))
    @static if VERSION < v"1.6.0-DEV"
        @test string(PeriodicVertex3D(12, (1,0,-1))) == "PeriodicVertex{3}(12, (1,0,-1))"
        @test string(PeriodicVertex2D[(1, (0,0)), (1, (-1,0))]) == "PeriodicVertex{2}[(1, (0,0)), (1, (-1,0))]"
        @test string(PeriodicVertex1D(2, (1,))) == "PeriodicVertex{1}(2, (1,))"
    else
        @test string(PeriodicVertex3D(12, (1,0,-1))) == "PeriodicVertex3D(12, (1,0,-1))"
        @test string(PeriodicVertex2D[(1, (0,0)), (1, (-1,0))]) == "PeriodicVertex2D[(1, (0,0)), (1, (-1,0))]"
        @test string(PeriodicVertex1D(2, (1,))) == "PeriodicVertex1D(2, (1,))"
    end
    @test ofs(PeriodicVertex(1, (2,1))) == SVector{2,Int}(2,1)
end

@testset "PeriodicEdge" begin
    @test_throws LoopException PeriodicEdge2D(1, 1, (0, 0))
    @test_throws LoopException PeriodicEdge(2, PeriodicVertex3D(2))
    @test_throws LoopException PeriodicEdge{0}(3, 3, Int[])
    @test_throws LoopException PeriodicEdge(1, 1, ())
    try
        PeriodicEdge(1, PeriodicVertex2D(1))
        @test false
    catch e
        @test e isa LoopException
        io = IOBuffer()
        showerror(io, e)
        @test occursin("LoopException", String(take!(io)))
    end
    # There should be no default zero offset to prevent programmers from forgetting the offset
    @test_throws MethodError PeriodicEdge1D(1, 2)
    @test PeriodicEdge2D(1, 2, (0,0)) != PeriodicEdge3D(1, 2, (0,0,0))
    @test PeriodicEdge(2, PeriodicVertex2D(1, (1,3))) == PeriodicEdge2D(2, 1, (1,3))
    @test PeriodicEdge(1, 1, SVector{1,Int}(2)) == PeriodicEdge1D(1, 1, SVector{1,Int}(2))
    @test PeriodicEdge(2, 1, ()) == PeriodicEdge{0}((2, 1, SVector{0,Int}()))
    @test PeriodicEdge((1, 2, (1,0,0))) == PeriodicEdge3D((1, 2, (1,0,0)))
    @test convert(PeriodicEdge, (3, PeriodicVertex{1}(2, (1,)))) == convert(PeriodicEdge1D, (3, PeriodicVertex{1}(2, (1,))))
    @test PeriodicEdge1D[(1,2,SVector{1,Int}(3))] == PeriodicEdge[(1,2,(3,))] == [PeriodicEdge(1, 2, (3,))]
    @static if VERSION < v"1.6.0-DEV"
        @test string(PeriodicEdge3D[(1, 2, (1,0,0))]) == "PeriodicEdge{3}[(1, 2, (1,0,0))]"
        @test string(PeriodicEdge(3,4, (0,0))) == "PeriodicEdge{2}(3, 4, (0,0))"
        @test string(PeriodicEdge{1}(1, 1, (2,))) == "PeriodicEdge{1}(1, 1, (2,))"
    else
        @test string(PeriodicEdge3D[(1, 2, (1,0,0))]) == "PeriodicEdge3D[(1, 2, (1,0,0))]"
        @test string(PeriodicEdge(3,4, (0,0))) == "PeriodicEdge2D(3, 4, (0,0))"
        @test string(PeriodicEdge{1}(1, 1, (2,))) == "PeriodicEdge1D(1, 1, (2,))"
    end
    @test reverse(reverse(PeriodicEdge(1, 1, (2,3,1)))) == PeriodicEdge3D(1, 1, (2,3,1))
    @test reverse(PeriodicEdge(1, 2, (1,-2))) == PeriodicEdge(2, 1, (-1,2))
    @test reverse(PeriodicEdge(3, 3, (1,0))) == PeriodicEdge(3, 3, (-1,0))
    @test src(PeriodicEdge(1, 2, (-1, 0))) == 1
    @test dst(PeriodicEdge(2, 2, (1,))) == 2
    @test ofs(PeriodicEdge(1, 1, (0,1,3))) == SVector{3,Int}(0,1,3)
    @test PeriodicEdge2D(1, 1, (2,3)) < PeriodicEdge2D(1, 1, (3,0)) < PeriodicEdge2D(1, 2, (3,0)) < PeriodicEdge2D(2, 2, (3,0))
    @test ndims(PeriodicEdge(1, 2, ())) == 0ndims(PeriodicEdge(1, 2, ())) == 0
    @test ndims(PeriodicEdge(1, 1, (2,0,0,0))) == 4
end

@testset "Empty graphs" begin
    @test_throws MethodError PeriodicGraph() # Missing dimension
    g = PeriodicGraph3D()
    @test ne(g) == nv(g) == 0
    @test ndims(g) == 3
    @test !rem_edge!(g, PeriodicEdge3D(1, 2, (0,0,0)))
    @test_throws ArgumentError rem_vertex!(g, 0)
    @test_throws ArgumentError rem_vertex!(g, 1)
    @test isempty(rem_vertices!(g, Int[]))
    @test isempty(edges(g))
    @test isempty(vertices(g))
    gg = vertex_permutation(g, [])
    @test ne(gg) == nv(gg) == 0 && isempty(edges(gg))
    @test g == gg
    gg, vlist = induced_subgraph(g, Int[])
    @test isempty(vlist)
    @test ne(gg) == nv(gg) == 0 && isempty(edges(gg))
    @test g == gg
    gg = offset_representatives!(g, Int[])
    @test ne(gg) == nv(gg) == 0 && isempty(edges(gg))
    @test ne(g) == nv(g) == 0 && isempty(edges(g))
    @test g == gg
    @test_throws DimensionMismatch swap_axes!(g, [2,1])
    gg = swap_axes!(g, Int[1,2,3])
    @test ne(gg) == nv(gg) == 0 && isempty(edges(gg))
    @test ne(g) == nv(g) == 0 && isempty(edges(g))
    @test g == gg
    @test isempty(connected_components(g))
    @test edgetype(g) == PeriodicEdge3D
    @test eltype(g) == PeriodicVertex3D
    @test edges(g) isa PeriodicGraphs.PeriodicEdgeIter{3}
    @test_throws BoundsError coordination_sequence(g, 1, 0)
    @test_throws BoundsError coordination_sequence(g, 1, 0) # once width is set
    @test !has_vertex(g, 0)
    @test !has_vertex(g, 1)
end

@testset "Basic graph construction and modification" begin
    @test_throws MethodError PeriodicGraph(2) # missing dimension
    g = PeriodicGraph3D()

    @testset "add_vertex!, rem_vertex!" begin
        @test add_vertex!(g)
        @test ne(g) == 0 && nv(g) == 1
        @test g == PeriodicGraph3D(1)
        @test add_vertex!(g)
        @test ne(g) == 0 && nv(g) == 2
        @test g == PeriodicGraph3D(2)
        @test_throws ArgumentError rem_vertex!(g, 0)
        @test_throws ArgumentError rem_vertex!(g, 3)
        @test rem_vertex!(g, 1)
        @test ne(g) == 0 && nv(g) == 1
        @test g == PeriodicGraph3D(1)
        @test rem_vertex!(g, 1)
        @test ne(g) == 0 && nv(g) == 0
        @test g == PeriodicGraph3D(0)
    end

    @testset "add_vertices!, rem_vertices!" begin
        @test isempty(edges(g))
        @test isempty(vertices(g))
        @test_throws ArgumentError rem_vertices!(g, Int[1])
        @test isempty(rem_vertices!(g, Int[]))
        @test add_vertices!(g, 0) == 0
        @test nv(g) == ne(g) == 0 && g == PeriodicGraph3D()
        @test add_vertices!(g, 1) == 1
        @test ne(g) == 0 && nv(g) == 1 && g == PeriodicGraph3D(1)
        @test only(vertices(g)) == 1
        @test rem_vertices!(g, [1]) == Int[]
        @test add_vertices!(g, 5) == 5
        @test ne(g) == 0 && nv(g) == 5 && g == PeriodicGraph3D(5)
        @test vertices(g) == 1:5
        @test issetequal(rem_vertices!(g, [2,3]), [1,4,5])
        @test ne(g) == 0 && nv(g) == 3 && g == PeriodicGraph3D(3)
        @test vertices(g) == 1:3
        @test rem_vertex!(g, 2)
        @test ne(g) == 0 && nv(g) == 2 && g == PeriodicGraph3D(2)
        @test vertices(g) == [1,2]
        @test rem_vertices!(g, [1,2]) == Int[]
        @test add_vertices!(g, 5) == 5
        @test rem_vertices!(g, [2,4], true) == [1,3,5]
        @test g == PeriodicGraph3D(3)
        @test rem_vertices!(g, [1,2,3], true) == Int[]
        @test g == PeriodicGraph3D()
    end

    @testset "Single edge" begin
        @test isempty(edges(g))
        @test isempty(vertices(g))
        @test !rem_edge!(g, 1, PeriodicVertex(2, (0,0,0)))
        @test g == PeriodicGraph3D()
        @test !add_edge!(g, PeriodicEdge(1, 1, (-1,0,0)))
        @test g == PeriodicGraph3D()
        @test !has_edge(g, 1, 1)
        @test !has_edge(g, PeriodicEdge3D(1, 2, (0,0,0)))
        @test add_vertex!(g)
        @test !has_edge(g, 1, 1)
        @test !has_edge(g, PeriodicEdge3D(1, 1, (0,1,1)))
        @test add_edge!(g, 1, PeriodicVertex(1, (-1,0,0)))
        @test has_edge(g, 1, 1)
        @test !add_edge!(g, PeriodicEdge(1, 1, (-1,0,0)))
        @test has_edge(g, 1, 1)
        @test has_edge(g, PeriodicEdge3D(1, 1, (1,0,0)))
        @test has_edge(g, 1, PeriodicVertex(1, (-1,0,0)))
        @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 1, (1,0,0))])
        @test g != PeriodicGraph3D(1)
        @test rem_edge!(g, PeriodicEdge(1, 1, (-1,0,0)))
        @test g == PeriodicGraph3D(1)
        @test !has_edge(g, 1, 1)
        @test !has_edge(g, 1, PeriodicVertex(1, (1,0,0)))
        @test !has_edge(g, PeriodicEdge3D(1, 1, (-1,0,0)))
        @test add_edge!(g, 1, PeriodicVertex(1, (2,-1,0)))
        @test !add_edge!(g, PeriodicEdge(1, 1, (-2,1,0)))
        @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 1, (-2,1,0))])
        @test has_edge(g, 1, 1)
        @test has_edge(g, PeriodicEdge3D(1, 1, (2,-1,0)))
        @test !has_edge(g, PeriodicEdge3D(1, 1, (1,0,0)))
        @test rem_edge!(g, 1, PeriodicVertex(1, (2,-1,0)))
        @test g == PeriodicGraph3D(1)
        @test add_edge!(g, PeriodicEdge(1, 1, (1,0,0)))
        @test rem_vertex!(g, 1)
        @test g == PeriodicGraph3D(0)
        @test add_vertices!(g, 2) == 2
        @test !has_edge(g, 1, 2)
        @test g == PeriodicGraph3D(2)
        @test add_edge!(g, PeriodicEdge(1, 2, (0,0,0)))
        @test has_edge(g, 1, 2)
        @test has_edge(g, PeriodicEdge(1, 2, (0,0,0)))
        @test !has_edge(g, PeriodicEdge(1, 2, (1,0,0)))
        @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 2, (0,0,0))])
        @test g != PeriodicGraph3D(2)
        @test rem_edge!(g, PeriodicEdge(1, 2, (0,0,0)))
        @test !has_edge(g, 1, 2)
        @test g == PeriodicGraph3D(2)
        @test add_edge!(g, 1, PeriodicVertex(2, (-1,0,0)))
        @test g == PeriodicGraph3D(PeriodicEdge3D[(1, 2, (-1,0,0))])
        @test g == PeriodicGraph3D(PeriodicEdge3D[(2, 1, (1,0,0))])
        @test g != PeriodicGraph3D(PeriodicEdge3D[(1, 2, (1,0,0))])
        @test has_edge(g, 1, 2)
        @test has_edge(g, 1, PeriodicVertex(2, (-1,0,0)))
        @test has_edge(g, PeriodicEdge(2, 1, (1,0,0)))
        @test rem_vertex!(g, 1)
        @test g == PeriodicGraph3D(1)
        @test !has_edge(g, 1, 2)
        @test !has_edge(g, 1, 1)
        @test add_vertex!(g)
        @test nv(g) == 2
        @test add_edge!(g, PeriodicEdge(1, 2, (0,1,1)))
        @test ne(g) == 1
        @test rem_vertices!(g, [2]) == [1]
        @test ne(g) == 0
        @test nv(g) == 1
    end
end

@testset "String graph construction" begin
    @test_throws MethodError length(PeriodicGraphs.KeyString{Int}("3 1 2 0 0 0"))
    @test PeriodicGraph("0 1 2 1 3") == PeriodicGraph(PeriodicEdge{0}[(1, 2, ()), (1, 3, ())])
    @test PeriodicGraph("3") == PeriodicGraph3D("3") == PeriodicGraph3D()
    @test ne(PeriodicGraph("3  1 1  1 0 0  1 1 -1 0 0  1 1 3 0 0")) == 2
    @test nv(PeriodicGraph("2  11 12 0 0")) == 12
    @test collect(PeriodicGraphs.KeyString{Int}("3 2 1 0 1 -1")) == [3, 2, 1, 0, 1, -1]
    @test_throws ArgumentError PeriodicGraph2D("")
    @test_throws ArgumentError PeriodicGraph("")
    @test_throws ArgumentError PeriodicGraph3D("a")
    @test_throws ArgumentError PeriodicGraph2D("2 a")
    @test_throws DimensionMismatch PeriodicGraph2D("4 1 2 0 0 0 1")
    @test_throws DimensionMismatch PeriodicGraph3D("2  1 2  0 0 0   1 3  0 1 0")
    @test_throws ArgumentError PeriodicGraph("1 2 1")
    @test_throws ArgumentError PeriodicGraph{4}("4 1")
    @test_throws ArgumentError PeriodicGraph("2 1 1 0 1 0")
    @test_throws LoopException PeriodicGraph("3 1 1 0 0 0")
    @test string(PeriodicGraph("2 1 2 0 0 2 3 1 0")) == "2 1 2 0 0 2 3 1 0"
    @test string(PeriodicGraph("1    23 4  5    10  9 43")) == "1 4 23 -5 9 10 -43"
    @test string(PeriodicGraph{5}()) == "5"
end

@testset "Basic graph representation and functions" begin
    @test zero(PeriodicGraph{6}) == PeriodicGraph{6}()
    io = IOBuffer()
    show(io, PeriodicGraph("2 1 2 0 0 2 3 1 0"))
    g1 = eval(Meta.parse(String(take!(io))))
    @test g1 isa PeriodicGraph2D
    @test string(g1) == "2 1 2 0 0 2 3 1 0"
    g2 = PeriodicGraph(g1)
    g3 = deepcopy(g1)
    @test g1 == g2 == g3
    rem_vertex!(g1, 1)
    @test g1 != g2 == g3
    @test add_edge!(g2, PeriodicEdge(1, 2, (-1,0)))
    @test add_edge!(g2, 1, PeriodicVertex(1, (1,1)))
    @test add_edge!(g2, PeriodicEdge(1, 3, (0,0)))
    @test find_edges(g2, 1, 1) == PeriodicVertex2D[(1, (-1,-1)), (1, (1,1))]
    @test find_edges(g2, 1, 2) == PeriodicVertex2D[(2, (-1,0)), (2, (0,0))]
    @test find_edges(g2, 3, 2) == PeriodicVertex2D[(2, (-1,0))]
    @test inneighbors(g2, 2) == outneighbors(g2, 2)
    @test has_vertex(g1, 1)
end

@testset "Complex graphs modifications" begin
    edgs = PeriodicEdge3D[(1, 1, (1,0,0)),  (1, 2, (0,0,0)), (1, 3, (0,0,1)),
                          (1, 3, (0,1,-1)), (2, 3, (0,0,0)), (2, 4, (0,-1,0))]
    g = PeriodicGraph3D(edgs)
    @test nv(g) == 4
    @test ne(g) == length(edges(g)) == length(edgs)
    @test collect(edges(g)) == edgs
    insert!(edgs, 5, PeriodicEdge(1, 4, (-1,-1,0)))
    @test add_edge!(g, edgs[5])
    @test ne(g) == length(edges(g)) == length(edgs)
    @test collect(edges(g)) == edgs
    @test rem_edge!(g, edgs[4])
    deleteat!(edgs, 4)
    @test ne(g) == length(edges(g)) == length(edgs)
    @test collect(edges(g)) == edgs
    @test g[[1,2,3,4]] == g
    gg = g[[1,4,2,3]]
    @test gg != g
    @test gg[[1,3,4,2]] == g
    gg = g[[3,1]]
    @test nv(gg) == 2
    @test collect(edges(gg)) == PeriodicEdge3D[(1, 2, (0,0,-1)), (2, 2, (1,0,0))]
    @test_throws ArgumentError g[[3,3]]
    @test_throws ArgumentError offset_representatives!(gg, [2])
    gcopy = deepcopy(g)
    @test offset_representatives!(g, [0,0,0,0]) == gcopy == g
    ggg = offset_representatives!(gg, [(1,-1,1),0])
    @test ggg == gg
    @test collect(edges(gg)) == PeriodicEdge3D[(1, 2, (1,-1,0)), (2, 2, (1,0,0))]
    gcopy = deepcopy(gg)
    offset_representatives!(gg, [(2,1,-3), (2,1,-3)])
    @test gg == gcopy
    @test_throws DimensionMismatch swap_axes!(g, [1,3])
    ggg = swap_axes!(gg, (1,2,3))
    @test gcopy == gg == ggg
    ggg = swap_axes!(gg, (2,1,3))
    @test gcopy != gg == ggg
    @test collect(edges(gg)) == PeriodicEdge3D[(1, 2, (-1,1,0)), (2, 2, (0,1,0))]
    @test edges(gg) < edges(gcopy) # lexicographical ordering on the list of edges
end

@testset "Neighbors" begin
    g = PeriodicGraph3D(1)
    @test add_vertices!(g, 2) == 2
    @test add_edge!(g, PeriodicEdge3D(2, 3, (0, 0, 1)))
    @test only(neighbors(g, 2)) == PeriodicVertex3D(3, (0, 0, 1))
    @test only(neighbors(g, 3)) == PeriodicVertex3D(2, (0, 0, -1))
    @test !has_edge(g, 2, 2)
    @test add_edge!(g, PeriodicEdge3D(2, 2, (1, 0, 0)))
    @test only(neighbors(g, 3)) == PeriodicVertex3D(2, (0, 0, -1))
    @test length(neighbors(g, 2)) == 3
    @test has_edge(g, 2, 2)
    @test neighborhood(g, 1, 1) == pushfirst!(collect(neighbors(g, 1)), PeriodicVertex3D(1, (0, 0, 0)))
    @test coordination_sequence(g, 1, 3) == zeros(Int, 3)
    @test coordination_sequence(g, 3, 6) == [1, 2, 4, 4, 4, 4]

    g = PeriodicGraph{0}(4)
    @test add_edge!(g, PeriodicEdge{0}(1, 2, ()))
    @test add_edge!(g, PeriodicEdge(1, 3, ()))
    @test add_edge!(g, 2, PeriodicVertex(4, ()))
    @test coordination_sequence(g, 1, 10) == [2, 1, 0, 0, 0, 0, 0, 0, 0, 0,]

    g = PeriodicGraph1D(3)
    @test add_edge!(g, PeriodicEdge1D(1, 1, (1,)))
    @test add_edge!(g, PeriodicEdge(1, 2, (0,)))
    @test add_edge!(g, 2, PeriodicVertex(3, (0,)))
    @test add_edge!(g, 3, PeriodicVertex1D(3, (1,)))
    @test coordination_sequence(g, 1, 5) == coordination_sequence(g, 3, 5) == [3, 5, 6, 6, 6]
    @test coordination_sequence(g, 2, 5) == [2, 4, 6, 6, 6]

    g = PeriodicGraph{4}(1)
    @test add_edge!(g, PeriodicEdge{4}(1, 1, (1,0,0,0)))
    @test add_edge!(g, PeriodicEdge(1, 1, (0,1,0,0)))
    @test add_edge!(g, 1, PeriodicVertex(1, (0,0,1,0)))
    @test add_edge!(g, 1, PeriodicVertex{4}(1, (0,0,0,1)))
    @test coordination_sequence(g, 1, 30) == [8i*(i^2+2)/3 for i in 1:30] # sequence A008412 of the OEIS
end

@testset "Graph reduction" begin
    g = PeriodicGraph3D("3  1 2  0 0 0  1 3  0 1 0   2 3  1 0 -1
                            4 4  0 0 1  4 5  0 -1 0  5 6  0 0 0")
    c = cellgraph(g)
    p = periodiccellgraph(g)
    @test nv(g) == nv(p) == nv(c)
    @test issubset(edges(c), edges(p))
    for e in edges(p)
        @test has_edge(g, src(e), dst(e))
    end
    @test connected_components(p) == connected_components(g)
end

@testset "Periodic vertex hash" begin
    ## dim = 0
    @test PeriodicGraphs.hash_position(SVector{0,Int}()) == 0
    ## dim = 1
    @test PeriodicGraphs.hash_position(SVector{1,Int}(0)) == 0
    @test PeriodicGraphs.hash_position(SVector{1,Int}(-1)) == 1
    @test PeriodicGraphs.hash_position(SVector{1,Int}(1)) == 2
    @test PeriodicGraphs.hash_position(SVector{1,Int}(-2)) == 3
    @test PeriodicGraphs.hash_position(SVector{1,Int}(2)) == 4
    @test PeriodicGraphs.hash_position(SVector{1,Int}(3)) == 6
    ## dim = 2
    n::Int = 4
    seen = falses((2n-1)^2)
    for i in -n+1:n-1, j in -n+1:n-1
        x = PeriodicGraphs.hash_position(SVector{2,Int}(i,j))+1
        @test !seen[x]
        seen[x] = true
    end
    @test all(seen)
    append!(seen, falses((2n+1)^2-(2n-1)^2))
    for i in (-n,n), j in -n:n
        x = PeriodicGraphs.hash_position(SVector{2,Int}(i,j))+1
        @test !seen[x]
        seen[x] = true
    end
    for i in -n+1:n-1, j in (-n,n)
        x = PeriodicGraphs.hash_position(SVector{2,Int}(i,j))+1
        @test !seen[x]
        seen[x] = true
    end
    @test all(seen)
    ## dim = 3
    seen = falses((2n-1)^3)
    for i in -n+1:n-1, j in -n+1:n-1, k in -n+1:n-1
        x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1
        @test !seen[x]
        seen[x] = true
    end
    @test all(seen)
    append!(seen, falses((2n+1)^3-(2n-1)^3))
    for i in (-n,n), j in -n:n, k in -n:n
        x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1
        @test !seen[x]
        seen[x] = true
    end
    for i in -n+1:n-1, j in (-n,n), k in -n:n
        x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1
        @test !seen[x]
        seen[x] = true
    end
    for i in -n+1:n-1, j in -n+1:n-1, k in (-n,n)
        x = PeriodicGraphs.hash_position(SVector{3,Int}(i,j,k))+1
        @test !seen[x]
        seen[x] = true
    end

    @test all(seen)
    g = PeriodicGraph3D([PeriodicEdge3D(1, 1, (0, 0, 1)),
                         PeriodicEdge3D(1, 1, (1, 1, 0))])
    @test neighborhood(g, 1, 1) == PeriodicVertex3D[(1, (0, 0, 0)),
            (1, (-1, -1, 0)), (1, (0, 0, -1)), (1, (0, 0, 1)), (1, (1, 1, 0))]
    @test neighborhood(g, 1, 2) == PeriodicVertex3D[(1, (0, 0, 0)), (1, (-1, -1, 0)),
            (1, (0, 0, -1)), (1, (0, 0, 1)), (1, (1, 1, 0)), (1, (-2, -2, 0)),
            (1, (-1, -1, -1)), (1, (-1, -1, 1)), (1, (0, 0, -2)), (1, (1, 1, -1)),
            (1, (0, 0, 2)), (1, (1, 1, 1)), (1, (2, 2, 0))]
end

@testset "Edge iteration" begin
    g::PeriodicGraph3D = PeriodicGraph3D([PeriodicEdge3D(2, 3, (-1, 0, 0)),
                                          PeriodicEdge3D(2, 3, (0, 1, 0)),
                                          PeriodicEdge3D(4, 4, (-1, 1, 0)),
                                          PeriodicEdge3D(1, 3, (0, 0, 1)),
                                          PeriodicEdge3D(3, 2, (0, 0, 0))])
    expected = PeriodicEdge3D[(1, 3, (0, 0, 1)), (2, 3, (-1, 0, 0)),
                              (2, 3, (0, 0, 0)), (2, 3, (0, 1, 0)),
                              (4, 4, (1, -1, 0))]
    @test collect(edges(g)) == expected
    @test PeriodicEdge3D(4, 4, (1, -1, 0)) ∈ edges(g)
    add_vertex!(g)
    @test collect(edges(g)) == expected
    @test length(edges(g)) == ne(g)
end

@testset "Vertex removal" begin
    g::PeriodicGraph3D = PeriodicGraph3D([PeriodicEdge3D(1, 3, (0, 0, 1)),
                                          PeriodicEdge3D(2, 3, (0, 1, 0)),
                                          PeriodicEdge3D(2, 3, (0, -1, 0)),
                                          PeriodicEdge3D(3, 2, (0, 0, 0)),
                                          PeriodicEdge3D(4, 4, (1, 1, 0))])
    gg = deepcopy(g)
    @test rem_vertices!(g, Int[]) == [1, 2, 3, 4]
    @test g == gg
    @test rem_vertices!(g, [2, 4], false) == [1, 3]
    @test nv(g) == 2
    @test ne(g) == 1
    @test neighbors(g, 1) == PeriodicVertex3D[(2, (0, 0, 1))]
    @test neighbors(g, 2) == PeriodicVertex3D[(1, (0, 0, -1))]
    @test rem_vertex!(gg, 2)
    @test rem_vertex!(gg, 2)
    @test g == gg

    g1 = PeriodicGraph("2
            1 1  1 0
            1 2  0 0
            1 2  0 1
            1 3  0 0
            2 3  1 0
            2 4  0 1
            3 4  0 0
            2 5  0 0
            4 5  1 1")
    g2 = deepcopy(g1)
    @test rem_vertices!(g1, [2,4], false) == [1, 5, 3]
    @test g1 == PeriodicGraph("2 1 1 1 0 1 3 0 0")
    @test rem_vertices!(g2, [2,4], true) == [1, 3, 5]
    gref = PeriodicGraph("2 1 1 1 0 1 2 0 0")
    @test add_vertex!(gref)
    @test g2 == gref
end

function equivalent_dict(d1, d2)
    keys(d1) == keys(d2) || return false
    for k in keys(d1)
        s1 = Set(Set(x) for x in d1[k])
        s2 = Set(Set(x) for x in d2[k])
        s1 == s2 || return false
    end
    return true
end

@testset "Dimensionality" begin
    @test dimensionality(PeriodicGraph{0}(5)) == Dict(0 => [collect(1:5)])
    g::PeriodicGraph3D = PeriodicGraph3D([PeriodicEdge3D(1, 3, (0, 0, 1)),
                                          PeriodicEdge3D(2, 3, (0, 1, 0)),
                                          PeriodicEdge3D(2, 3, (0, -1, 0)),
                                          PeriodicEdge3D(4, 4, (1, 1, 0))])
    @test equivalent_dict(dimensionality(g), Dict(1 => connected_components(g)))
    @test add_edge!(g, PeriodicEdge(3, 2, (0,0,0)))
    @test equivalent_dict(dimensionality(g), Dict(1 => connected_components(g)))
    @test add_edge!(g, PeriodicEdge(2, 2, (0,0,1)))
    @test equivalent_dict(dimensionality(g), Dict(1 => [[4]], 2 => [[1,2,3]]))
    @test rem_edge!(g, PeriodicEdge(1, 3, (0,0,1)))
    @test equivalent_dict(dimensionality(g), Dict(0 => [[1]], 1 => [[4]], 2 => [[2,3]]))
    @test rem_edge!(g, 3, PeriodicVertex(2, (0,0,0)))
    @test equivalent_dict(dimensionality(g), Dict(0 => [[1]], 1 => [[4]], 2 => [[2,3]]))
    @test rem_edge!(g, PeriodicEdge(2, 3, (0,1,0)))
    @test equivalent_dict(dimensionality(g), Dict(0 => [[1]], 1 => [[2,3],[4]]))
    @test add_edge!(g, 2, PeriodicVertex3D(3, (1,0,-1)))
    @test equivalent_dict(dimensionality(g), Dict(0 => [[1]], 1 => [[4]], 2 => [[2,3]]))
    @test add_edge!(g, 2, PeriodicVertex(3, (2,0,-1)))
    @test equivalent_dict(dimensionality(g), Dict(0 => [[1]], 1 => [[4]], 3 => [[2,3]]))
end

@testset "Dimension change" begin
    g = PeriodicGraph("3
            1 2 0 0 0
            2 3 1 0 0
            3 4 -1 0 0
            4 1 0 0 1
            4 4 0 0 1")
    @test equivalent_dict(dimensionality(g), Dict(1 => [[1,2,3,4]]))
    gg = change_dimension(PeriodicGraph1D, g)
    @test string(gg) == "1 1 2 0 1 4 -1 2 3 0 3 4 0 4 4 1"
    @test change_dimension(PeriodicGraph2D, gg) == change_dimension(PeriodicGraph2D, g) ==
        PeriodicGraph("2 1 2 0 0 1 4 -1 0 2 3 0 0 3 4 0 0 4 4 1 0")
    @test_throws DimensionMismatch change_dimension(PeriodicGraph{0}, g)

    s = [14,28,-17,-34,12]
    d, λ = extended_gcd(s)
    @test d == 1
    @test reduce(+, s .* λ) == d
end


text:
using SafeTestsets, Test

@safetestset "CanonicalSymplecticStructures" begin
    include("test_CanonicalSymplecticStructures.jl")
end

@safetestset "FirstPoincareInvariants" begin
    include("FirstPoincareInvariants/test_FirstPoincareInvariants.jl")
end

@safetestset "SecondPoincareInvariants" begin
    include("SecondPoincareInvariants/test_SecondPoincareInvariants.jl")
end


text:
function test_l2_constraint(backend::Backend, T, eps)
  println("-- Testing L2 constraint on $(typeof(backend)){$T}...")
  # this simulates a convolutional filter and applies
  # the l2 constraint to it
  n_filters = 5
  coef = 0.2
  param = rand(T, 2,3,4,n_filters) - 0.5
  param_after = zeros(T, size(param))
  param_blob = make_blob(backend, param)

  cons = L2Cons(coef)
  constrain!(backend, cons, param_blob)
  copy!(param_after, param_blob)
  param_after = reshape(param_after, size(param))
  for f=1:n_filters
    norm2 = vecnorm(param_after[:, :, :, f])
    @test norm2 <= coef + eps
  end

  # this is the same as above but for fully connected weights
  n_input = 10
  n_out   = 12
  param = rand(T, n_input,n_out) - 0.5
  param_after = zeros(T, size(param))
  param_blob = make_blob(backend, param)

  cons = L2Cons(coef)
  constrain!(backend, cons, param_blob)
  copy!(param_after, param_blob)
  param_after = reshape(param_after, size(param))
  for f=1:n_out
    norm2 = vecnorm(param_after[:, f])
    @test norm2 <= coef + eps
  end

  # The case for bias
  len = 10
  param = rand(T, len) - 0.5
  param_after = similar(param)
  param_blob = make_blob(backend, param)
  cons = L2Cons(coef)
  constrain!(backend, cons, param_blob)
  copy!(param_after, param_blob)
  @test vecnorm(param_after) <= coef+eps
end

function test_l2_constraint(backend::Backend)
  test_l2_constraint(backend, Float32, 1e-5)
  test_l2_constraint(backend, Float64, 1e-5)
end

if test_cpu
  test_l2_constraint(backend_cpu)
end
if test_gpu
  test_l2_constraint(backend_gpu)
end


text:
 function checkzerovar(m::LinearMixedModel)
    ff = deepcopy(m.f)
    # set fixed term=true, reterm=false in rhs of formula

    fixterms=map(x-> !(Meta.isexpr(x,:call) && x.args[1] == :|), m.mf.terms.terms)

    retrms = filter(x->Meta.isexpr(x,:call) && x.args[1] == :|, m.mf.terms.terms)
    eachtrms=[aa.args[2] for aa in retrms]
    grps = Any[t.args[3] for t in retrms]
    ugrps = unique(grps)

    #check the numer of terms in each group
    trmsleng=zeros(length(eachtrms))

    for k in 1:length(eachtrms)
      if isa(eachtrms[k],Int64)
        trmsleng[k]=1
      elseif  isa(eachtrms[k],Symbol)
          trmsleng[k]=2
      elseif  eachtrms[k].args[1] == :+ && eachtrms[k].args[2] == 0
        trmsleng[k]=length(eachtrms[k].args)-2
      else
        trmsleng[k]=length(eachtrms[k].args)
      end
    end

    #find λ and their lowe bounds , compare them,
    # set true if boundary conditions are met or set false otherwise
    est=map(MixedModels.θ,m.λ)
    lowbd=map(MixedModels.lower,m.λ)

    if all(trmsleng .== 1.0)

      if length(ugrps) == length(grps)
        remainterms1=map(all,Any[l .> k  for (l,k) in zip(est,lowbd)])

      else # length(ugrps) < length(grps)
        aa=[l .> k  for (l,k) in zip(est,lowbd)]
        remainterms1=aa[1]
        for i in 2:length(aa)
         remainterms1=[remainterms1,aa[i]]
        end
      end

        # set reterm be false (will be excluded) if  boundary conditions are not met
      fixterms[fixterms .== false] = remainterms1
      # chaneg rhs of formula by excluding reterms for which boundary conditions are not met
      ffm=m.mf.terms.terms[convert(Array{Bool,1},fixterms)]

      if all(remainterms1)
          ff=ff
      elseif length(ffm)==0
          warn("Est. of variance components of all terms are on the boundary")
          ff=y~1
      elseif (length(ffm)==1 && !any(remainterms1))
          warn("Est. of variance components of all terms are on the boundary")
          ff=eval(y~$(ffm[1]))
      elseif (length(ffm)==1 && any(remainterms1))
          ff=eval(y~$(ffm[1]))
      # Note that
      # if there are more than 1 term, the first element of formula.rhs.args is :+ with
      # head function :call
      else
          a=Expr(:call,tuple ( [:+,ffm]...)...)
          ff=eval(y~$a)
      end

          #returm modofied formula
      return(ff)

    else

      ffmt=deepcopy(m.mf.terms.terms)
      numfixterm= sum(fixterms .== true)

      if length(ugrps) == length(grps)

          aa=[l .> k  for (l,k) in zip(est,lowbd)]

          for i in 1:length(aa)
              if length(aa[i])==1

                if aa[i][1]==true
                  fixterms[numfixterm+i]=true
                else #  aa[i][1]==false
                  fixterms[numfixterm+i]=false
                end

         #  especially for model m5
              elseif length(aa[i])==3

                   if aa[i][1] == true && aa[i][3] == true
                     fixterms[numfixterm+i]=true
                   elseif aa[i][1] == false && aa[i][3] == true
                     fixterms[numfixterm+i]=true
                     ffmt[2]= :(0+time | id)
                   elseif aa[i][1] == true && aa[i][3] == false
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :(1 | id)
                   else
                     fixterms[numfixterm+i]=false
                   end
         #  especially for model m7

              elseif length(aa[i])==6

                   if aa[i][1] == true && aa[i][4] == true  && aa[i][6] == true
                     fixterms[numfixterm+i]=true

                   elseif aa[i][1] == false && aa[i][4] == true  && aa[i][6] == true
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :((0+time + time2) | id)
                   elseif aa[i][1] == true && aa[i][4] == false  && aa[i][6] == true
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :((1+ time2) | id)
                   elseif aa[i][1] == true && aa[i][4] == true  && aa[i][6] == false
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :((1+ time) | id)

                   elseif aa[i][1] == false && aa[i][4] == false  && aa[i][6] == true
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :((0+ time2) | id)
                   elseif aa[i][1] == true && aa[i][4] == false  && aa[i][6] == false
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :(1 | id)
                   elseif aa[i][1] == false && aa[i][4] == true  && aa[i][6] == false
                     fixterms[numfixterm+i]=true
                     ffmt[2] = :((0+ time) | id)

                    else
                      fixterms[numfixterm+i]=false
                    end

               end
           end
        end

        ffm=ffmt[convert(Array{Bool,1},fixterms)]

        a=Expr(:call,tuple ( [:+,ffm]...)...)
        ff=eval(y~$a)

      return(ff)

    end

end
 ## Return a block in the Zt matrix from one term.
function ztblk(m::Matrix,v)
    nr,nc = size(m)
    nblk = maximum(v)
    NR = nr*nblk                        # number of rows in result
    cf = length(m) < typemax(Int32) ? int32 : int64 # conversion function
    SparseMatrixCSC(NR,nc,
                    cf(cumsum(vcat(1,fill(nr,(nc,))))), # colptr
                    cf(vec(reshape([1:NR],(nr,int(nblk)))[:,v])), # rowval
                    vec(m))            # nzval
end
ztblk(m::Matrix,v::PooledDataVector) = ztblk(m,v.refs)

##  Conditional Likelihood
function condll(m::LinearMixedModel)
    n,p,q,k = size(m)
    sig = std(m)[end][1]
    cll=0.0
    for i in 1:n
        cll += logpdf(Normal(m.μ[i],sig),m.y[i])
    end
    cll
end

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
include("activation.jl")

type LogisticRegression
    n_input::Int
    n_output::Int
    η::Float64  # learning rate
    maxiter::Int
    error_threshold::Float64
    W::AA
    b::AA  # mean along axis 2
    X::AA
    Ŷ::AA  # Output of LogisticRegression Layer

    function LogisticRegression(n_input::Int,
                                   n_output::Int,
                                   learning_rate::Float64,
                                   maxiter::Int,
                                   error_threshold::Float64 = 1e-3,
                                   W::AA = zeros(n_output, n_input),
                                   b::AA = zeros(n_output))
        X = zeros(0, 0)
        Ŷ = zeros(0, 0)
        new(n_input, n_output, learning_rate, maxiter, error_threshold,
            W, b, X, Ŷ)
    end
end

forward{T}(W::AA{T, 2}, x::AA{T, 1}, b::AA{T, 1}) = sigmoid(W*x + b)

∇bias(y::Float64, ŷ::Float64) = (ŷ - y) * ŷ * (1 - ŷ)

∇weights(x::AA, y::Float64, ŷ::Float64) = (ŷ - y) * ŷ * (1 - ŷ) * x

function ∇bias(y::AA, ŷ::AA)
    assert(size(y) == size(ŷ))
    ∇b = zeros(size(y))
    for i in 1:size(y, 1)
        ∇b[i] = ∇bias(y[i], ŷ[i])
    end
    return ∇b
end

function ∇weights(x::AA, y::AA, ŷ::AA)
    assert(size(y) == size(ŷ))
    ∇W = zeros(size(y, 1), size(x, 1))
    for i in 1:size(y, 1)
        ∇W[i, :] = ∇weights(x, y[i], ŷ[i])
    end
    return ∇W
end

function forward{T}(W::AA{T, 2}, X::AA{T, 2}, b::AA{T, 1})
    N = size(X, 2)
    Ŷ = zeros(size(W, 1), N)
    for i in 1:N
        Ŷ[:, i] = forward(W, view(X, :, i), b)
    end
    return Ŷ
end

function forward!(layer::LogisticRegression, X::AA)
    layer.X = X
    layer.Ŷ = forward(layer.W, X, layer.b)
end

function backward!(layer::LogisticRegression, Y::AA)
    """
    Parameters
    ----------
    Y : AbstractArray
    Expected output represented as one-hot vectors
    """

    for i in 1:size(Y, 2)
        x = view(layer.X, :, i)
        y = view(Y, :, i)
        ŷ = view(layer.Ŷ, :, i)

        layer.W -= layer.η * ∇weights(x, y, ŷ)
        layer.b -= layer.η * ∇bias(y, ŷ)
    end
end

function fit!(layer::LogisticRegression, X::AA, Y::AA)
    assert(size(X, 2) == size(Y, 2))  # TODO rewrite with error()

    # TODO show progress
    for i in 1:layer.maxiter
        Ŷ = forward!(layer, X)
        backward!(layer, Y)
        if sum((Y - Ŷ) .^ 2) < layer.error_threshold
            return layer
        end
    end
    return layer
end

predict(layer::LogisticRegression, X::AA) = forward(layer.W, X, layer.b)

max_stars_repo_name: IshitaTakeshi/NeuralNetwork

max_issues_repo_name: IshitaTakeshi/NeuralNetwork

max_forks_repo_name: IshitaTakeshi/NeuralNetwork


text:
@testset "Test load_iris" begin
    xtmp, ytmp, vars = load_iris()

    @test typeof(xtmp) <: Matrix{Float64}
    @test typeof(ytmp) <: Vector{String}
    @test typeof(vars) <: Vector{String}
    @test size(xtmp) == (4, 150)
    @test length(ytmp) == 150
    @test length(vars) == size(xtmp, 1)
    @test mean(xtmp, dims=2) ≈ [5.843333333333333333, 3.05733333333333333, 3.758000, 1.199333333333333]
    @test mean(xtmp[:,1:50], dims=2) ≈ [5.006, 3.428, 1.462, 0.246]
    @test mean(xtmp[:,51:100], dims=2) ≈ [5.936, 2.77, 4.26, 1.326]
end

@testset "Test load_sin" begin
    xtmp, ytmp, vars = load_sin()

    @test typeof(xtmp) <: Vector{Float64}
    @test typeof(ytmp) <: Vector{Float64}
    @test typeof(vars) <: Vector{String}
    @test length(xtmp) == length(ytmp) == 15
    @test length(vars) == 2
end

@testset "Test load_line" begin
    xtmp, ytmp, vars = load_line()

    @test typeof(xtmp) <: Vector{Float64}
    @test typeof(ytmp) <: Vector{Float64}
    @test typeof(vars) <: Vector{String}
    @test length(xtmp) == length(ytmp) == 11
    @test length(vars) == 2
end

@testset "Test load_poly" begin
    xtmp, ytmp, vars = load_poly()

    @test typeof(xtmp) <: Vector{Float64}
    @test typeof(ytmp) <: Vector{Float64}
    @test typeof(vars) <: Vector{String}
    @test length(xtmp) == length(ytmp) == 50
    @test length(vars) == 2
end

@testset "Test load_spiral" begin
    xtmp, ytmp, vars = load_spiral()
    @test typeof(xtmp) <: Matrix{Float64}
    @test typeof(ytmp) <: Vector{Int}
    @test typeof(vars) <: Vector{String}
    @test length(xtmp[:, 1]) == length(ytmp) == 194
    @test length(vars) == 3
end


text:
# durbin_watson.jl
# Durbin-Watson test for autocorrelation
#
# Copyright (C) 2017   Benjamin Born
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export DurbinWatsonTest

struct DurbinWatsonTest <: HypothesisTest
    xmat::Array{Float64}  # regressor matrix
    n::Int                # number of observations
    DW::Float64           # test statistic
    p_compute::Symbol     # determines how p-values are computed
end

"""
    DurbinWatsonTest(X::AbstractArray, e::AbstractVector; p_compute::Symbol = :ndep)

Compute the Durbin-Watson test for serial correlation in the residuals of a regression model.

`X` is the matrix of regressors from the original regression model and `e` the vector of
residuals. Note that the Durbin-Watson test is not valid if `X` includes a lagged dependent
variable. The test statistic is computed as
```math
DW = \\frac{\\sum_{t=2}^n (e_t - e_{t-1})^2}{\\sum_{t=1}^n e_t^2}
```
where `n` is the number of observations.

By default, the choice of approach to compute p-values depends on the sample size
(`p_compute=:ndep`). For small samples (n<100), Pan's algorithm (Farebrother, 1980) is
employed. For larger samples, a normal approximation is used (Durbin and Watson, 1950). To
always use Pan's algorithm, set `p_compute=:exact`. `p_compute=:approx` will always use the
normal approximation.

Default is a two-sided p-value for the alternative hypothesis of positive or negative
serial correlation. One-sided p-values can be requested by calling
`pvalue(x::DurbinWatsonTest; tail=)` with the options `:left` (negative serial correlation)
and `:right` (positive serial correlation).

# References

  * J. Durbin and G. S. Watson, 1951, "Testing for Serial Correlation in Least Squares
    Regression: II", Biometrika, Vol. 38, No. 1/2, pp. 159-177,
    [http://www.jstor.org/stable/2332325](http://www.jstor.org/stable/2332325).
  * J. Durbin and G. S. Watson, 1950, "Testing for Serial Correlation in Least Squares
    Regression: I", Biometrika, Vol. 37, No. 3/4, pp. 409-428,
    [http://www.jstor.org/stable/2332391](http://www.jstor.org/stable/2332391).
  * R. W. Farebrother, 1980, "Algorithm AS 153: Pan's Procedure for the Tail Probabilities
    of the Durbin-Watson Statistic", Journal of the Royal Statistical Society, Series C
    (Applied Statistics), Vol. 29, No. 2, pp. 224-227,
    [http://www.jstor.org/stable/2986316](http://www.jstor.org/stable/2986316).

# External links

  * [Durbin-Watson test on Wikipedia](https://en.wikipedia.org/wiki/Durbin–Watson_statistic)
"""
function DurbinWatsonTest(xmat::AbstractArray{T}, e::AbstractArray{T};
    p_compute::Symbol = :ndep) where T<:Real

    n = length(e)
    DW = sum(diff(e) .^2) / sum(e .^2)

    DurbinWatsonTest(xmat, n, DW, p_compute)
end

testname(::DurbinWatsonTest) = "Durbin-Watson autocorrelation test"
population_param_of_interest(x::DurbinWatsonTest) =
    ("sample autocorrelation parameter", "0", 1 - x.DW / 2)
default_tail(test::DurbinWatsonTest) = :both

function show_params(io::IO, x::DurbinWatsonTest, ident)
    println(io, ident, "number of observations:     ", x.n)
    println(io, ident, "DW statistic:               ", x.DW)
end

"""
    pan_algorithm(a::AbstractArray, x::Float64, m::Int, n::Int)

Compute exact p-values for the Durbin-Watson statistic using Pan's algorithm (Farebrother,
1980).

`a` is the vector of non-zero Eigenvalues of ``(I-(X(X'X)^{-1}X'))A`` (see Durbin and
Watson, 1971, p. 2), `x` is the value of the Durbin-Watson statistic, `m` the number of
elements in `a`, and `n` the number of approximation terms (see Farebrother, 1980, eq. 5).

# References

  * J. Durbin and G. S. Watson, 1971, "Testing for Serial Correlation in Least Squares
  Regression: III", Biometrika, Vol. 58, No. 1, pp. 1-19,
  [http://www.jstor.org/stable/2334313](http://www.jstor.org/stable/2334313).
  * R. W. Farebrother, 1980, "Algorithm AS 153: Pan's Procedure for the Tail Probabilities
  of the Durbin-Watson Statistic", Journal of the Royal Statistical Society, Series C
  (Applied Statistics), Vol. 29, No. 2, pp. 224-227,
  [http://www.jstor.org/stable/2986316](http://www.jstor.org/stable/2986316).

"""
function pan_algorithm(a::AbstractArray, x::Float64, m::Int, n::Int)

    ν = findfirst(ai -> ai >= x, a)
    if ν == 0 || ν === nothing
        return 1.0
    elseif ν == 1
        return 0.0
    else
        k = 1
        ν = ν - 1
        h = m - ν

        if ν <= h
            d  = 2; h  = ν; k  = - k; j1 = 0; j2 = 2; j3 = 3; j4 = 1
        else
            d  = - 2; ν  = ν + 2; j1 = m - 2; j2 = m - 1; j3 = m + 1; j4 = m
        end

        pin = pi / (2n)
        sum = (k + 1) / 2
        sgn = k / n
        n2  = 2n - 1

        # first integral
        for  f1 = h - 2 * floor(Int,h/2) : -1 : 0
            for f2 = j2:d:ν
                sum1 = a[j4]
                if f2 == 0
                    prod = x
                else
                    prod = a[f2]
                end
                u = 0.5 * (sum1 + prod)
                v = 0.5 * (sum1 - prod)
                sum1 = 0.0
                for  i = n2:-2:1
                    y = u - v * cos(i * pin)
                    num = y - x
                    prod = 1.0
                    for k in [(1:j1)' (j3:m)']
                        prod *= num / (y - a[k])
                    end
                    sum1 += sqrt(abs(prod))
                end
                sgn = -sgn
                sum += sgn * sum1
                j1 += d
                j3 += d
                j4 += d
            end

            # second integral
            if d == 2
                j3 = j2
            else
                j1 = j2
            end
            j2 = 0
            ν  = 0
        end
        return sum
    end

end

function pvalue(x::DurbinWatsonTest; tail=:both)

    exact_problem_flag = 0
    if (x.p_compute == :ndep && x.n <= 100) || x.p_compute == :exact
        # p-vales based on Pan's algorithm (see Farebrother, 1980)

        # the following setup is, e.g, described in Durbin and Watson (1971)
        A = diagm(-1 => -ones(x.n - 1), 0 => 2 * ones(x.n), 1 => -ones(x.n - 1))
        A[1, 1] = 1
        A[x.n, x.n] = 1
        EV_temp = sort(real(eigvals(A - x.xmat / (x.xmat' * x.xmat) * (x.xmat'*A))))
        EV = EV_temp[EV_temp .> 1e-10]
        p_temp = pan_algorithm(EV, x.DW, length(EV), 15)

        if p_temp < 0.0 || p_temp > 1.0
            # println(p_temp)
            @warn("Exact p-values outside [0,1]. Approximate p-values reported instead.")
            exact_problem_flag = 1
        end
    end

    if exact_problem_flag == 1 || (x.p_compute == :ndep && x.n > 100
        ) || x.p_compute == :approx
        # p-values based on normal approximation (see Durbin and Watson, 1950)

        # the following derivations follow Durbin and Watson (1951, p. 164)
        X = x.xmat
        k = size(X, 2)
        inv_XX = inv(X' * X)

        AX = zeros(x.n, k)
        AX[[1, x.n], :] = X[[1, x.n], :] - X[[2, x.n - 1], :]
        for i = 2:(x.n - 1)
            AX[i, :] = - X[i - 1, :] + 2 * X[i, :] - X[i + 1, :]
        end

        temp_mat = X' * AX * inv_XX
        P = 2 * (x.n - 1) - tr(temp_mat) # first term: tr(A)
        Q = 2 * (3 * x.n - 4) - 2 * tr(AX' * AX * inv_XX) + tr(temp_mat^2)
        # first term: tr(A^2)

        dw_mean = P / (x.n - k)
        dw_var = 2 / ((x.n - k) * (x.n - k + 2)) * (Q - P * dw_mean)

        p_temp = cdf(Normal(dw_mean, sqrt(dw_var)), x.DW)
    end

    if tail == :both
        2 * min(p_temp, 1 - p_temp)
    elseif tail == :right
        p_temp
    elseif tail == :left
        1 - p_temp
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end

end


text:
using GLVisualize, GLAbstraction, Reactive, GeometryTypes, Colors, GLWindow

if !isdefined(:runtests)
    window = glscreen()
    addprocs(1)
end

description = """
Example showing off how to run computations in a different process and then
visualizing them.
"""

const workerid = workers()[]

@everywhere begin
    using GeometryTypes
    function spin(dir, α)
        x, y = dir
        Point(x*cos(α) - y*sin(α), x*sin(α) + y*cos(α))
    end
    function fractal_step!(
            a, b, depth, angles,
            result = Point2f0[a], levels = Float32[depth] # tmp too not allocate
        )
        depth == 0 && return result, levels, b
        N = length(angles)
        diff = (b - a)
        len = norm(diff) / (N - 1)
        nth_segment = normalize(diff)
        for n = 1:N
            b = a + spin(nth_segment * len, angles[n]) # rotate with current angle
            _, _, b = fractal_step!(a, b, depth-1, angles, result, levels) # recursion step
            if n < N
                push!(result, b)
                push!(levels, depth)
            end
            nth_segment = normalize(b - a)
            a = b
        end
        result, levels, b
    end
    function generate_fractal(angles, depth = 5)
        tmp = zeros(Point2f0, length(angles))
        angles = map(deg2rad, angles)
        result, levels, b = fractal_step!(Point2f0(0,0), Point2f0(300,0), round(Int, depth), angles)
        push!(result, b)
        push!(levels, depth)
        result, depth .- levels
    end
end

editarea, viewarea = x_partition(window.area, 30.0)
editscreen = Screen(
    window, area = editarea,
    color = RGBA{Float32}(0.1f0, 0.1f0, 0.1f0, 1f0),
    stroke = (1f0, RGBA{Float32}(0.13f0, 0.13f0, 0.13f0, 13f0))
)
viewscreen = Screen(window, area=viewarea, color = RGBA(0.1f0, 0.1f0, 0.1f0, 1f0))

controls = [
    :angles => Vec4f0(0.0, 80.0, -140.0, 80.0),
    :colormap => map(RGBA{Float32}, colormap("Blues", 10)),
    :thickness => 0.5f0,
    :iterations => 6,
]

menu, controls = GLVisualize.extract_edit_menu(controls, editscreen, true)
_view(menu, editscreen, camera = :fixed_pixel)

const channel2 = Channel{Tuple{Vector{Point2f0}, Vector{Float32}}}(1)
v0 = (Point2f0[0,0,0,0], Float32[0,0,0,0])
put!(channel2, v0)
line_level = foldp(v0, controls[:angles], controls[:iterations]) do v0, angles, iter
    if isready(channel2)
        points = take!(channel2)
        task = @async put!(
            channel2,
            remotecall_fetch(generate_fractal, workerid, angles, iter)
        )
        return points
    end
    v0
end

_view(visualize(
    map(first, line_level), :lines,
    intensity = map(last, line_level),
    thickness = controls[:thickness],
    color_map = controls[:colormap],
    color_norm = Vec2f0(0,10),
), viewscreen, camera=:orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


text:
################################################################################
# CircularWall
################################################################################

mutable struct CircularWall{T}
    x::T
    y::T
    r::T
end

function CircularWall()
    x = 0.0
    y = 0.0
    r = 0.0
    return CircularWall{typeof(x)}(x,y,r)
end

################################################################################
# Lane
################################################################################

mutable struct Lane{T}
    id::Int
    name::Symbol
    wall::Vector{Wall}
    circle::Vector{CircularWall{T}}
    start::StartingArea{T}
    centerline::Function
end

function Lane()
    name = :lane_0
    w = Wall([0., 0.], [0., 0.], [0., 0.])
    c = CircularWall()
    start = StartingArea()
    return Lane{typeof(c.x)}(0, name, [w], [c], start, x -> x)
end

function Lane(id::Int, name::Symbol, wall::Vector{Wall}, circle::Vector{CircularWall}, start::StartingArea{T}, centerline) where {T}
    return Lane{T}(id, name, wall, circle, start, centerline)
end


text:
function fscore(golds::Vector{UnitRange{Int}}, preds::Vector{UnitRange{Int}})
    g = Set(golds)
    p = Set(preds)
    inter = intersect(g, p)
    correct = length(inter)

    prec = correct / length(p)
    recall = correct / length(g)
    fvalue = 2*recall*prec/(recall+prec)
    println("precision: $(prec)")
    println("recall: $(recall)")
    println("f-value: $(fvalue)")
end

function accuracy(golds::Vector{Int}, preds::Vector{Int})
    @assert length(golds) == length(preds)
    correct = 0
    total = 0
    for i = 1:length(golds)
        golds[i] == preds[i] && (correct += 1)
        total += 1
    end
    correct / total
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using Rhapsodie
using DelimitedFiles
using EasyFITS
using InterpolationKernels

#include("test_data_simulation.jl")
#-----------------------------------------------
# Loading the dataset parameters
par=readdlm("test_results/Parameters.txt")

DSIZE=Int64(par[1]);
NTOT=Int64(par[2]);
Nframe=Int64(par[3]);
Nrot=Int64(par[4]);
Nangle=NTOT÷(Nframe*4)
Center=par[5:6];
DerotAng=Vector{Float64}();
Epsilon=Vector{Tuple{Array{Float64,1},Array{Float64,1}}}();

for iter=1:NTOT
    ind=div(iter-1, NTOT/4)+1
    push!(Epsilon,([0. ,0. ],par[end-1:end]));
    push!(DerotAng, deg2rad(-par[7]));
end

psf_center=readdlm("../data/PSF_centers_Airy.txt");

load_parameters((DSIZE, 2*DSIZE, NTOT), Nframe, Nrot, Nangle, Center, (psf_center[1:2], psf_center[3:4]), Epsilon, DerotAng)

for tau in [0.03, 0.03]#, 0.07, 0.1, 0.15, 0.25, 0.5]
#----------------------------------------------
# Loading the pre-processed data (i.e. bad pixels have been interpolated
    Load_Data("test_results/DATA_$tau-$DSIZE.fits",
                         "test_results/WEIGHT_$tau-$DSIZE.fits")
	Sdim=length(Rhapsodie.dataset)

    DATA=zeros(get_par().cols[1], get_par().cols[2], Sdim,2);
    WEIGHT=zeros(get_par().cols[1], get_par().cols[2], Sdim,2);

    ker= CatmullRomSpline(Float64, Flat)
    input_size=(get_par().rows[1], get_par().rows[2]÷2);
    output_size= get_par().cols[1:2];

# Pre processing
    for i=1:Sdim
        T1=TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, inv(Rhapsodie.Trans_Table[i][1]))
        T2=TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, inv(Rhapsodie.Trans_Table[i][2]))

        I1=T1*Rhapsodie.dataset[i].data[:,1:end÷2]
        I2=T2*Rhapsodie.dataset[i].data[:,end÷2+1:end]
        W1=T1*Rhapsodie.dataset[i].weights[:,1:end÷2]
        W2=T2*Rhapsodie.dataset[i].weights[:,end÷2+1:end]

        DATA[:,:,i,1]=I1;
        DATA[:,:,i,2]=I2;
        WEIGHT[:,:,i,1]=W1;
        WEIGHT[:,:,i,2]=W2;
    end

#-----------------------------------------------------
# Reconstruction
# Non-Linear Separable inverse method
    #MI=NonLinear_Method(DATA, WEIGHT);
    #write(MI, "test_results/Results_Separable_NonLinear_$tau-$DSIZE.fits")

# Linear Separable inverse method
    ML=Linear_Method(DATA, WEIGHT);
    crop!(ML)
    write(crop(ML), "test_results/Results_Separable_Linear_$tau-$DSIZE.fits")

# Double Ratio
    DR=Double_Ratio(DATA);
    write(DR, "test_results/Results_Separable_DoubleRatio_$tau-$DSIZE.fits")

# Double Difference
    DD=Double_Difference(DATA);
    write(DD, "test_results/Results_Separable_DoubleDifference_$tau-$DSIZE.fits")

    empty!(Rhapsodie.dataset);
end

max_stars_repo_name: LaurenceDenneulin/Rhapsodie

max_issues_repo_name: LaurenceDenneulin/Rhapsodie

max_forks_repo_name: LaurenceDenneulin/Rhapsodie


text:
function sortstring(str, ordering)
    order = Dict( map( reverse, enumerate(unique(ordering)) ) )
    N = length(order)+1
    specialorder = c -> haskey(order,c) ? order[c] : N
    join( sort([c for c in str], by=specialorder ) )
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using JTetra
using Base.Test

println("==== Test De Rham Complex ====")

mesh = getTetraMesh("points.txt", "tetras.txt")

G = getNodalGradientMatrix(mesh)
C = getCurlMatrix(mesh)
D = getDivergenceMatrix(mesh)

# @printf("||G||     = %.3e\n", norm(G,1))
# @printf("||C||     = %.3e\n", norm(C,1))
# @printf("||D||     = %.3e\n", norm(D,1))
@printf("||C * G|| = %.3e\n", norm(C*G,1))
@printf("||D * C|| = %.3e\n", norm(D*C,1))

@test norm(C*G,1) < 1e-12
@test norm(D*C,1) < 1e-12

print("Passed!\n")

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type flash
	flash()=begin
		PP=outers.PP
		NComp=outers.NComp
		[
			:(LI = Levelpercent_Initial),
			:(OutletLiquid.T = Temperature_Initial),
			:(OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial)),
		],
		[
			"Initial level Percent","Initial Outlet Liquid Temperature","Initial Outlet Liquid Composition Normalized",
		],
		new(
			VesselVolume(),
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components",
				:Lower=>1
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Mol Weight",
				:Hidden=>true
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Initial liquid height in Percent",
				:Default=>0.70
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Initial Liquid Temperature",
				:Default=>330
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Initial Composition",
				:Default=>0.10
			)),(NComp)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed Stream",
				:PosX=>0,
				:PosY=>0.48,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid outlet stream",
				:PosX=>0.43,
				:PosY=>1,
				:Symbol=>"_{out}^{Liquid}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0.43,
				:PosY=>0,
				:Symbol=>"_{out}^{Vapour}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat Duty",
				:PosX=>1,
				:PosY=>0.81,
				:Protected=>true,
				:Symbol=>"Q_{in}"
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the Vessel",
				:Protected=>true
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup",
				:Protected=>true
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup",
				:Protected=>true
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy Holdup in the Vessel",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume",
				:Protected=>true
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Vapourization fraction",
				:Symbol=>"\\phi",
				:Protected=>true
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Pressure Ratio",
				:Symbol=>"P_{ratio}",
				:Protected=>true
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop",
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P",
				:Protected=>true
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Equilibrium pressure on the liquid surface",
				:Protected=>true,
				:Symbol=>"\\Delta P_{eq}"
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Static head at the bottom of the tank",
				:Protected=>true,
				:Symbol=>"P_{static}^{Liquid}"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature Indicator",
				:PosX=>1,
				:PosY=>0.39,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure Indicator",
				:PosX=>1,
				:PosY=>0.21,
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Level Indicator",
				:PosX=>1,
				:PosY=>0.59,
				:Protected=>true
			)),
			[
				:(diff(TotalHoldup)=Inlet.F*Inlet.z - OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z),
				:(diff(E) = Inlet.F*Inlet.h - OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ),
				:(TotalHoldup = LiquidHoldup*OutletLiquid.z + VapourHoldup*OutletVapour.z),
				:(E = LiquidHoldup*OutletLiquid.h + VapourHoldup*OutletVapour.h - OutletLiquid.P*Vtotal),
				:(sum(OutletLiquid.z)=1.0),
				:(sum(OutletLiquid.z)=sum(OutletVapour.z)),
				:(OutletVapour.F = Inlet.F * vfrac),
				:(vL = PP.LiquidVolume(OutletLiquid.T, Peq, OutletLiquid.z)),
				:(vV = PP.VapourVolume(OutletVapour.T, Peq, OutletVapour.z)),
				:(PP.LiquidFugacityCoefficient(OutletLiquid.T, Peq, OutletLiquid.z)*OutletLiquid.z = PP.VapourFugacityCoefficient(OutletVapour.T, Peq, OutletVapour.z)*OutletVapour.z),
				:(OutletVapour.T = OutletLiquid.T),
				:(OutletVapour.P = Peq),
				:(Pstatic = PP.LiquidDensity(OutletLiquid.T, Peq, OutletLiquid.z) * _base_1.Gconst * Level),
				:(OutletLiquid.P = Peq + Pstatic),
				:(OutletLiquid.P = Inlet.P - Pdrop),
				:(OutletLiquid.P = Inlet.P * Pratio),
				:(Vtotal = LiquidHoldup * vL + VapourHoldup * vV),
				:(TI * "K" = OutletLiquid.T),
				:(PI * "atm" = OutletLiquid.P),
				:(LI*Vtotal= Vfilled),
				:(LiquidHoldup * vL = Vfilled),
			],
			[
				"Component Molar Balance","Energy Balance","Molar Holdup","Energy Holdup","Mol fraction normalisation","Mol fraction normalisation","Vaporization Fraction","Liquid Volume","Vapour Volume","Chemical Equilibrium","Thermal Equilibrium","Mechanical Equilibrium for the Vapour Phase","Static Head","Mechanical Equilibrium for the Liquid Phase","Pressure Drop","Pressure Ratio","Geometry Constraint","Temperature indicator","Pressure indicator","Level indicator","Liquid Level",
			],
			[:PP,:NComp,:Mw,:Levelpercent_Initial,:Temperature_Initial,:Composition_Initial,],
			[:Inlet,:OutletLiquid,:OutletVapour,:InletQ,:TotalHoldup,:LiquidHoldup,:VapourHoldup,:E,:vL,:vV,:vfrac,:Pratio,:Pdrop,:Peq,:Pstatic,:TI,:PI,:LI,]
		)
	end
	_base_1::VesselVolume
	PP::DanaPlugin
	NComp::DanaInteger
	Mw::Array{molweight}
	Levelpercent_Initial::positive
	Temperature_Initial::temperature
	Composition_Initial::Array{fraction}
	Inlet::stream
	OutletLiquid::liquid_stream
	OutletVapour::vapour_stream
	InletQ::power
	TotalHoldup::Array{mol}
	LiquidHoldup::mol
	VapourHoldup::mol
	E::energy
	vL::volume_mol
	vV::volume_mol
	vfrac::positive
	Pratio::positive
	Pdrop::press_delta
	Peq::pressure
	Pstatic::pressure
	TI::control_signal
	PI::control_signal
	LI::control_signal
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export flash
function set(in::flash)
	Mw=PP.MolecularWeight()

end
function initial(in::flash)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function setEquationFlow(in::flash)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
end
function atributes(in::flash,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Flash"
	fields[:Brief]="Model of a Dynamic Flash Vessel."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;

== SPECIFY ==
* the Inlet stream;
* the outlet flows: OutletVapour.F and OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the Flash model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the Flash temperature (OutletLiquid.T);
* Initial_Level : the Flash liquid level (Level);
* Initial_Composition : (NoComps) OutletLiquid compositions.
"
	drive!(fields,_)
	return fields
end
flash(_::Dict{Symbol,Any})=begin
	newModel=flash()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
struct TaskManager
    queue::Channel{Task}
    task::Base.RefValue{Task}
end

TaskManager() = TaskManager(Channel{Task}(Inf), Ref{Task}())

function launch!(m::TaskManager)
    m.task[] = @async for (i, t) in enumerate(m.queue)
        try
            schedule(t)
            wait(t)
        catch exception
            @error "Error from $i-th task" exception
        end
    end
end

get_ready!(m::TaskManager) = isassigned(m.task) || launch!(m)

function enqueue!(m::TaskManager, t::Task)
    get_ready!(m)
    push!(m.queue, t)
    return
end

const taskmanager = Ref{TaskManager}()

function _taskmanager()
    if !isassigned(taskmanager)
        taskmanager[] = TaskManager()
    end
    return taskmanager[]
end

enqueue(t::Task) = enqueue!(_taskmanager(), t)


text:
module NaiveDataset

import ILKBackend

export BinDataset, openBinDataset, readVector, readPose

mutable struct BinDataset

    istream
    eof::Bool
    buf
end

function openBinDataset(filename::String)
        i = open(filename, "r")
        buf = Array{Float32, 1}(undef,128);
        return BinDataset(i, false, buf)
end

function readVector(ds::BinDataset, size::UInt16)
    readBuff(ds, size)
    buf = Array{Float32, 1}(undef, size)
    buf[:] = ds.buf[1:size]
    return buf
end

function readPose(ds::BinDataset)
    if ! readBuff(ds, UInt16(12))
        return None
    end
    pose = ILKBackend.pose()
    R = ILKBackend.rotationView(pose);
    p = ILKBackend.positionView(pose);
    p[1] = ds.buf[1];
    p[2] = ds.buf[2];
    p[3] = ds.buf[3];

    # expect the rotation matrix elements in row-major order, in the buffer
    for r in 1:3
        for c in 1:3
            R[r,c] = ds.buf[r*3 + c];
        end
    end
    return pose
end

sizeOfFloatInDataset = sizeof(Float32)
bytes = Array{UInt8, 1}(undef, 64)

function readBuff(dataset::BinDataset, howmany::UInt16)
    if howmany == 0
        return
    end

    bytesCount  = howmany*sizeOfFloatInDataset
    actualCount = readbytes!(dataset.istream, bytes, bytesCount )

    if actualCount==0
        dataset.eof = true
        return False
    end
    if actualCount<bytesCount
        print("Warning, unexpected number of bytes read. Assuming EOF")
        dataset.eof = true
        return False
    end

    # Decode the raw bytes into floats
    for f in 0:howmany-1
        buffer = IOBuffer(bytes[f*sizeOfFloatInDataset+1:(f+1)*sizeOfFloatInDataset+1])
        dataset.buf[f+1] = read(buffer, Float32)
    end
    dataset.eof = eof(dataset.istream)

    return true
end

end

max_stars_repo_name: ESROCOS/ilk-compiler

max_issues_repo_name: ESROCOS/ilk-compiler

max_forks_repo_name: ESROCOS/ilk-compiler


text:
# Autogenerated wrapper script for Libcerf_jll for x86_64-w64-mingw32
export libcerf

JLLWrappers.@generate_wrapper_header("Libcerf")
JLLWrappers.@declare_library_product(libcerf, "libcerf.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcerf,
        "bin\\libcerf.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "228.summary-ranges.jl" begin
    @test summary_ranges([0,2,3,4,6,8,9]) == ["0","2->4","6","8->9"]
    @test summary_ranges([0,1,2,4,5,7]) == ["0->2","4->5","7"]
    @test summary_ranges(Int[]) == String[]
    @test summary_ranges([-1]) == ["-1"]
end


text:
function composer(fns::Vector{GradientHook})
    return if length( fns ) == 0
        identity
    else
        reduce( ∘ , reverse( fns ) )
    end
end
composer(fns::Nothing) = identity

struct IdentityHook <: GradientHook; end
(a::IdentityHook)( x ) = x
@adjoint (a::IdentityHook)( x ) = x, Δ -> Δ
Base.show(io::IO, l::IdentityHook) = print(io, "Identity Gradient Transform")

struct L1Hook{F} <: GradientHook
    λ::F
end
ΔL1(Δ, θ, λ) = ( Δ .+ ( λ .* sign.( Δ ) ) ) .* ( abs.( θ ) .> λ )
( a::L1Hook )( x ) = x
@adjoint ( a::L1Hook )( θ ) = θ, Δ -> ( nothing, ΔL1( Δ, θ, a.λ ) )
Base.show(io::IO, l::L1Hook) = print(io, "L1 Penalty (λ = $(l.λ))")

struct L2Hook{F} <: GradientHook
    λ::F
end
ΔL2(Δ, θ, λ) = Δ .+ ( λ .* sum(abs2, θ ) )
( a::L2Hook )( x ) = x
@adjoint ( a::L2Hook )( θ ) = θ, Δ -> ( nothing, ΔL2( Δ, θ, a.λ ) )
Base.show(io::IO, l::L2Hook) = print(io, "L2 Penalty (λ = $(l.λ))")

struct LipschitzHook{F} <: GradientHook
    λ::F
end
LipschitzRegularization(Δ, θ, λ) = Δ .+ ( λ .* sum(abs2, θ ) )
( a::LipschitzHook )( x ) = x
@adjoint ( a::LipschitzHook )( θ ) = θ, Δ -> ( nothing, LipschitzRegularization( Δ, θ, a.λ ) )
Base.show(io::IO, l::LipschitzHook) = print(io, "Lipschitz Penalty (λ = $(l.λ))")

struct ReverseGradientHook <: GradientHook; end
( a::ReverseGradientHook )( x ) = x
@adjoint ( a::ReverseGradientHook )( θ ) = θ, Δ -> (nothing, -Δ)
Base.show(io::IO, l::ReverseGradientHook) = print(io, "Reverse Gradient")

struct ClampHook{F} <: GradientHook
    ϵ::F
end
( a::ClampHook )( x ) = x
@adjoint ( a::ClampHook )( θ ) = θ, Δ -> ( nothing, clamp.( Δ, -a.ϵ, a.ϵ ) )
Base.show(io::IO, a::ClampHook) = print(io, "Gradient Clamp (Δ ∈ [-$(a.ϵ),$(a.ϵ)])")


text:
using BenchmarkTools, FinancialMonteCarlo, ReverseDiff

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2
p = 0.3;
lam = 5.0;
lamp = 30.0;
lamm = 20.0;
mc = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = KouProcess(sigma, lam, p, lamp, lamm, Underlying(S0, d));
Model = KouProcess(x[1], lam, p, lamp, lamm, Underlying(S0, d));

#f(x) = pricer(KouProcess(x[1],lam,p,lamp,lamm,Underlying(S0,d)),ZeroRate(x[2],x[3],x[4]),mc,EuropeanOption(x[5],K));
f(x) = pricer(KouProcess(x[1], lam, p, lamp, lamm, Underlying(S0, d)), ZeroRate(x[2], r, d), mc, EuropeanOption(T, K));
#x=Float64[sigma,S0,r,d,T]
x = Float64[sigma, S0]
g = x -> ReverseDiff.gradient(f, x);
y0 = g(x);
@btime f(x);
@btime g(x);

#f1_(x) = pricer(BlackScholesProcess(x[1],Underlying(S0,d)),ZeroRate(x[2],x[3],x[4]),mc,AmericanOption(x[5],K));
f1_(x) = pricer(KouProcess(x[1], lam, p, lamp, lamm, Underlying(S0, d)), ZeroRate(x[2], r, d), mc, AmericanOption(T, K));
g_ = x -> ReverseDiff.gradient(f1_, x);
#y0=g(x);
@btime f1_(x);
#@btime g_(x);


text:
############################################################################################
"""
$(TYPEDEF)

Sampling Container that determines tuning updates.

# Fields
$(TYPEDFIELDS)
"""
abstract type SamplingPhase end
" Find/converges to the typical set (estimators suffer from initial but ultimately transient biases)"
struct Warmup <: SamplingPhase end
"Explore typical set and tune MCMC parameter(initial bias rapidly vanishes and the estimators become much more accurate.)"
struct Adaptionˢˡᵒʷ <: SamplingPhase end
"Tune final stepsize for fixed Covariance Matrix."
struct Adaptionᶠᵃˢᵗ <: SamplingPhase end
"Sample parameter with tuned MCMC (gradually reducing the precision error of the MCMC estimators towards zero.)"
struct Exploration <: SamplingPhase end

############################################################################################
"""
$(TYPEDEF)

Default configuration for tuning iterations in each cycle.

# Fields
$(TYPEDFIELDS)
"""
struct ConfigTuningWindow
    "MCMC Phase tuning window lengths, i.e: 5 = 5 repeats of second window in phasenames"
    window::Vector{Int64}
    "(Increasing) Length of windows, i.e.: if window=3, buffer=10 -> total window length: 10-20-30"
    buffer::Vector{Int64}
    "Name of phasenames, currently supported: Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(), Exploration()"
    phasenames::Vector{SamplingPhase}
    function ConfigTuningWindow(
        window::Vector{Int64},
        buffer::Vector{Int64},
        phasenames::Vector{SamplingPhase}
        )
        @argcheck size(window, 1) == size(buffer, 1) "Window and Buffer size different"
        @argcheck (size(window, 1) + 1) == size(phasenames, 1) "Window and phasename (without Exploration) size different"
        return new(
            window, buffer, phasenames
            )
    end
end

############################################################################################
"""
$(TYPEDEF)

Information about current SamplingPhase.

# Fields
$(TYPEDFIELDS)
"""
struct PhaseTune{T<:Tuple}
    "Boolean if current iteration needs update."
    update::Updater
    "Counts current iteration in phase."
    iter::Iterator
    "MCMC Phases ~ counter = current cyle in slices/name/iterations."
    counter::Iterator
    "Vector of MCMC iterations at each phase."
    slices::Vector{Int64}
    "Name of Sampling phases."
    name::T
    "Counts total iterations."
    iterations::Vector{Int64}
    ## MCMC Windows
    "Counts cycles in adaption phases Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(),
    i.e., 1-5-1 means 1 window init, 5 adapt, 1 exploration."
    window::Vector{Int64}
    "Counts iteration in each cycle for adaption phases Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(),
    i.e., 50-25-50 means 1 time 50, 5 times 25*i, i in 1:5, 1 time 50."
    buffer::Vector{Int64}
    function PhaseTune(;
        window=[1, 5, 1],
        buffer=[50, 25, 50],
        phasenames=[Warmup(), Adaptionˢˡᵒʷ(), Adaptionᶠᵃˢᵗ(), Exploration()],
    )
        ArgCheck.@argcheck length(window) == length(buffer) == (length(phasenames) - 1) "Size of window, buffer, or phasenames are not matching."
        ## Initiate mutable update container
        update = Updater(false)
        #!NOTE: start with 0, so first propose!() call will update iter.current to 1, and proposal will be first filled with 1st column.
        iter = Iterator(0)
        #!NOTE: start with 1 as Warmup matrix is initiated before propose! steps
        counter = Iterator(1)
        ## Initiate Warmup and Adaption Cycles
        phase_name = SamplingPhase[]
        phase_iterations = Int64[]
        for i in eachindex(window)
            for j in Base.OneTo(window[i])
                push!(phase_iterations, convert(Int64, j * buffer[i]))
                push!(phase_name, phasenames[i])
            end
        end
        ## Include Exploration Cycle
        push!(phase_name, phasenames[end])
        push!(phase_iterations, 1)
        ## Calculate cumsum of cylces for chainₜ iterations
        phase_slices = cumsum(phase_iterations)
        ## Change type of phase_name for concrete struct type
        phase_name_tuple = Tuple(phase_name)
        return new{typeof(phase_name_tuple)}(
            update,
            iter,
            counter,
            phase_slices,
            phase_name_tuple,
            phase_iterations,
            window,
            buffer,
        )
    end
end

############################################################################################
"""
$(SIGNATURES)
Update current MCMC Phase.

# Examples
```julia
```

"""
function update!(phase::PhaseTune, iter::Int64)
    @unpack slices = phase
    ## Check if updated last iteration
    if phase.update.current
        ## If update == true, set updating to false if updated last iteration
        update!(phase.update)
        ## Set back to 0 and add 1 afterwards
        phase.iter.current = 0
    end
    ## If update now, change Sampling phase
    if phase.counter.current < length(slices) && slices[phase.counter.current] == iter
        update!(phase.update)
        update!(phase.counter)
    end
    ## Update current iteration
    update!(phase.iter)
    return nothing
end

############################################################################################
# Export
export SamplingPhase, ConfigTuningWindow, Warmup, Adaptionˢˡᵒʷ, Adaptionᶠᵃˢᵗ, Exploration, PhaseTune


text:
HyperBox(a::Array{T,2}) where T = HyperBox{T}([Interval(a[1,i],a[2,i]) for i=1:size(a,2)])

## Access
## ======
getindex(b::HyperBox, i::Int) = b.intervals[i]
setindex!(b::HyperBox{T}, i::Int, v::Interval{T}) where T = b.intervals[i] = v
ndims(b::HyperBox) = length(b.intervals)
dims(b::HyperBox) = collect(1:ndims(b)) # make this iterable?
hasdim(b::HyperBox, i::Int) = i <= ndims(b)

## Splitting
## =========
mid(b::HyperBox{T}) where T = T[mid(b[dim]) for dim in dims(b)]

"""
    split_box(b::HyperBox{T}, split_point::Vector{T}) where T

Split box into 2^d equally sized boxes by cutting down middle of each axis.
"""
function split_box(b::HyperBox{T}, split_point::Vector{T}) where T
  @assert(length(split_point) == ndims(b))
  splits = [split_box(b[i],split_point[i]) for i = 1:ndims(b)] # Split intervals
  intervals_set = prodsubboxes(splits)
  HyperBox[HyperBox([intervals...]) for intervals in intervals_set]
end

"""
    partial_split_box(b::HyperBox{T}, split_point::Dict{Int, T}) where T

Split a box along not all dimensions. `split_point` maps integer dimensions to point in that dimension to split.
"""
function partial_split_box(b::HyperBox{T}, split_point::Dict{Int, T}) where T
  @assert(length(keys(split_point)) <= ndims(b))
  intervals_set = product_boxes(b, split_point)
  HyperBox[HyperBox([intervals...]) for intervals in intervals_set]
end

## Print
## =====
string(b::HyperBox) = b.intervals
print(io::IO, b::HyperBox) = print(io, string(b))
show(io::IO, b::HyperBox) = print(io, string(b))
showcompact(io::IO, b::HyperBox) = print(io, string(b))