text:
using AbstractStorage, FolderStorage, JSON, Serialization, Test

base = get(ENV, "FOLDERSTORAGE_TESTDIR", ".")
@info "running tests in $base"

@testset "mkpath" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    @test isdir("foo")
    rm(c)
end

@testset "rm folder" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    rm(c)
    @test !isdir(joinpath(base,"foo"))
end

@testset "rm file" begin
    c = Folder(joinpath(base,"foo"))
    write(c, "bar", rand(2))
    @test isfile(joinpath(base,"foo","bar"))
    rm(c, "bar")
    @test !isfile(joinpath(base,"foo","bar"))
    @test isdir(joinpath(base,"foo"))
    rm(c)
end

@testset "copy" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    write(c, "o", rand(10))
    d = copy(c)
    @test !isdir(d)
    mkpath(d)
    @test isempty(readdir(d))
    rm(c)
    rm(d)
end

@testset "deepcopy" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    x = rand(10)
    write(c, "o", x)
    d = deepcopy(c)
    @test isdir(d)
    @test read!(d, "o", Vector{Float64}(undef, 10)) == x
    rm(c)
end

@testset "isfile" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    x = rand(10)
    write(c, "o", x)
    @test isfile(c, "o")
    @test !isfile(c, "j")
    rm(c)
end

@testset "filesize" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    x = rand(10)
    write(c, "o", x)
    @test filesize(c, "o") == 80
    rm(c)
end

@testset "write, canonical" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    x = rand(10)
    write(c, "o", x)
    @test read!(c.foldername*"/o", Vector{Float64}(undef, 10)) ≈ x
    rm(c)
end

@testset "read!, canonical" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    x = rand(10)
    write(c.foldername*"/o", x)
    @test read!(c, "o", Vector{Float64}(undef, 10)) == x
    rm(c)
end

struct Foo
    x::Int
    y::Float64
end
@testset "write, serialized" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    x = [Foo(1,2.0)]
    write(c, "o", x)
    io = open(c.foldername*"/o")
    _x = deserialize(io)
    close(io)
    @test _x[1].x == x[1].x
    @test _x[1].y ≈ x[1].y
    rm(c)
end

@testset "read!, serialized" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    io = open(c.foldername*"/o", "w")
    x = [Foo(1,2.0)]
    serialize(io, x)
    close(io)
    _x = read!(c, "o", Vector{Foo}(undef, 1))
    @test _x[1].x == x[1].x
    @test _x[1].y ≈ x[1].y
end

@testset "read/write string" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    write(c, "bar", "hello")
    @test read(c, "bar", String) == "hello"
    rm(c)
end

@testset "implicit pathing" begin
    c = Folder(joinpath(base,"foo"))
    mkpath(c)
    write(c, "bar/baz", "hello")
    @test read(joinpath(c.foldername,"bar/baz"), String) == "hello"
    rm(c)
end

@testset "json" begin
    c = Container(Folder, JSON.parse(json(Folder(joinpath(base,"foo"), nretry=5))))
    if isabspath(normpath(joinpath(base,"foo")))
        @test c.foldername == normpath(joinpath(base,"foo"))
    else
        @test c.foldername == normpath(joinpath(pwd(),base,"foo"))
    end
    @test c.nretry == 5
    rm(c)
end

@testset "minimal dictionary" begin
    c = Container(Folder, JSON.parse(json(Folder(joinpath(base,"foo"), nretry=5))))
    _c = minimaldict(c)
    @test _c["foldername"] == normpath(joinpath(pwd(),base,"foo"))
    @test length(_c) == 1
end


text:
module Pervasives
using MLStyle.MatchCore
using MLStyle.Infras
using MLStyle.Extension
using MLStyle.Err
using MLStyle.Toolz: ($)

import MLStyle.Infras: mk_gapp_pattern

@use GADT

function mk_pat_by(f)
    (tag, case, mod) -> body ->
     @format [f, tag, case, body] quote
        f(tag, case) ? body : failed
     end
end

const strict_eq_types = Union{Int, Nothing}

def_pattern(Pervasives,
        predicate = x -> x isa strict_eq_types,
        rewrite = mk_pat_by(===)
)

def_pattern(Pervasives,
        predicate = x -> x isa Union{Number, AbstractString, AbstractChar, QuoteNode},
        rewrite   = (tag, case, mod) ->
        let f = isimmutable(case) ? mk_pat_by(===) : mk_pat_by(isequal)
            f(tag, case, mod)
        end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head == :(||),
        rewrite = (tag, case, mod) ->
            let pats = [mk_pattern(tag, arg, mod) for arg in case.args]
                reduce(patternOr, pats)
            end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head == :(&&),
        rewrite = (tag, case, mod) ->
            let pats = [mk_pattern(tag, arg, mod) for arg in case.args]
                reduce(∘, pats)
            end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head == :(&),
        rewrite = (tag, case, mod) -> begin
                @assert length(case.args) == 1 "invalid ref pattern."
                var = case.args[1]
                body -> @format [tag, var, body] quote
                   tag == var ? body : failed
                end
        end
)

def_pattern(Pervasives,
        predicate = x -> x isa Symbol && x == :(_),
        rewrite = (_, _, _) -> identity
)

def_pattern(Pervasives,
        predicate = x -> x isa Symbol && x == :nothing,
        rewrite = mk_pat_by(===)
)

islower(s)::Bool = !isempty(s) && islowercase(s[1])
isupper(s)::Bool = !isempty(s) && isuppercase(s[1])

def_pattern(Pervasives,
        predicate = x -> x isa Symbol && islower ∘ string $ x,
        rewrite = (tag, case, mod) -> body ->
        @format [case, tag, body] quote
            let case = tag
                body
            end
        end
)

# Not decided yet about capitalized symbol's semantics, for generic enum is impossible in Julia.
def_pattern(Pervasives,
        predicate = x -> x isa Symbol && isupper ∘ string $ x,
        rewrite = (tag, case, mod) ->
        if used(:Enum, mod)
            mk_app_pattern(tag, case, [], mod)
        elseif used(:UppercaseCapturing, mod)
            body ->
            @format [case, tag, body] quote
                let case = tag
                    body
                end
            end
        else
            @syntax_err "You should use extension `Enum` or `UppercaseCapturing` to specify the behaviour of $(string(case))."
        end
)

def_pattern(Pervasives,
    predicate = x -> x isa Expr && x.head === :curly,
    rewrite = (tag, case, mod) -> mk_gapp_pattern(tag, [], case, [], mod)
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head == :if,
        rewrite   = (_, case, mod) ->
        # TODO: perform syntax validation here.
        let cond = case.args[1]
            body -> @format [cond, body] quote
                cond ? body : failed
            end
        end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head == :function,
        rewrite  = (tag, case, mod) ->
        let n = length(case.args)
            if n === 1
                fn = case.args[1]
            else
                fn = case
            end
            body -> @format [body, fn, tag] quote
                fn(tag) ? body : failed
            end
        end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head === :tuple,
        rewrite  = (tag, case, mod) ->
        let pat_elts = case.args,
            n = length(pat_elts),
            TARGET = mangle(mod),
            IDENTS = [mangle(mod) for _ in 1:n]
            map(1:n) do i
                IDENT = IDENTS[i]
                elt = pat_elts[i]
                function (body)
                    Expr(:block, Expr(:(=), IDENT, :($TARGET[$i])), body)
                end ∘ mk_pattern(IDENT, elt, mod)
            end |> match_elts ->
            let match_elts =  reduce(∘, match_elts, init=identity)
                (@typed_as NTuple{n, Any}) ∘ match_elts
            end
        end
)

def_pattern(Pervasives,
    predicate = x -> x isa Expr && x.head == :vect,
    rewrite   = (tag, case, mod) -> begin
        elts = case.args
        ordered_seq_match(tag, elts, mod)
    end
)

struct QuotePattern
        value
end

function mk_expr_template(expr :: Expr)
    if expr.head == :($)
        return expr.args[1]
    end
    rec = mk_expr_template
    Expr(:call, :Expr, rec(expr.head), filter(x -> x !== _discard, map(rec, expr.args))...)
end

struct Discard end
const _discard = Discard()

function mk_expr_template(expr :: Symbol)
        QuoteNode(expr)
end

function mk_expr_template(expr :: QuoteNode)
        QuotePattern(expr.value)
end

function mk_expr_template(expr :: LineNumberNode)
    _discard
end

function mk_expr_template(expr)
        expr
end

def_pattern(Pervasives,
       predicate = x -> x isa Expr && x.head == :quote,
       rewrite  = (tag, case, mod) -> begin
        expr = case.args[1]
        expr = mk_expr_template(expr)
        mk_pattern(tag, expr, mod)
       end
)

def_pattern(Pervasives,
       predicate = x -> x isa QuotePattern,
       rewrite  = (tag, case, mod) -> begin
        expr = case.value
        expr = mk_expr_template(expr)
        TARGET = mangle(mod)
        VALUE = mangle(mod)
        access_value(body) =
            @format [body, TARGET, VALUE] quote
                let VALUE = TARGET.value
                    body
                end
            end
        (@typed_as QuoteNode) ∘  access_value ∘ mk_pattern(VALUE ,expr, mod)
       end
)

def_app_pattern(Pervasives,
    predicate = (hd_obj, args) -> hd_obj === (:),
    rewrite = (tag, hd_obj, args, mod) -> begin
        pat = Expr(:call, hd_obj, args...)
        body -> @format [pat, tag, body] quote
            tag in pat ? body : failed
        end
    end
)

def_app_pattern(Pervasives,
    predicate = (hd_obj, args) -> hd_obj === Expr && !isempty(args),
    rewrite = (tag, hd_obj, args, mod) ->
    let TARGET     = mangle(mod)
        length(args) === 1 ?
        let arg = args[1]
            (arg isa Expr && arg.head === :...) ?
            let _          = (@assert length(arg.args) === 1),
                lst        = mangle(mod),
                arg        = arg.args[1],
                perf_match = mk_pattern(lst, arg, mod)
                exprargs_to_arr(body) =
                    @format [body, lst, TARGET] quote
                        lst = [TARGET.head, TARGET.args...]
                        body
                    end
                (@typed_as Expr) ∘ exprargs_to_arr ∘ perf_match
            end :
            let HEAD = mangle(mod),
                perf_match = mk_pattern(HEAD, arg, mod)
                bind_head(body) =
                    @format [body, HEAD, TARGET] quote
                        !isempty(TARGET.args) ? failed :
                        let HEAD = TARGET.head
                            body
                        end
                    end
                (@typed_as Expr) ∘ bind_head ∘ perf_match
            end
        end :
        let HEAD       = mangle(mod),
            ARGS       = mangle(mod)

            head_pat = args[1]
            args_pat = args[2:end]

            assign_attrs(body) =
                @format [body, HEAD, ARGS, TARGET] quote
                    let (HEAD, ARGS) = (TARGET.head, TARGET.args)
                        body
                    end
                end

            (@typed_as Expr) ∘ assign_attrs ∘ mk_pattern(HEAD, head_pat, mod) ∘ ordered_seq_match(ARGS, args_pat, mod)
        end
    end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head == :(::),
        rewrite = (tag, case, mod) ->
                let args   = (case.args..., ),
                    TARGET = mangle(mod)
                    function f(args :: NTuple{2, Any})
                        pat, t = args
                        (@typed_as t) ∘ mk_pattern(TARGET, pat, mod)
                    end

                    function f(args :: NTuple{1, Any})
                        t = args[1]
                        @typed_as t
                    end
                    f(args)
                end,
        qualifiers = Set([internal])
)

def_app_pattern(Pervasives,
        predicate = (hd_obj, args) -> hd_obj === Dict,
        rewrite = (tag, hd_obj, args, mod) ->
        let TARGET = mangle(mod)

            foldr(args, init=identity) do kv, last
                if !(isa(kv, Expr) && kv.head === :call && (@eval mod $(kv.args[1])) === Pair)
                    throw $
                    SyntaxError("Dictionary destruct must take patterns like Dict(<expr> => <pattern>, ...)")
                end
                let (k, v)    = kv.args[2:end],
                    IDENT     = mangle(mod),
                    match_elt = mk_pattern(IDENT, v, mod)
                    function (body)
                        @format [IDENT, TARGET, get, k, body] quote
                            IDENT = get(TARGET, k) do
                                failed
                            end
                            IDENT === failed ? failed : body
                        end
                    end ∘ match_elt ∘ last
                end
            end |> match_kvs ->
            (@typed_as Dict) ∘ match_kvs
        end
)

# arbitray ordered sequential patterns match
function ordered_seq_match(tag, elts, mod)

        TARGET = mangle(mod)
        NAME = mangle(mod)
        T = mangle(mod)
        function check_generic_array(body)
            @format [AbstractArray, NAME, TARGET, body, T, tag] quote

                @inline __L__ function NAME(TARGET :: AbstractArray{T, 1}) where {T}
                    body
                end

                @inline __L__ function NAME(_)
                    failed
                end

                NAME(tag)

            end
        end
        length(elts) == 0 ?
        (
            check_generic_array ∘
            function (body)
                @format [isempty, body, tag] quote
                    isempty(tag) ? body : failed
                end
            end
        )              :
        begin
            atleast_element_count = 0
            unpack_begin          = nothing
            unpack_end            = 0
            unpack                = []
            foreach(elts) do elt
                if elt isa Expr && elt.head === :...
                    if unpack_begin === nothing
                        unpack_begin = atleast_element_count + 1
                    else
                        throw $
                        SyntaxError("Sequential unpacking can only be performed once at most.")
                    end
                    push!(unpack, elt.args[1])
                else
                    atleast_element_count = atleast_element_count + 1
                    IDENT = mangle(mod)
                    index = unpack_begin === nothing ?
                        begin
                            :($TARGET[$atleast_element_count])
                        end                       :
                        begin
                            let exp = :($TARGET[end - $unpack_end])
                                unpack_end = unpack_end + 1
                                exp
                            end
                        end
                    if elt === :_
                        push!(unpack, identity)
                    else
                        perf_match = mk_pattern(IDENT, elt, mod)
                        push!(
                            unpack,
                            let IDENT = IDENT, index = index
                                function (body)
                                    @format [IDENT, body, index] quote
                                        IDENT = index
                                        body
                                    end
                                end ∘ perf_match
                            end
                        )
                    end
                end
            end
            if unpack_begin !== nothing
                IDENT = mangle(mod)
                check_len = body -> @format [body, TARGET, atleast_element_count, length] quote
                    length(TARGET) >= atleast_element_count ? body : failed
                end
                elt = unpack[unpack_begin]
                if elt === :_
                    unpack[unpack_begin] = identity
                else
                    unpack[unpack_begin] = function (body)
                        @format [body, IDENT, TARGET, unpack_begin, unpack_end, length] quote
                            IDENT = view(TARGET, unpack_begin: (length(TARGET) - unpack_end))
                            body
                        end
                    end ∘ mk_pattern(IDENT, elt, mod)
                end
            else
                check_len = body -> @format [body, TARGET, length, atleast_element_count] quote
                    length(TARGET) === atleast_element_count ? body : failed
                end
            end
            check_generic_array ∘ check_len ∘ foldr(∘, unpack)

        end
end

struct _ManyToken end
struct _DoToken end

export Many, Do
Many   = _ManyToken()
Do     = _DoToken()

function allow_assignment(expr :: Expr)
        head = expr.head == :kw ? :(=) : expr.head
        Expr(head, expr.args...)
end

function allow_assignment(expr)
        expr
end

def_app_pattern(Pervasives,
        predicate = (hd_obj, args) -> hd_obj === Do,
        rewrite = (_, _, args, _) ->
        let action = Expr(:block, map(allow_assignment, args)...)
                body -> @format [body, action] quote
                        (@inline __L__ function ()
                                action
                                body
                        end)()
                end
        end
)

def_app_pattern(Pervasives,
        predicate = (hd_obj, args) -> hd_obj === QuoteNode && length(args) == 1,
        rewrite = (tag, _, args, mod) ->
        begin
            VAR = mangle(mod)

            function (body)
                @format [tag, VAR, body, failed, QuoteNode] quote
                    if tag isa QuoteNode
                        VAR = tag.value
                        body
                    else
                        failed
                    end
                end
            end ∘
            mk_pattern(VAR, args[1], mod)

        end
)

def_app_pattern(Pervasives,
    predicate = (hd_obj, args) -> hd_obj === Many,
    rewrite   = (tag, hd_obj, args, mod) ->
    let inner = args[1],
        ITER_VAR = mangle(mod),
        TEST_VAR = mangle(mod),
        iter_check = mk_pattern(ITER_VAR, inner, mod)(nothing)

        @assert length(args) === 1 "syntax form should be `Many(pat)`."
        body -> @format [ITER_VAR, TEST_VAR, tag, iter_check, body] quote
            TEST_VAR = true
            for ITER_VAR in tag
                if iter_check !== nothing
                    TEST_VAR = false
                    break
                end
            end
            TEST_VAR ? body : failed
        end
    end
)

def_pattern(Pervasives,
        predicate = x -> x isa Expr && x.head === :where,
        rewrite = (tag, case, mod) ->
        begin
            if !used(:GADT, mod)
                @syntax_err "GADT extension hasn't been enabled. Try `@use GADT` and run your codes again."
            end
            # Not sure about if there's any other references of `where`,
            # but GADT is particularly important,
            # Tentatively, we use `where` for GADT support only.
            @match case begin
                :($hd($(tl...)) where {$(forall...)}) => mk_gapp_pattern(tag, forall, hd, tl, mod)
                _ =>
                    @syntax_err "Unknown usage of `where` in pattern region. Currently `where` is used for only GADT syntax."
            end
        end
)

function mk_gapp_pattern(tag, forall, hd, tl, use_mod)
        @match hd begin
            ::Symbol && if isempty(forall) end => mk_app_pattern(tag, hd, tl, use_mod)
            :($(ctor :: Symbol){$(spec_vars...)}) || ctor :: Symbol && Do(spec_vars = [])=>
                begin
                    if isdefined(use_mod, ctor)
                        let ctor = getfield(use_mod, ctor)
                            for (def_mod, desc) in Infras.GAPP_DESTRUCTORS
                                if qualifier_test(desc.qualifiers, use_mod, def_mod) && desc.predicate(spec_vars, ctor, tl)
                                    return desc.rewrite(tag, forall, spec_vars, ctor, tl, use_mod)
                                end
                            end
                        end
                    end
                    info = string(:($ctor{$(spec_vars...)}($(tl...))))
                    throw(PatternUnsolvedException("invalid usage or unknown application case $info."))
                end
        end
end

end

ext: jl


text:
const XIntV = Union{XInt, Nothing}

_vec(n::Integer) = Vector{Limb}(undef, n)
@inline vec(x::XInt) = x.v::Vector{Limb}

lenvec(x::XSigned) = len(x), vec(x)

vec!(x::XInt, n::Integer) =
    if is_short(x)
        _vec(n)
    else
        vec!(vec(x), n)
    end

function vec!(xv::Vector, n::Integer)
    len = length(xv)
    if n > len
        Base._growend!(xv, n-len)
    end
    xv
end

vec!(::Nothing, n::Integer) = _vec(n)

# copy x into r, allocating for r when r===nothing
# assumes !is_short(x)
_copy!(r::XIntV, x::XSigned) =
    if r === x
        x
    else
        xl, xv = lenvec(x)
        _XInt(flipsign(xl, _sign(x)), _copy!(vec!(r, xl), xv, 1:xl))
    end

function _copy!(r::Vector, xv::Vector, span::UnitRange)
    a, b = first(span), last(span)
    if r !== xv
        unsafe_copyto!(r, a, xv, a, b-a+1)
    end
    r
end

function _copy!(r::Vector, xv::SVec, span::UnitRange)
    for i=span
        @inbounds r[i] = xv[i]
    end
    r
end

@inline function assign!(r::XIntV, vals::T...) where T
    v = vec!(r, length(vals))
    for i = eachindex(vals)
        @inbounds v[i] = vals[i] % Limb
    end
    v
end

@inline XInt!(r::XIntV, z::Short) =
    fits(z) ? _XInt(z % SLimb) :
              _XInt(flipsign(one(SLimb), z), assign!(r, z))

# make a XInt into r, with value -z
@inline XInt_neg!(r::XIntV, z::Limb) =
    fits(z) ? _XInt(-z % SLimb) :
              _XInt(-one(SLimb), assign!(r, z))

@inline XInt!(r::XIntV, z::Limb, flip::SLimb) =
    flip < 0 ?
        XInt_neg!(r, z) :
        XInt!(r, z)

@inline XInt!(r::XIntV, z::SLimbW, extra::SLimb=0) =
    # @assert extra >= 0
    if slimbmin < z <= slimbmax
        _XInt(z % SLimb)
    else
        XInt_big!(r, z, extra)
    end

@noinline function XInt_big!(r::XIntV, z::SLimbW, extra::SLimb)
    # @assert extra >= 0
    zz = abs(z)
    z1 = zz % Limb
    z2 = (zz >>> BITS_PER_LIMB) % Limb
    if iszero(z2)
        rv = vec!(r, 1+extra)
        @inbounds rv[1] = z1
        _XInt(sign(z) % SLimb, rv)
    else
        rv = vec!(r, 2+extra)
        @inbounds rv[1] = z1
        @inbounds rv[2] = z2
        _XInt(flipsign(SLimb(2), z), rv)
    end
end

XInt!(r::XIntV, z::LimbW, extra::SLimb=0) =
    # @assert extra >= 0
    if z <= slimbmax
        _XInt(z % SLimb)
    else
        z1 = z % Limb
        z2 = (z >>> BITS_PER_LIMB) % Limb
        if iszero(z2)
            rv = vec!(r, 1+extra)
            @inbounds rv[1] = z1
            _XInt(SLimb(1), rv)
        else
            rv = vec!(r, 2+extra)
            @inbounds rv[1] = z1
            @inbounds rv[2] = z2
            _XInt(SLimb(2), rv)
        end
    end

@inline function normalize(v::Vector, l::Integer)
    # @assert l <= length(v)
    while l > 0
        iszero(@inbounds v[l]) || break
        l -= 1
    end
    l
end

function add1!(r::XIntV, x::XSigned, y::Short)
    # @assert !is_short(x)
    samesign = signbit(x) == signbit(y)
    if samesign
        _add1!(r, x, abs(y) % Limb)
    else
        _sub1!(r, x, abs(y) % Limb)
    end
end

function _add1!(r::XIntV, x::XSigned, y::Limb)
    xl = len(x)
    rv = vec!(r, xl+1)
    c = MPN.add_1(rv, x=>xl, y)
    @inbounds rv[xl+1] = c
    rl = xl + !iszero(c)
    _XInt(flipsign(rl, _sign(x)), rv)
end

function _sub1!(r::XIntV, x::XSigned, y::Limb)
    xl = len(x)
    # we know that abs(x) > y, as y comes from a SLimb, whose absolute value
    # is < |typemmin(SLimb)| == limb1min (upper bit of y is not set)
    if xl == 1
        xv = vec(x)
        # xv[1] >= limb1min by specification
        rv1 = @inbounds xv[1] - y
        XInt!(r, rv1, _sign(x))
    else
        # if xl > 1, we see that x-y can't be an immediate integer
        rv = vec!(r, xl)
        MPN.sub_1(rv, x=>xl, y)
        rl = xl-iszero(@inbounds rv[xl])
        _XInt(flipsign(rl, _sign(x)), rv)
    end
end

_widen(x::XSigned) = widen(short(x))
_promote(x::XInt) = x
_promote(x::BitInteger) = Static(x)

@inline add!(r::XIntV, x, y) = add!(r, _promote(x), _promote(y))

@inline function add!(r::XIntV, x::XSigned, y::XSigned)
    xshort = is_short(x)
    yshort = is_short(y)
    if xshort & yshort
        XInt!(r, _widen(x) + _widen(y))
    elseif xshort
        iszero(x) ? XInt(y) : add1!(r, y, short(x))
    elseif yshort
        iszero(y) ? XInt(x) : add1!(r, x, short(y))
    else
        addbig!(r, x, y)
    end
end

# fast path for sum
# this allocate in x even if x == 0
# (i.e. this doesn't return the other arg y in this case, so no aliasing)
@inline function addbig_positive!(x::XInt, y::XInt)
    if iszero(x.x)
        yl, yv = lenvec(y)
        xv = _vec(yl+2) # +1 for later invocation which requests it, and +1 for margin
        unsafe_copyto!(xv, 1, yv, 1, yl)
        return _XInt(yl, xv)
    end
    xl, xv = lenvec(x)
    yl, yv = lenvec(y)
    rv = xv # result in x, even if we have to swap the arguments
    if xl < yl
        xl, yl = yl, xl
        xv, yv = yv, xv
    end
    rv = vec!(rv, xl+1)
    c = MPN.add(rv, xv=>xl, yv=>yl) # c ∈ (0, 1)
    @inbounds rv[xl+1] = c
    _XInt(xl + c % SLimb, rv)
end

# NOTE: this is still unfortunately roughly 2x slower than MPZ.add! for BigInt
# we recover a lot of perfs if we inline instead, but that's a big function...
@inline function addbig!(r::XIntV, x::XSigned, y::XSigned)
    xl, xv = lenvec(x)
    yl, yv = lenvec(y)
    xz, yz = _sign(x), _sign(y)
    if xl < yl
        # we don't just exchange the variables within the main function to
        # avoid type-unstability
        addbig!(r, yl, yv, yz, xl, xv, xz)
    else
        addbig!(r, xl, xv, xz, yl, yv, yz)
    end
end

@noinline function addbig!(r::XIntV, xl, xv, xz, yl, yv, yz)
    if xz == yz
        rv = vec!(r, xl+1)
        c = MPN.add(rv, xv=>xl, yv=>yl) # c ∈ (0, 1)
        @inbounds rv[xl+1] = c
        _XInt(flipsign(xl + c % SLimb, xz), rv)
    elseif xl >= yl+2 # >= 3
        # we are sure that the result will have length xl or xl-1
        rv = vec!(r, xl)
        MPN.sub(rv, xv=>xl, yv=>yl)
        rl = xl - iszero(@inbounds rv[xl])
        _XInt(flipsign(rl, xz), rv)
    else
        local xvl, yvl
        flip = false
        if xl == yl
            # same length, need to compare the content; we could use:
            # MPN.cmp(x=>xl, y=>yl) % Int
            # but we can as well do it manually, to record an upper bound on rl
            while xl > 0
                xvl = @inbounds xv[xl]
                yvl = @inbounds yv[xl]
                if xvl != yvl
                    if xvl < yvl
                        flip = true
                    end
                    break
                end
                xl -= 1
            end
            if xl == 0
                return _XInt(0)
            elseif xl == 1
                if flip
                    return XInt!(r, yvl-xvl, yz)
                else
                    return XInt!(r, xvl-yvl, xz)
                end
            end
            yl = xl
        else # xl == yl+1
            xvl = @inbounds xv[xl]
            yvl = zero(Limb)
        end
        if flip
            _addbig!(r, yl, yv, yz, yvl, xl, xv, xz, xvl)
        else
            _addbig!(r, xl, xv, xz, xvl, yl, yv, yz, yvl)
        end
    end
end

@inline function _addbig!(r, xl, xv, xz, xvl, yl, yv, yz, yvl)
    # here, xvl > yvl && xl >= 2
    rvl = xvl - yvl # value of rv[xl] "so far"
    @inbounds if rvl == 1
        rv1, c = Base.sub_with_overflow(xv[1], yv[1])
        if c && fits(rv1) # above limbs could cancel out, check if result is immediate
            i = 2
            local ri::Limb
            while i < xl
                ri, c2 = Base.sub_with_overflow(xv[i], yv[i])
                ri, c = ri-c,
                        c2 | c & iszero(ri)
                iszero(ri) || break
                # at this point (ri==0), we can see that c has the same value as c2,
                # but letting c = c2 above is not enough, as we re-use c later for
                # the "real" computation, so it must be set according to the general
                # case
                i += 1
            end
            # i == xl || ri != 0
            if i == xl && c
                return _XInt(flipsign(rv1 % SLimb, xz)) # we know that rv1 fits
            else
                rv = vec!(r, xl)
                rv[1] = rv1
                fill!(@view(rv[2:i-1]), zero(Limb))
                if i == xl # !c, as the case c==true was handled above
                    rv[xl] = rvl
                    rl = xl
                else
                    rv[i] = ri
                    # we would like to use native MPN.sub_n, but we can't pass to it
                    # a pre-existing carry
                    c = MPN.sub_n(rv, xv, yv, i+1:xl-1, c)
                    rv[xl] = rvl-c
                    rl = normalize(rv, xl)
                end
                return _XInt(flipsign(rl, xz), rv)
            end
        end
    end
    @assert xl == yl || xl == yl+1
    rv = vec!(r, xl)
    c = MPN.sub_n(rv, xv, yv, yl)
    if xl > yl
        rv[xl] = @inbounds(xv[xl]-c)
    end
    rl = normalize(rv, xl)
    _XInt(flipsign(rl, xz), rv)
end

add!(x::XInt, y) = add!(x, x, _promote(y))

neg!(x::XInt) = _XInt(-x.x, x.v)

sub!(r::XIntV, x, y::XInt) = add!(r, _promote(x), neg!(y))
# we don't implement neg!(s::Static), as s might contain unsigned, or typemin(SLimb)
sub!(r::XIntV, x::XInt, y::BitInteger) = neg!(add!(r, neg!(x), Static(y)))

sub!(x::XInt, y) = sub!(x, x, _promote(y))

is_limb(x::XSigned) = is_short(x) # we don't bother for non-short with one limb
is_limb(x::ShortMax) = true
limb(x::XSigned) = short(x)
limb(x::SLimbMax) = x % SLimb
limb(x::ULimbMax) = x % Limb
_widen(x::ShortMax) = SLimbW(x)
_sign(x::Limb) = SLimb(1)
_sign(x::SLimb) = sign(x)

mul!(x::XInt, y::Union{ShortMax,XInt}) = mul!(x, x, y)

@inline function mul!(r::XIntV, x::XInt, y::Union{ShortMax,XInt})
    xshort = is_short(x)
    yshort = is_limb(y)
    if xshort & yshort
        XInt!(r, _widen(x) * _widen(y))
    elseif xshort
        mul1!(r, y, short(x))
    elseif yshort
        mul1!(r, x, limb(y))
    else
        mulbig!(r, x, y)
    end
end

function mul1!(r::XIntV, x::XInt, y::Short)
    # @assert !is_short(x) # y can be slimbmin
    iszero(y) && return XInt(0)
    isone(y) && return x
    y == -1 && return neg!(x)
    sz = flipsign(_sign(x), _sign(y))
    xl = len(x)
    rv = vec!(r, xl + 1)
    high = MPN.mul_1(rv, x, abs(y) % Limb) # rv === vec(x) is supported
    @inbounds rv[xl+1] = high
    _XInt(flipsign(xl + !iszero(high), sz), rv)
end

function mulbig!(r::XIntV, x::XInt, y::XInt)
    # @assert !is_short(x) && !is_short(y)
    xl, xv = lenvec(x)
    yl, yv = lenvec(y)
    sz = flipsign(_sign(x), _sign(y))
    if xl < yl
        xl, yl = yl, xl
        xv, yv = yv, xv
    end

    tv = rv = vec!(r, xl+yl)
    # check for aliasing
    if rv === xv
        # MPN can't do mpn_mul inplace; rv is already big enough to hold the result,
        # so copy xv to a temporary, which can be discarded afterwards
        tv = tmpvec(xl)
        unsafe_copyto!(tv, 1, xv, 1, xl)
        xv = tv
        if rv === yv
            # assume xl == yl; make yv egal again to xv, the square will be computed
            yv = xv
        end
    elseif rv === yv
        tv = tmpvec(yl)
        unsafe_copyto!(tv, 1, yv, 1, yl)
        yv = tv
    end
    high = MPN.mul(rv, xv => xl, yv => yl)
    tv !== rv && freetmpvec(tv)
    _XInt(flipsign(xl+yl-iszero(high), sz), rv)
end

# set r to ~x == -(x+1)
com!(r::XIntV, x::XInt=r) =
    if is_short(x)
        XInt!(r, ~x.x) # could be non-short for x.x == slimbmax
    else
        neg!(add1!(r, x, one(SLimb)))
    end

lshift!(x::XInt, c::Integer) = lshift!(x, x, c)

# from base/operators.jl
@inline function lshift!(r::XIntV, x::XInt, c::Integer)
    typemin(Cint) <= c <= typemax(Cint) && return lshift!(r, x, c % Cint)
    (x >= 0 || c >= 0) && return zero(XInt)
    XInt(-1)
end

@inline lshift!(r::XIntV, x::XInt, c::Unsigned) =
    c <= typemax(Cuint) ? lshift!(r, x, c % Cuint) : zero(XInt)

@inline lshift!(r::XIntV, x::XInt, c::Cint) =
    c >= 0 ? lshift!(r, x, unsigned(c)) : rshift!(r, x, unsigned(-c))

@inline function lshift!(r::XIntV, x::XInt, c::Cuint)
    z = x.x
    iszero(z) && return x # NOTE: needs updating if a reduce argument is added
    if is_short(x)
        iszero(c) && return x
        if z > 0
            leading_zeros(z) > c &&
                return _XInt(z << c)
        else
            if leading_ones(z) > c
                s = z << c
                s !== slimbmin &&
                    return _XInt(z << c)
                # NOTE: it seems to be almost 2x faster for the general case to not
                # return the (known) result here for s==slimbmin, but rather to let
                # this special case be handled by lshift_small!
            end
        end
        lshift_small!(r, z, c)
    else
        lshift_big!(r, x, c)
    end
end

@noinline function lshift_small!(r::XIntV, z::SLimb, c::Cuint)
    zu = abs(z) % Limb
    offset = _divLimb(c) % SLimb
    cnt = _modLimb(c) % SLimb
    rl = 1+offset
    if iszero(cnt)
        rv = vec!(r, rl)
        @inbounds rv[rl] = zu
    else
        rv = vec!(r, rl+1)
        w = widen(zu) << cnt
        @inbounds rv[rl] = w % Limb
        high = @inbounds rv[rl+1] = (w >> BITS_PER_LIMB) % Limb
        rl += !iszero(high)
    end
    fill!(@view(rv[1:offset]), Limb(0))
    _XInt(flipsign(rl, z), rv)
end

@noinline function lshift_big!(r::XIntV, x::XInt, c::Cuint)
    r === x && c == 0 && return x
    xl = len(x)
    xv = vec(x)
    offset = _divLimb(c) % SLimb
    cnt = _modLimb(c) % SLimb
    rl = xl+offset
    if iszero(cnt)
        rv = vec!(r, rl)
        # TODO: when r===x and rv requires allocation anyway in vec! above,
        # maybe it would be more efficient to prepend! zeros ?
        unsafe_copyto!(rv, offset+1, xv, 1, xl)
    else
        rv = vec!(r, rl+1)
        high = MPN.lshift(@view(rv[offset+1:rl]), xv=>xl, cnt)
        @inbounds rv[rl+1] = high
        rl += !iszero(high)
    end
    fill!(@view(rv[1:offset]), Limb(0))
    _XInt(flipsign(rl, x.x), rv)
end

rshift!(x::XInt, c::Integer) = rshift!(x, x, c)

# from base/operators.jl
@inline function rshift!(r::XIntV, x::XInt, c::Integer)
    typemin(Cint) <= c <= typemax(Cint) && return rshift!(r, x, c % Cint)
    (x >= 0 || c < 0) && return zero(XInt)
    XInt(-1)
end

@inline rshift!(r::XIntV, x::XInt, c::Cint) =
    c >= 0 ? rshift!(r, x, unsigned(c)) : lshift!(r, x, unsigned(-c))

@inline function rshift!(r::XIntV, x::XInt, c::Cuint)
    z = x.x
    is_short(x) && return _XInt(z >> c)
    (iszero(z) || iszero(c)) && return x
    rshift_big!(r, x, c)
end

# For non-negative integers, right-shift is as obvious as it gets.
# For negative numbers, we use as usual the formula ~y = -(y+1) to emulate 2-complement
# integers; assuming x < 0 and c > 0, letting u = -x and v = u - 1 >= 0, we have
# x == -(v + 1) == ~v, from which we infer
# x >> c == ~v >> c == ~(v >> c) == -(v >> c + 1)
# Now, we don't want to compute v >> c, only u >> c. There are two cases:
# 1. u is a multiple of 2^c, i.e. u = k*2^c and u >> c == k
#    then v == u-1 == (k-1)*2^c + (2^c-1), therefore v >> c == k-1
#    hence: x >> c == -(k-1 + 1) == -k == -(u >> c) == -(-x >> c)
# 2. otherwise, u = k*2^c + r where 1 <= r < 2^c
#    then v == u-1 == k*2^c + r-1 where 0 <= r-1 < 2^c, therefore v >> c == u >> c
#    hence: x >> c == -(-x >> c + 1)
# Also, this can be derived more direcly from x >> c == fld(x, 2^c) where x = -k*2^c - r
@noinline function rshift_big!(r::XIntV, x::XInt, c::Cuint)
    xl = len(x)
    xv = vec(x)
    offset = _divLimb(c) % SLimb
    cnt = _modLimb(c) % SLimb
    rl = xl-offset
    neg = signbit(x.x)
    rl <= 0 && return XInt(-neg)

    rest = false # records whether there is any set bit among the first c in -x (when x < 0)
    if neg
        rest = any(!=(zero(Limb)), @view(xv[1:offset]))
    end

    if  rl <= 2 # special case, to avoid allocating for nothing
        xvend = @inbounds xv[xl] % LimbW
        if rl == 2
            xvend = xvend << BITS_PER_LIMB | @inbounds xv[xl-1] % LimbW
        end
        u = xvend >> cnt
        if neg & !rest
            rest |= trailing_zeros(xvend) < cnt
        end
        if rest && u == typemax(LimbW)
            # happens e.g. with x = -(XInt(2)^192-1) and c == 64 (with sizeof(Limb) == 8)
            _XInt(SLimb(-3), assign!(r, 0, 0, 1))
        else
            u += rest
            if u <= slimbmax
                _XInt(flipsign(u % SLimb, x.x))
            elseif u <= typemax(Limb)
                _XInt(flipsign(one(SLimb), x.x),
                      assign!(r, u % Limb))
            else
                _XInt(flipsign(SLimb(2), x.x),
                      assign!(r, u % Limb , (u >> BITS_PER_LIMB) % Limb))
            end
        end
    else
        rv = vec!(r, rl + neg)
        if iszero(cnt)
            unsafe_copyto!(rv, 1, xv, offset+1, rl)
        else
            discarded = MPN.rshift(rv=>rl, @view(xv[offset+1:xl]), cnt)
            if neg
                rest |= !iszero(discarded)
            end
            rl -= iszero(@inbounds rv[rl])
        end
        if rest
            d = MPN.add_1(rv, rv=>rl, Limb(1))
            @inbounds rv[rl+1] = d
            rl += !iszero(d)
        end
        _XInt(flipsign(rl, x.x), rv)
    end
end

@inline and!(r::XIntV, x, y) = and!(r, _promote(x), _promote(y))

@inline and!(x::XInt, y::Union{ShortMax,XSigned}) = and!(x, x, y)

@inline function and!(r::XIntV, x::Union{ShortMax,XSigned}, y::Union{ShortMax,XSigned})
    xshort = is_limb(x)
    yshort = is_limb(y)
    if xshort & yshort
        XInt!(r, limb(x) & limb(y)) # if y is unsigned, result of & is also unsigned and
                                    # XInt! does the right thing
    elseif xshort
        iszero(x) ? _XInt(0) : and_small!(r, y, limb(x))
    elseif yshort
        iszero(y) ? _XInt(0) : and_small!(r, x, limb(y))
    else
        and_big!(r, x, y)
    end
end

function and_small!(r::XIntV, x::XInt, y::Limb)
    xv = vec(x)
    x1 = @inbounds xv[1]
    if ispos(x)
        XInt!(r, x1 & y)
    else
        # x & y = ~(-x - 1) & y
        XInt!(r, ~(x1-1) & y)
    end
end

@noinline function and_small!(r::XIntV, x::XSigned, y::SLimb)
    # @assert !iszero(x) && !iszero(y) && !is_short(x) && is_short(y)
    # y is short, but can also be typemin(SLimb)
    xl, xv = lenvec(x)
    yy = y % Limb
    x1 = @inbounds xv[1]
    if ispos(x)
        if y > 0
            # x & y <= min(x, y) == y, so the result is short
            _XInt(x1 & yy)
        else # x > 0, y < 0
            u = x1 & yy
            res = _copy!(r, x)
            rv = vec(res)
            @inbounds rv[1] = u
            # @assert_reduced x # if x is not reduced, the result could be immediate
            res
        end
    elseif y > 0 # x < 0
        # x & y = ~(-x - 1) & y; result >= 0, and <= y, so is immediate
        _XInt((~(x1-1) & yy) % SLimb)
    else # x < 0, y < 0
        # let u = -x: -u & y = ~(u-1) & y = ~(u-1 | ~y) = -(u-1 | ~y) - 1
        # the result is < 0, and it's easy to see that it's also <= min(x, y), so
        # it can't be an immediate integer (provided x isn't)

        if iszero(x1)
            # if x1 == 0, we see that a 2-complement representation of x would also
            # have a null first limb, so x & y == x
            XInt(x)
        else
            # x1 > 0, so x-1 differs from x only in the 1st limb
            r1 = (x1-1) | ~y
            if r1 == typemax(Limb)
                # +1 overflows, so we need to allocate one more limb, just in case
                rv = vec!(r, xl+1)
                @inbounds rv[1] = zero(Limb)
                carry = true
                @inbounds for i=2:xl
                    u = xv[i]
                    if carry && u == typemax(Limb)
                        rv[i] = zero(Limb)
                    else
                        rv[i] = u+carry
                        carry = false
                    end
                end
                @inbounds rv[xl+1] = carry
                rl = xl + carry
                _XInt(-rl, rv)
            else
                rv = vec!(r, xl)
                _copy!(rv, xv, 2:xl)
                @inbounds rv[1] = r1+1
                _XInt(flipsign(xl, _sign(x)), rv)
            end
        end
    end
end

@noinline function and_big!(r::XIntV, x::XSigned, y::XSigned)
    # @assert !iszero(x) && !iszero(y) && !is_short(x) && !is_short(y)
    if ispos(x)
        if ispos(y)
            xl, xv = lenvec(x)
            yl, yv = lenvec(y)
            if xl < yl
                xv, yv = yv, xv
                xl, yl = yl, xl
            end
            rl = yl
            local rvl # stores the last limb in r
            @inbounds while rl > 0
                rvl = xv[rl] & yv[rl]
                iszero(rvl) || break
                rl -= 1
            end
            if rl <= 1
                XInt!(r, rvl)
            else
                rv = vec!(r, rl)
                # MPN.and_n(rv, xv, yv, rl)
                for idx = 1:rl
                    # we already computed the value of the last iteration, but
                    # it seems no more expensive to just recompute it in this loop
                    @inbounds rv[idx] = xv[idx] & yv[idx]
                end
                _XInt(rl, rv)
            end
        else # x > 0, y < 0
            x, y = y, x
            @goto mixed
        end
    elseif ispos(y) # x < 0
        @label mixed
        # result is positive
        # x & y = ~(-x - 1) & y
        xl, xv = lenvec(x)
        yl, yv = lenvec(y)
        if yl <= xl
            # we don't know the resulting rl in advance; on the high limbs of -x,
            # we don't know yet the effect of `-1`. So we first scan from the low
            # end: as long as x's limbs are 0, they become 0xF...F after `-1`, and
            # 0 after `~`, so r's corresponding limbs are 0.
            i = 1
            local xi
            while i <= yl
                xi = @inbounds xv[i]
                iszero(xi) || break
                i += 1
            end
            if yl < i
                _XInt(0)
            else
                rl = yl
                local rvl
                # we scan from the high end to compute rl: above i, the `-1` had no effect,
                # and we search for the hightest non-zero limb in r
                while rl > i
                    rvl = @inbounds ~xv[rl] & yv[rl]
                    iszero(rvl) || break
                    rl -= 1
                end
                if rl == i
                    # all the limbs of r above and below i are 0
                    rvl = @inbounds ~(xi-1) & yv[i]
                    if iszero(rvl)
                        rl = 0
                    end
                end
                if rl <= 1
                    XInt!(r, rvl)
                else
                    rv = vec!(r, rl)
                    for j = 1:i-1
                        @inbounds rv[j] = zero(Limb)
                    end
                    if i < rl
                        @inbounds rv[i] = ~(xv[i]-1) & yv[i] # rvl if i == rl
                        for j = i+1:rl-1
                            @inbounds rv[j] = ~xv[j] & yv[j]
                        end
                    end
                    @inbounds rv[rl] = rvl
                    _XInt(rl, rv)
                end
            end
        else # yl > xl
            # no need to scan, we know that rl == yl
            rv = vec!(r, yl)
            xok = false
            @inbounds for i = 1:xl
                xi = xv[i]
                yi = yv[i]
                if xok
                    rv[i] = ~xi & yi
                else
                    xok = !iszero(xi)
                    rv[i] = ~(xi-1) & yi
                end
            end
            for i = xl+1:yl
                @inbounds rv[i] = yv[i]
            end
            _XInt(yl, rv)
        end
    else # x < 0, y < 0
        # result is negative
        # -u & -v = ~(u-1) & ~(v-1) = ~(u-1 | v-1) = -(u-1 | v-1) - 1
        xl, xv = lenvec(x)
        yl, yv = lenvec(y)
        if xl < yl
            xv, yv = yv, xv
            xl, yl = yl, xl
        end
        rl = xl
        rv = vec!(r, rl+1)
        xok, yok, oneok = false, false, false
        for i=1:xl
            xi = @inbounds xv[i]
            yi = i <= yl ? @inbounds(yv[i]) : zero(Limb)
            if !xok
                xok = !iszero(xi) # carry won't propagate anymore
                xi -= 1
            end
            if !yok
                yok = !iszero(yi)
                yi -= 1
            end
            u = xi | yi
            if !oneok
                oneok = u != typemax(Limb)
                u += 1
            end
            @inbounds rv[i] = u
        end
        if !oneok
            rl += 1
            @inbounds rv[rl] = one(Limb)
        end
        _XInt(-rl, rv)
    end
end

## ior (|)

@inline ior!(r::XIntV, x, y) = ior!(r, _promote(x), _promote(y))
@inline ior!(x::XInt, y::Union{ShortMax,XSigned}) = ior!(x, x, y)

@inline function ior!(r::XIntV, x::Union{ShortMax,XSigned}, y::Union{ShortMax,XSigned})
    xshort = is_limb(x)
    yshort = is_limb(y)
    if xshort & yshort
        x1 = limb(x)
        y1 = limb(y)
        r1 = x1 | y1
        if x1 >= 0 && y1 >= 0 # result is >= 0
            XInt!(r, unsigned(r1))
        else # result is < 0
            XInt!(r, signed(r1))
        end
    elseif xshort
        iszero(x) ? XInt(y) : ior_small!(r, y, limb(x))
    elseif yshort
        iszero(y) ? XInt(x) : ior_small!(r, x, limb(y))
    else
        ior_big!(r, x, y)
    end
end

@noinline function ior_small!(r::XIntV, x::XSigned, y::Short)
    # @assert !iszero(x) && !iszero(y) && !is_short(x) && is_short(y)
    # y is small
    xl, xv = lenvec(x)
    yy = y % Limb
    x1 = @inbounds xv[1]
    if ispos(x)
        if y > 0
            res = _copy!(r, x)
            rv = vec(res)
            @inbounds rv[1] = x1 | yy
            res
        else # x > 0, y < 0
            # when fits(y), y is short, so more than its sign bit is set,
            # and x | y >= y (result is short); but when y doesn't come
            # from a XInt, y==slimbmin is possible, so result is not
            # guarranteed to be short
            # _XInt(x1 | yy) # when fits(y)
            XInt!(r, (x1 % SLimb) | y)
        end
    elseif y > 0 # x < 0
        # x | y = ~(-x - 1) | y = ~((-x - 1) & ~y) = -((-x - 1) & ~y) - 1
        carry = iszero(x1)
        r1 = (x1-1) & ~yy + 1
        # a small analysis shows that, given y > 0, the `+1` can't carry over upper limbs

        if xl == 1
            # then !carry as -x > 0
            XInt_neg!(r, r1)
        else
            x2 = @inbounds xv[2]
            r2 = x2 - carry
            carry &= iszero(x2)
            if xl == 2 && iszero(r2)
                XInt_neg!(r, r1)
            else
                rv = vec!(r, xl)
                @inbounds rv[1] = r1
                @inbounds rv[2] = r2
                @inbounds for i=3:xl
                    xi = xv[i]
                    rv[i] = xi - carry
                    carry &= iszero(xi)
                end
                rl = xl - iszero(@inbounds rv[xl])
                _XInt(-rl, rv)
            end
        end
    else # x < 0, y < 0
        # x | y = ~(-x-1) | y # similar to x > 0, y < 0
        # _XInt(~(x1-1) | yy)
        XInt!(r, (~(x1-1) | yy) % SLimb)
    end
end

@noinline function ior_big!(r::XIntV, x::XSigned, y::XSigned)
    # @assert !iszero(x) && !iszero(y) && !is_short(x) && !is_short(y)
    # we keep here (and in and_big!) some type instabilities if
    # typeof(x) != typeof(y), but the overhead is not that bad, and might
    # not be concerning enough to warrant type-stabilizing this function
    # (like we did for addbig!)
    if ispos(x)
        if ispos(y) # r > 0
            xl, xv = lenvec(x)
            yl, yv = lenvec(y)
            if xl < yl
                xv, yv = yv, xv
                xl, yl = yl, xl
            end
            rv = vec!(r, xl)
            for i=1:yl
                @inbounds rv[i] = xv[i] | yv[i]
            end
            for i=yl+1:xl
                @inbounds rv[i] = xv[i]
            end
            _XInt(xl, rv)
        else # x > 0, y < 0 : r < 0
            x, y = y, x
            @goto mixed
        end
    elseif ispos(y) # x < 0 : r < 0
        @label mixed
        # x | y = ~(-x - 1) | y = ~((-x - 1) & ~y) = -((-x - 1) & ~y) - 1
        xl, xv = lenvec(x)
        yl, yv = lenvec(y)
        if yl + 2 <= xl
            # then -x-1 still has at least one more limb than y, which is not cancelled
            # when &-ed with ~y; then rl could be xl, or xl-1 when xv[l] == 1 and
            # xv[l-1] == 0
            rv = vec!(r, xl)
            cm = cp = true # carries for -1 and +1
            @inbounds for i=1:yl
                xi = xv[i]
                yi = yv[i]
                ri = (xi - cm) & ~yi + cp
                cm &= iszero(xi)
                cp &= iszero(ri) # note: !cp ==> !cm, i.e. cm can't last longer than cp
                rv[i] = ri
            end
            # at this point, we know for sure that !cp
            # bewond y's limbs, ~y is all ones, so we just write x limbs, carrying out cm
            for i=yl+1:xl
                xi = xv[i]
                rv[i] = xi - cm
                cm &= iszero(xi)
            end
            rl = xl - iszero(@inbounds rv[xl])
            _XInt(-rl, rv)
        else # yl >= xl - 1
            # rl could be anything, e.g.  -(2^128) | (2^128-1) == -1
            # so the idea is to scan from the end to find the highest non-zero limb
            # but for this, we need to know where the carry from `-1` stops: this is the
            # first non-zero position in xv, call it k; we therefore need to first scan
            # from the beginning to find k, and then from the end with the
            # formula xi & ~yi; if it happens that this is all zero beyond k, we would
            # then need to scan again from the beginning to find rl, before allocating;
            # this is a rare case, so we might not be concerned with having to scan
            # twice from the beginning; but for simplicity, in the first scan from the
            # beginning, we also record the highest non-zero ri upto k (this is a bit
            # more expensive than just finding the first non-zero xi, but in general
            # this loop won't have many iterations)
            cm = cp = true # carries for -1 and +1
            i = 0 # the first non-zero position in xv
            rl = 0
            @inbounds while cm
                i += 1
                xi = xv[i]
                yi = i > yl ? zero(Limb) : yv[i]
                ri = (xi - cm) & ~yi + cp
                cm &= iszero(xi)
                cp &= iszero(ri)
                if !iszero(ri)
                    rl = i
                end
            end
            k = i
            # beyond i, the formula is xi & ~yi: we scan from the end to find the first non-0
            i = xl
            @inbounds while i > k
                yi = i > yl ? zero(Limb) : yv[i]
                iszero(xv[i] & ~yi) || break
                i -= 1
            end
            # so far, rl points at the highest non-zero limb upto k in the result
            if i != k # not all zero above k
                rl = i
            end
            # we are now ready to store the result
            if rl == 1
                r1 = @inbounds (xv[1]-1) & ~yv[1] + 1
                XInt_neg!(r, r1)
            else
                rv = vec!(r, rl)
                cm = cp = true # carries for -1 and +1
                @inbounds for i=1:rl
                    xi = xv[i]
                    yi = i > yl ? zero(Limb) : yv[i]
                    ri = (xi - cm) & ~yi + cp
                    cm &= iszero(xi)
                    cp &= iszero(ri)
                    rv[i] = ri
                end
                _XInt(-rl, rv)
            end
        end
    else # x < 0, y < 0 : r < 0
        # -u | -v = ~(u-1) | ~(v-1) = ~(u-1 & v-1) = -(u-1 & v-1) - 1
        # we do a scan, similar to what we do for x < 0, y > 0 and yl >= xl - 1
        xl, xv = lenvec(x)
        yl, yv = lenvec(y)
        if xl < yl
            xv, yv = yv, xv
            xl, yl = yl, xl
        end

        cx, cy, cp = true, true, true
        i = 0
        rl = 0
        while (cx | cy) && i < yl
            i += 1
            xi = @inbounds xv[i]
            yi = @inbounds yv[i]
            u = (xi-cx) & (yi-cy) + cp
            if !iszero(u)
                rl = i
            end
            cx &= iszero(xi)
            cy &= iszero(yi)
            cp &= iszero(u) # cp = cx & cy would also work
        end
        k = i
        i = yl
        @inbounds while i > k
            iszero(xv[i] & yv[i]) || break
            i -= 1
        end
        if i > k
            rl = i
        end
        if rl == 1
            r1 = @inbounds (xv[1]-1) & (yv[1]-1) + 1
            XInt_neg!(r, r1)
        else
            rv = vec!(r, rl)
            cx, cy, cp = true, true, true
            @inbounds for i=1:rl
                xi = xv[i]
                yi = yv[i]
                rv[i] = u = (xi-cx) & (yi-cy) + cp
                cx &= iszero(xi)
                cy &= iszero(yi)
                cp &= iszero(u)
            end
            _XInt(-rl, rv)
        end
    end
end

@inline xor!(r::XIntV, x, y) = xor!(r, _promote(x), _promote(y))
@inline xor!(x::XInt, y::Union{ShortMax,XSigned}) = xor!(x, x, y)

@inline function xor!(r::XIntV, x::Union{ShortMax,XSigned}, y::Union{ShortMax,XSigned})
    xshort = is_limb(x)
    yshort = is_limb(y)
    if xshort & yshort
        x1 = limb(x)
        y1 = limb(y)
        XInt!(r, _widen(x1) ⊻ _widen(y1))
    elseif xshort
        iszero(x) ? XInt(y) : xor_small!(r, y, limb(x))
    elseif yshort
        iszero(y) ? XInt(x) : xor_small!(r, x, limb(y))
    else
        xor_big!(r, x, y)
    end
end

@noinline function xor_small!(r::XIntV, x::XSigned, y::Limb)
    xl, xv = lenvec(x)
    x1 = @inbounds xv[1]
    if ispos(x)
        if xl == 1
            XInt!(r, x1 ⊻ y)
        else
            res = _copy!(r, x)
            rv = vec(res)
            @inbounds rv[1] = x1 ⊻ y
            res
        end
    else
        # x ⊻ y = ~(~x ⊻ y) = ~((-x-1) ⊻ y) = -((-x-1) ⊻ y) - 1
        cm = iszero(x1)
        r1 = (x1-1) ⊻ y + 1
        cp = iszero(r1) # if cm, then !cp (as y != 0)
        if !cp && xl == 1
            # xl == 1 ==> !cm
            XInt_neg!(r, r1)
        elseif xl == 2
            x2 = @inbounds xv[2]
            r2 = x2 - cm + cp
            if iszero(r2)
                if cp
                    _XInt(SLimb(-3), assign!(r, r1, r2, one(Limb)))
                else
                    XInt_neg!(r, r1)
                end
            else
                _XInt(SLimb(-2), assign!(r, r1, r2))
            end
        else
            rv = vec!(r, xl+1)
            @inbounds rv[1] = r1
            @inbounds for i=2:xl
                xi = xv[i]
                rv[i] = xi - cm + cp
                cm &= iszero(xi)
                cp &= xi == typemax(Limb)
            end
            @inbounds rv[xl+1] = cp
            _XInt(-(xl + !iszero(cp)), rv)
        end
    end
end

@noinline function xor_small!(r::XIntV, x::XSigned, y::SLimb)
    # @assert !iszero(x) && !iszero(y) && !is_short(x) && is_short(y)
    xl, xv = lenvec(x)
    yy = y % Limb
    x1 = @inbounds xv[1]
    if ispos(x)
        if y > 0
            res = _copy!(r, x)
            rv = vec(res)
            @inbounds rv[1] = x1 ⊻ yy
            res
        else # x > 0, y < 0
            # x ⊻ y = ~(x ⊻ ~y) = -(x ⊻ ~y) - 1
            # if x is one limb, then the high bit of its limb is set,
            # but the high bit of ~y is not, so the result can't be immediate
            rv = vec!(r, xl+1) # ?TODO: avoid +1 in the majority of cases where its not needed?
            @inbounds rv[1] = r1 = x1 ⊻ ~yy + 1
            c = iszero(r1)
            for i=2:xl
                xi = @inbounds xv[i]
                @inbounds rv[i] = xi + c
                c &= xi == typemax(Limb)
            end
            @inbounds rv[xl+1] = c
            _XInt(-(xl + !iszero(c)), rv)
        end
    elseif y > 0 # x < 0
        # x ⊻ y = ~(~x ⊻ y) = ~((-x-1) ⊻ y) = -((-x-1) ⊻ y) - 1
        cm = iszero(x1)
        r1 = (x1-1) ⊻ yy + 1
        cp = iszero(r1) # if cm, then !cp
        if !cp && xl == 1
            # xl == 1 ==> !cm
            XInt_neg!(r, r1)
        else
            rv = vec!(r, xl+1)
            @inbounds rv[1] = r1
            @inbounds for i=2:xl
                xi = xv[i]
                rv[i] = xi - cm + cp
                cm &= iszero(xi)
                cp &= xi == typemax(Limb)
            end
            @inbounds rv[xl+1] = cp
            rl =
                if cp
                    xl + 1
                elseif iszero(@inbounds rv[xl])
                    # can happen because of cm, but then !iszero(rv[xl-2])
                    xl - 1
                else
                    xl
                end
            _XInt(-rl, rv)
        end
    else # x < 0, y < 0 (result >= 0)
        # x ⊻ y = ~x ⊻ ~y = (-x - 1) ⊻ ~y
        r1 = (x1-one(Limb)) ⊻ ~yy
        if xl == 1
            # res is immediate if x1 == limb1min
            XInt!(r, r1)
        else
            rv = vec!(r, xl)
            @inbounds rv[1] = r1
            c = iszero(x1)
            local ri
            for i=2:xl
                xi = @inbounds xv[i]
                @inbounds ri = rv[i] = xi-c
                c &= iszero(xi)
            end
            _XInt(xl - iszero(ri), rv)
        end
    end
end

@noinline function xor_big!(r::XIntV, x::XSigned, y::XSigned)
    # @assert !iszero(x) && !iszero(y) && !is_short(x) && !is_short(y)
    xl, xv = lenvec(x)
    yl, yv = lenvec(y)
    if ispos(x)
        if ispos(y)
            if xl < yl
                xl, yl = yl, xl
                xv, yv = yv, xv
            end
            if xl == yl
                local i = xl
                local rvl
                while i > 0
                    xi = @inbounds xv[i]
                    yi = @inbounds yv[i]
                    rvl = xi ⊻ yi
                    iszero(rvl) || break
                    i -= 1
                end
                if i == 0
                    XInt(0)
                elseif i == 1
                    XInt!(r, rvl)
                else
                    rv = vec!(r, i)
                    rl = i
                    @inbounds rv[rl] = rvl
                    i -= 1
                    while i > 0
                        xi = @inbounds xv[i]
                        yi = @inbounds yv[i]
                        @inbounds rv[i] = xi ⊻ yi
                        i -= 1
                    end
                    _XInt(rl, rv)
                end
            else
                rv = vec!(r, xl)
                for i=1:yl
                    @inbounds rv[i] = xv[i] ⊻ yv[i]
                end
                for i=yl+1:xl
                    @inbounds rv[i] = xv[i]
                end
                _XInt(xl, rv)
            end
        else # x > 0, y < 0
            xl, yl = yl, xl
            xv, yv = yv, xv
            @goto mixed
        end
    elseif ispos(y) # x < 0
        @label mixed
        # x ⊻ y = ~(~x ⊻ y) = ~((-x-1) ⊻ y) = -((-x-1) ⊻ y) - 1
        if yl > xl
            rv = vec!(r, yl + 1)
            cp = cm = true
            @inbounds for i=1:xl
                xi = xv[i]
                yi = yv[i]
                ri = rv[i] = (xi-cm) ⊻ yi + cp
                cm &= iszero(xi)
                cp &= iszero(ri)
            end
            @inbounds for i=xl+1:yl
                ri = rv[i] = yv[i] + cp
                cp &= iszero(ri)
            end
            @inbounds rv[yl+1] = cp
            _XInt(-(yl + !iszero(cp)), rv)
        elseif yl + 2 <= xl
            rv = vec!(r, xl + 1)
            cp = cm = true
            @inbounds for i=1:yl
                xi = xv[i]
                yi = yv[i]
                ri = rv[i] = (xi-cm) ⊻ yi + cp
                cm &= iszero(xi)
                cp &= iszero(ri)
            end
            @inbounds for i=yl+1:xl
                xi = xv[i]
                rv[i] = xi-cm+cp
                cm &= iszero(xi)
                cp &= iszero(ri)
            end
            @inbounds rv[xl+1] = cp
            _XInt(-(xl + !iszero(cp)), rv)
        else # xl == yl || xl == yl+1
            x1 = @inbounds xv[1]
            y1 = @inbounds yv[1]
            r1 = (x1-1) ⊻ y1 + 1
            cm = iszero(x1)
            cp = iszero(r1)
            # when fits(r1), all upper limbs could be 0 and result could be immediate,
            # so we scan from the low end to check that out
            i = 2
            local ri
            @inbounds while i <= yl
                xi = xv[i]
                yi = yv[i]
                ri = (xi-cm) ⊻ yi + cp
                cm &= iszero(xi)
                cp &= iszero(ri)
                iszero(ri) || break
                i += 1
            end
            if i > yl && xl > yl
                ri = @inbounds xv[xl] - cm + cp
                cp &= iszero(ri)
                if iszero(ri)
                    i += 1
                else
                    cp = false
                end
            end
            if !cp && i > xl
                XInt_neg!(r, r1)
            else
                rv = vec!(r, xl + (i > xl))
                @inbounds rv[1] = r1
                fill!(@view(rv[2:i-1]), zero(Limb))
                if i <= xl # then !cp
                    @inbounds rv[i] = ri
                    rl = i
                    i += 1
                    @inbounds while i <= yl
                        xi = xv[i]
                        yi = yv[i]
                        ri = rv[i] = (xi-cm) ⊻ yi
                        cm &= iszero(xi)
                        if !iszero(ri)
                            rl = i
                        end
                        i += 1
                    end
                    if xl > yl && i <= xl
                        @inbounds ri = rv[i] = xv[i]-cm
                        if !iszero(ri)
                            rl = i
                        end
                    end
                    _XInt(-rl, rv)
                else
                    @inbounds rv[xl+1] = cp
                    _XInt(-SLimb(xl+1), rv)
                end
            end
        end
    else # x < 0, y < 0 (result >= 0)
        # x ⊻ y = ~x ⊻ ~y = (-x - 1) ⊻ (-y - 1)
        if xl < yl
            xl, yl = yl, xl
            xv, yv = yv, xv
        end
        @inbounds x1 = xv[1]
        @inbounds y1 = yv[1]
        r1 = (x1-1) ⊻ (y1-1)
        cx = iszero(x1)
        cy = iszero(y1)
        i = 2
        local ri
        while i <= xl
            @inbounds xi = xv[i]
            @inbounds yi = i <= yl ? yv[i] : zero(Limb)
            ri = (xi-cx) ⊻ (yi-cy)
            cx &= iszero(xi)
            cy &= iszero(yi)
            iszero(ri) || break
            i += 1
        end
        if i > xl
            XInt!(r, r1)
        else
            rv = vec!(r, xl)
            @inbounds rv[1] = r1
            fill!(@view(rv[2:i-1]), zero(Limb))
            @inbounds rv[i] = ri
            rl = i
            @inbounds for i=i+1:xl
                xi = xv[i]
                yi = i <= yl ? yv[i] : zero(Limb)
                ri = rv[i] = (xi-cx) ⊻ (yi-cy)
                cx &= iszero(xi)
                cy &= iszero(yi)
                if !iszero(ri)
                    rl = i
                end
            end
            _XInt(rl, rv)
        end
    end
end


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

deps = [
    "https://use.fontawesome.com/releases/v5.0.7/js/all.js",
    "https://raw.githubusercontent.com/piever/InteractResources/master/highlight/prism.css",
    "https://raw.githubusercontent.com/piever/InteractResources/master/highlight/prism.js",
    "https://cdn.bootcss.com/noUiSlider/11.1.0/nouislider.min.js",
    "https://cdn.bootcss.com/noUiSlider/11.1.0/nouislider.min.css",
]

for dep in deps
    download(dep, joinpath(_pkg_assets, splitdir(dep)[2]))
end

using NodeJS
npm_path = joinpath(_pkg_assets, "npm")
mkpath(npm_path)
cd(npm_path) do
    npm = NodeJS.npm_cmd()
    run(`$npm install -y`)
    run(`$npm install -y katex`)
end


text: dim(data::Array) = size(data)


text:
"""
    rule2_applicable(G::SimpleDiGraph, Y, Z, X, W)::Bool

Return whether rule 2 is applicable given graph `G` with Y, Z, X and W.
In other words, return true if ``(Y \\perp Z | X, W)_{G_{\\overline{x} \\underline{z}}}``.
"""
function rule2_applicable(G::AbstractGraph, Y, Z, X, W)::Bool
    G = without_incoming(G, Set(X))
    G = without_outgoing(G, Set(Z))
    return d_separated(G, Y, Z, Set(X))
end

sym2uppercase(s::Symbol) = Symbol(uppercase(string(s)))
name2node(G, s::SymbolicUtils.Sym) = name2node(G, sym2uppercase(s.name))
name2node(G, S::AbstractArray) = Set([name2node(G, s) for s in S])
name2node(G, s::Set) = s

function rule2(G::MetaDiGraph, y, z, x, w)
    Y = name2node(G, y)
    Z = name2node(G, z)
    X = name2node(G, x)
    W = name2node(G, w)
    out = rule2_applicable(G, Y, Z, X, W)
    return out
end

function rule2(G::MetaDiGraph)
    # Make sure to put a simplified version at the lhs
    # This even holds when not using the infix operator.
    # So, +(a, b), in nested expressions, is not always the same as +(b, a).
    # (https://github.com/JuliaSymbolics/SymbolicUtils.jl/issues/331).

    forms = [
        # Rule 2 (Pearl, 2009; eq. 3.32) where x is empty.
        @acrule(P(~y, +(~~w, Do(~z))) => rule2(G, ~y, ~z, ~~w, Set()) ? P(~y, +(~z, ~~w...)) : nothing),
        # Rule 2 (Pearl, 2009; eq. 3.32) where x and w are empty.
        @acrule(P(~y, Do(~z)) => rule2(G, ~y, ~z, Set(), Set()) ? P(~y, ~z) : nothing)
    ]

    return SU.Chain(forms)
end


text:
using StatsBase
import DSP
import DSP: Periodograms
import LombScargle
export estimate_period

"""
    estimate_period(v::Vector, method, t=0:length(v)-1; kwargs...)
Estimate the period of the signal `v`, with accompanying time vector `t`,
using the given `method`.

If `t` is an AbstractArray, then it is iterated through to ensure that it's
evenly sampled (if necessary for the algorithm).  To avoid this, you can pass
any `AbstractRange`, like a `UnitRange` or a `LinRange`, which are defined to be
evenly sampled.

## Methods requiring evenly sampled data

These methods are faster, but some are error-prone.

* `:periodogram` or `:pg`: Use the fast Fourier transform to compute a
   periodogram (power-spectrum) of the given data.  Data must be evenly sampled.

* `:multitaper` or `mt`: The multitaper method reduces estimation bias by using
  multiple independent estimates from the same sample. Data tapers are then
  windowed and the power spectra are obtained.  Available keywords follow:
  `nw` is the time-bandwidth product, and `ntapers` is the number of tapers.
  If `window` is not specified, the signal is tapered with `ntapers` discrete
  prolate spheroidal sequences with time-bandwidth product `nw`.
  Each sequence is equally weighted; adaptive multitaper is not (yet) supported.
  If `window` is specified, each column is applied as a taper. The sum of
  periodograms is normalized by the total sum of squares of `window`.

* `:autocorrelation` or `:ac`: Use the autocorrelation function (AC). The value
  where the AC first comes back close to 1 is the period of the signal. The
  keyword `L = length(v)÷10` denotes the length of the AC (thus, given the
  default setting, this method will fail if there less than 10 periods in the
  signal). The keyword `ϵ = 0.2` (`\\epsilon`) means that `1-ϵ` counts as "1" for the AC.

## Methods not requiring evenly sampled data

These methods tend to be slow, but versatile and low-error.

* `:lombscargle` or `:ls`: Use the Lomb-Scargle algorithm to compute a
  periodogram.  The advantage of the Lomb-Scargle method is that it does not
  require an equally sampled dataset and performs well on undersampled datasets.
  Constraints have been set on the period, since Lomb-Scargle tends to have
  false peaks at very low frequencies.  That being said, it's a very flexible
  method.  It is extremely customizable, and the keyword arguments that can be
  passed to it are given [in the documentation](https://juliaastro.github.io/LombScargle.jl/stable/index.html#LombScargle.plan).

*  `:zerocrossing` or `:zc`: Find the zero crossings of the data, and use the
  average difference between zero crossings as the period.  This is a naïve
  implementation, with only linear interpolation; however, it's useful as a
  sanity check.  The keyword `line` controls where the "crossing point" is.
  It deffaults to `mean(v)`.

For more information on the periodogram methods, see the documentation of
DSP.jl and LombScargle.jl.
"""
function estimate_period(v, method, t = 0:length(v)-1; kwargs...)
    @assert length(v) == length(t)

    even_methods  = (:periodogram, :pg, :multitaper, :mt, :autocorrelation, :ac)
    other_methods = (:lombscargle, :ls, :zerocrossing, :zc)
    if method ∉ even_methods && method ∉ other_methods
        error("Unknown method (`$method`) given to `estimate_period`.")
    end

    period = if method ∈ even_methods
                isevenlysampled(t) ||
                error("Your time data was not evenly sampled,
                and the algorithm `$method` requires evenly sampled data.")

                if method == :periodogram || method == :pg
                    _periodogram_period(v, t; kwargs...)
                elseif method == :multitaper || method == :mt
                    _mt_period(v, t; kwargs...)
                elseif method == :autocorrelation || method == :ac
                    _ac_period(v, t; kwargs...)
                end
            else
                if method == :lombscargle || method == :ls
                    _ls_period(v, t; kwargs...)
                elseif method == :zerocrossing || method == :zc
                    _zc_period(v, t; kwargs...)
                end
            end

    return period
end

function isevenlysampled(t::AbstractVector)
    for i in 2:length(t)-1
        if !(t[nextind(t, i)] - t[i] ≈ t[i] - t[prevind(t, i)])
            return false
            break
        end
    end
    return true
end

# AbstractRanges are defined to be evenly sampled.
isevenlysampled(::AbstractRange) = true

################################################################################
#                           Autocorrelation Function                           #
################################################################################
"""
    _ac_period(v, t; ϵ = 0.2, L = length(v)÷10)

Use the autocorrelation function (AC). The value where the AC first
comes back close to 1 is the period of the signal. The keyword
`L = length(v)÷10` denotes the length of the AC (thus, given the default
setting, this method will fail if there less than 10 periods in the signal).
The keyword `ϵ = 0.2` means that `1-ϵ` counts as "1" for the AC.
"""
function _ac_period(v, t; ϵ = 0.2, L = length(v)÷10)
    err = "The autocorrelation did not become close to 1."
    ac = autocor(v, 0:L)
    j = 0
    local_maxima = findall(i -> ac[i-1] < ac[i] ≥ ac[i+1], 2:length(ac)-1)
    isempty(local_maxima) && error("AC did not have any local maxima")
    # progressively scan the local maxima and find the first within 1-ϵ
    for i in local_maxima
        if ac[i] > 1-ϵ
            j = i
            break
        end
    end
    j == 0 && error("No local maximum of the AC exceeded 1-ϵ")
    # since now it holds that ac[j] is a local maximum within 1-ϵ:
    period = t[j+1] - t[1]
    return period
end

################################################################################
#                                 Periodograms                                 #
################################################################################

########################################
#          Basic Periodogram           #
########################################

"""
    _periodogram_period(v, t; kwargs...)

Use the fast Fourier transform to compute a periodogram (power-spectrum) of the
given data.  Data must be evenly sampled.
"""
function _periodogram_period(v, t; kwargs...)

    p = Periodograms.periodogram(v; fs = length(t)/(t[end] - t[1]), kwargs...)

    return 1 / Periodograms.freq(p)[findmax(Periodograms.power(p))[2]]

end

########################################
#        Multitaper periodogram        #
########################################

"""
    _mt_period(v, t;
                    nw = 4, ntapers = DSP.ceil(2nw)-1,
                    window = DSP.dpss(length(s), nw, ntapers), kwargs...)

The multitaper method reduces estimation bias by obtaining multiple independent
estimates from the same sample. Data tapers are then windowed and the power
spectra are obtained.
`nw` is the time-bandwidth product, and `ntapers` is the number of tapers.
If `window` is not specified, the signal is tapered with `ntapers` discrete
prolate spheroidal sequences with time-bandwidth product `nw`.
Each sequence is equally weighted; adaptive multitaper is not (yet) supported.
If `window` is specified, each column is applied as a taper. The sum of
periodograms is normalized by the total sum of squares of `window`.
"""
function _mt_period(v, t;
                        nw = 4,
                        ntapers::Integer = ceil(Int, 2nw)-1,
                        window = DSP.dpss(length(v), nw, ntapers),
                        kwargs...
                    )

    Sys.WORD_SIZE == 32 &&
        error("multitaper method doesn't work on 32-bit systems. Sorry!")

    p = Periodograms.mt_pgram(v;
                                 fs = length(t)/(t[end] - t[1]),
                                 nw = nw,
                                 ntapers = ntapers,
                                 window = window,
                                 kwargs...
                             )

    return 1 / Periodograms.freq(p)[findmax(Periodograms.power(p))[2]]
end

########################################
#       Lomb-Scargle periodogram       #
########################################

"""
    _ls_period(v, t;
                minimum_period::Real = 2 * (t[end] - t[1]) / length(t),
                maximum_period::Real = (t[end] - t[1]) / 1.5,
                kwargs... # see the documentation of LombScargle.jl for these
        )

Uses the Lomb-Scargle algorithm to compute a periodogram.  The advantage of the Lomb-Scargle method is that it does not require an equally sampled dataset, and it performs well on undersampled datasets.
Constraints have been set on the period, since Lomb-Scargle tends to have false peaks at very low frequencies.  That being said, it's a very flexible method.  It is extremely customizable, and the keyword arguments that can be passed to it are given [in the documentation](https://juliaastro.github.io/LombScargle.jl/stable/index.html#LombScargle.plan).
"""
function _ls_period(v, t;
                    minimum_period::Real = 4 * (t[end] - t[1]) / length(t),
                    maximum_period::Real = (t[end] - t[1]) / 1.5,
                    kwargs...
            )

    plan = LombScargle.plan(t, v;
                            minimum_frequency = 1/maximum_period,
                            maximum_frequency = 1/minimum_period,
                            kwargs...
                        )

    p = LombScargle.lombscargle(plan)
    return LombScargle.findmaxperiod(p)[1]
end

################################################################################
#                                Zero crossings                                #
################################################################################
"""
    _zc_period(v, t; line = mean(v))

Find the zero crossings of the data, and use the average difference between
zero crossings as the period.  This is a naïve implementation, with no
interpolation; however, it's useful as a sanity check.

The keyword `line` controls where the "zero point" is.

The implementation of the function was inspired by [this gist](https://gist.github.com/endolith/255291),
and has been modified for performance and to support arbitrary time grids.
"""
function _zc_period(v, t; line = mean(v))
    # This line might be a little opaque, so I'll provide some more detail.
    # The macro @. applies the broadcast operator to all operations in the
    # expression given to it, except for those which have an $ in front of
    # them.  It also does some optimization to fuse the broadcasts into a
    # single operation.
    # So, @. a * b + c is more like fma.(a, b, c) than a .* b .+ c.
    # The `@views` macro provides a view into an array without actually copying it,
    # so this allocates less.
    # What this is doing is finding all rising line crossings, which means all
    # pairs of consecutive points for which the first is below the line and the
    # second is above it.
    inds = findall(@. ≥(line, $@view(v[2:end])) & <(line, $@view(v[1:end-1])))
    mean(diff(t[inds]))
end


text:
module algfxn
    using GeometryTypes
    using Distributions
    using algT

    function dist(pt1, pt2)
        x2,y2 = pt2[1], pt2[2]
        x1,y1 = pt1[1], pt1[2]
        dist = sqrt( (x1-x2)^2 + (y1-y2)^2 )
        return dist
#        return min_euclidean(pt1), pt2)) #TODO fix this
    end

    function sampleFree(roomW, roomH, obstacles)
        pt_rand = Void
        sampledFree = false

        while !sampledFree
            xrand = rand(Uniform(0, roomW))
            yrand = rand(Uniform(0, roomH))
            pt = Point(xrand, yrand)
            if algfxn.isFreeState(pt, obstacles) # this does NOT check if point is on a line (e.g. a wall)
                pt_rand = Point(xrand, yrand)
                sampledFree = true
            end
        end

        return pt_rand
    end

    function steer(fromPt, toPt, connectRadius)
        dist = algfxn.dist(fromPt, toPt)
        # If pt_rand is too far from its nearest node, "truncate" it to be closer
        if dist <= connectRadius
            return toPt
        else
            x1,y1 = fromPt
            x2,y2 = toPt

            theta = atan2((y2-y1) , (x2-x1))
            newX = x1 + connectRadius * cos(theta)
            newY = y1 + connectRadius * sin(theta)

            newPt = Point(newX, newY)
            return newPt
        end
    end

    function decompRect(r::HyperRectangle) #GeometryTypes.HyperRectangle{2,Float64}
        corners = decompose(Point{2, Float64}, r)
        corners = [Point(pt) for pt in corners]
        lineBottom  = LineSegment(corners[1], corners[2])
        lineTop     = LineSegment(corners[3], corners[4])
        lineLeft    = LineSegment(corners[1], corners[3])
        lineRight   = LineSegment(corners[2], corners[4])
        lines = Vector{LineSegment}()
        push!(lines, lineTop, lineRight, lineBottom, lineLeft)
        return lines
    end

    function isFreeState(node::Point{2, Float64}, obsList::Vector{HyperRectangle})
        for obs in obsList
            if contains(obs, node)
                return false
            end
        end
        return true
    end

    function isFreeMotion(line::LineSegment, obstacles::Vector{HyperRectangle}, walls::Vector{LineSegment})
        return (!algfxn.isCollidingWalls(line, walls) && !algfxn.isCollidingObstacles(line, obstacles))
    end

    #helper1
    function isCollidingObstacles(line::LineSegment, obstacles::Vector{HyperRectangle})
        for obs in obstacles
            rectLines = decompRect(obs)
            for rectline in rectLines
                if intersects(line, rectline)[1]
                    return true
                end
            end
        end
        return false
    end

    #helper2
    function isCollidingWalls(line::LineSegment, walls::Vector{LineSegment})
        for wall in walls
            if intersects(line, wall)[1]
                return true
            end
        end
        return false
    end

#=     function findNearestNodes(nodestate, nodeslist, maxDist) =#
        #  # given maxDist, return all nodes within that distance of node
        #  # list is NOT sorted by distance
        #  nearestNodes = Vector{Tuple{algT.TreeNode, Float64}}()
        #  for n in nodeslist
            #  dist = algfxn.dist(nodestate, n.state)
            #  if dist < maxDist
                #  push!(nearestNodes, (n, dist))
            #  end
        #  end
        # return nearest
    #  end

    function nearestN(pt::Point{2,Float64}, nodeslist)
        #
        # This loops through all nodes  TODO
        nearestDist = 99999999;
        nearestNode = [];

        for n in nodeslist
            dist = algfxn.dist(pt, n.state)
            if dist < nearestDist
                nearestDist = dist
                nearestNode = n
            end
        end
        #print("$(nearestNode.state), which is type $(typeof(nearestNode)), and from $pt   ")
        return nearestNode
    end

    function inGoalRegion(pt, pt_goal)
        if pt == pt_goal
            return true
        else
            dist = algfxn.dist(pt, pt_goal)
            if dist <= 3 #HARDCODED REGION #TODO
                return true
            end
        end

        return false
    end

    function costPath(solPath)
        pathcost = 0
        for i in 2:length(solPath)
            curN  = solPath[i].state
            prevN = solPath[i-1].state
            pathcost += algfxn.dist(curN, prevN)
        end
        return pathcost
    end

    function findNode(nodeID, nodeslist)
        for n in nodeslist #TODO I reverted back to looping through alll the nodes; fix this
            if n.id == nodeID
                return n
                break
            end
        end
        #idlist = [node.id for node in nodeslist]
        #index = findfirst(idlist, nodeID)
        #return nodeslist[index]
    end

    ######## RRT* Specific Functions ##########

    function chooseParent(nearestN, newNode,  nodeslist, connectRadius, obstacles, walls)
                        #print("\n --- looking for parent ---\n")
        for n in nodeslist
            pt = n.state
            if ( isFreeMotion( LineSegment(pt, nearestN.state) ,  obstacles, walls) &&
                dist(pt, newNode.state) < connectRadius  &&
                n.cost + dist(pt, newNode.state) < nearestN.cost + dist(nearestN.state, pt)  )

                nearestN = n
            end
        end
        newNode.cost = nearestN.cost + dist(nearestN.state, newNode.state)
        newNode.parentID = nearestN.id
        return newNode, nearestN
    end

    function rewire(nodeslist, edgeslist, newNode, connectRadius, obstacles, walls)
        for j in length(nodeslist)
            p = nodeslist[j]
            newMove = LineSegment( Point(p.state), Point(newNode.state))

            if ( algfxn.isFreeMotion(newMove, obstacles, walls) &&
				p.id != newNode.parentID && dist(p.state, newNode.state) < connectRadius &&
                newNode.cost + dist(p.state, newNode.state) < p.cost )
                ###terriblecode just loops for now
                # remove edge

				#  i = findfirst(edgeslist) do y
                    #  y.startNode.id == p.parentID && y.endNode.id == p.id
				#  end

                # print("\n --- REWIRING ---\n")
                for i in length(edgeslist)
                    e = edgeslist[i]
                    if (e.startNode.id == p.parentID && e.endNode.id == p.id)
                        deleteat!(edgeslist, i)
                        # print("\n deletd an edge! $(e) \n")
                        break
                    end

                end

                p.parentID = newNode.id
                p.cost = newNode.cost + dist(p.state, newNode.state)
                nodeslist[i] = p #update node
                foo = algT.Edge(p, newNode)
                # add back new edge
				push!(edgeslist,foo)
            end
        end
        return nodeslist, edgeslist
    end

    #  function findEdge(startID, endID, edgeslist)
                #  #removeID = findEdge(p.startID, p.endID)
        #  for edge in edgeslist
            #  if edge.startID == startID && edge.endID == endID
                #  return edge
            #  end
        #  end
        #  return 0
    #  end

    ##################

end


text:
julia> using Gtk

julia> screen_size()
(3840, 1080)

julia>

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("gen.jl")

ENV["TEST_SRVR_ASYNC"] = "true"
include("srvr.jl")
include("clnt.jl")

include("memtransport_tests.jl")
include("filetransport_tests.jl")
include("utils_tests.jl")


text:
using HypergraphsEpidemics
using DataStructures
using SimpleHypergraphs
using Dates

households = generate_dataset("resources/gen_settings/Salerno/gen_param.json","resources/datasets/Salerno/")

# mean = sum(h -> h.num_components, households) / length(households)

start_simulation(
    "resources/datasets/Salerno/dataset.csv",
    "resources/datasets/Salerno/contacts.csv",
    "resources/gen_settings/Salerno/gen_param.json",
    "resources/datasets/Salerno/info.json"
)

df_generated_model, intervals, user2vertex, loc2he = HypergraphsEpidemics.generate_model_data(
                                            "resources/datasets/Salerno/contacts.csv",
                                            [:Time, :Id, :position], #Column list
                                            :Id, #userid column
                                            :position, #venue column
                                            :Time, #Check in column
                                            "yyyy-mm-ddTHH:MM:SS", #Date format
                                            Δ = Dates.Millisecond(86400000), # 24 hours
                                            δ = Dates.Millisecond(900000), # minutes
                                            datarow = 2 # Data starts at second row
                                            # limit = 10000 # limit
    )

length(keys(loc2he))

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
using DataStructures

function solve()
    n, m = [parse(Int, x) for x in split(readline())]
    ds = IntDisjointSets(n)
    for _=1:m
        a, b = [parse(Int, x) for x in split(readline())]
        union!(ds, a, b)
    end

    size = zeros(Int, n)
    for x=1:n
        root = find_root!(ds, x)
        size[root] += 1
    end
    maximum(size)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
"""

Particles constructor for DynamicHMC
Convert DynamcHMC samples to a Particle NamedTuple
* `posterior`: an array of NamedTuple consisting of mcmc samples
"""
function as_particles(posterior)
    d = Dict{Symbol, Union{Particles, Vector{Particles}}}()
    pnt = posterior[1]
    for parm in keys(pnt)
        if size(pnt[parm], 1) > 1
            d[parm] = Particles[]
            for i in 1:size(pnt[parm], 1)
                temp = Float64[]
                for post in posterior
                    push!(temp, post[parm][i])
                end
                push!(d[parm], Particles(temp))
            end
        else
            temp = Float64[]
            for post in posterior
                push!(temp, post[parm])
            end
            d[parm] = Particles(temp)
        end
    end
    return (; d...)
end

function nptoa3d(posterior)
    Np = length(vcat(posterior[1]...))
    Ns = length(posterior)
    a3d = Array{Float64,3}(undef,Ns,Np,1)
    for (i,post) in enumerate(posterior)
        temp = Float64[]
        for p in post
            push!(temp,values(p)...)
        end
        a3d[i,:,1] = temp'
    end
    parameter_names = getnames(posterior)
    return (a3d, parameter_names)
end

function getnames(post)
    nt = post[1]
    Np =length(vcat(nt...))
    parm_names = fill("",Np)
    cnt = 0
    for (k,v) in pairs(nt)
        N = length(v)
        if isa(v,Array)
            for i in 1:N
                cnt += 1
                parm_names[cnt] = string(k,"[",i,"]")
            end
        else
            cnt+=1
            parm_names[cnt] = string(k)
        end
    end
    return Symbol.(parm_names)
end

export
    Particles,
    as_particles


text:
#######################
# The Bernstein basis
#######################

using BasisFunctions: UnitInterval, PolynomialBasis
using BasisFunctions: hasderivative, hasantiderivative, ordering, support

const BernsteinInterval = UnitInterval
const BernsteinIndex = NativeIndex{:bernstein}

get_bernstein_nodes(T,n) = collect(LinRange{T}(0.0, 1.0, n))

struct Bernstein{T} <: PolynomialBasis{T}
    n :: Int
    nodes :: ScatteredGrid{T}

    function Bernstein(nodes::ScatteredGrid{T}) where {T}
        new{T}(length(nodes.points), nodes)
    end
    function Bernstein(ξ::Vector{T}) where {T}
        Bernstein(ScatteredGrid(ξ, BernsteinInterval{T}()))
    end
end

Bernstein(n::Int) = Bernstein(get_bernstein_nodes(Float64,n))
Bernstein{T}(n::Int) where {T} = Bernstein(get_bernstein_nodes(T,n))

# Convenience constructor: map the Bernstein basis to the interval [a,b]
Bernstein(ξ, a::Number, b::Number) = rescale(Bernstein(ξ), a, b)

nodes(b::Bernstein)  = b.nodes.points
nnodes(b::Bernstein) = b.n
degree(b::Bernstein) = b.n-1

BasisFunctions.native_index(b::Bernstein, idx) = BernsteinIndex(idx)
BasisFunctions.linear_index(b::Bernstein, idx) = BernsteinIndex(idx)
BasisFunctions.ordering(b::Bernstein) = Base.OneTo(nnodes(b))
BasisFunctions.support(b::Bernstein{T}) where {T} = BernsteinInterval{T}()

BasisFunctions.interpolation_grid(b::Bernstein{T}) where {T} = ScatteredGrid(get_bernstein_nodes(T,length(b)), BernsteinInterval{T}())

BasisFunctions.similar(::Bernstein, ::Type{T}, n::Int) where {T} = Bernstein{T}(n)
BasisFunctions.similar(::Bernstein, ::Type{T}, ξ::Vector{T}) where {T} = Bernstein(ξ)

Base.size(b::Bernstein) = (b.n,)

function _bernstein(b::Bernstein{T}, i::BernsteinIndex, n::BernsteinIndex, x::T) where {T}
    if i < 0 || i > n
        return zero(T)
    else
        if n == 0
            return one(T)
        else
            return _bernstein(b, i, n-1, x) * (1-x) + _bernstein(b, i-1, n-1, x) * x
        end
    end
end

_bernstein(b::Bernstein, i, n, x) = _bernstein(b, native_index(b,i), native_index(b,n), x)

function BasisFunctions.unsafe_eval_element(b::Bernstein, i, x)
    _bernstein(b, i-1, nnodes(b)-1, x)
end

function BasisFunctions.unsafe_eval_element_derivative(b::Bernstein, i, x)
    (nnodes(b)-1) * ( _bernstein(b, i-2, nnodes(b)-2, x) - _bernstein(b, i-1, nnodes(b)-2, x) )
end


text:
"""
    GaussianCopula{d,MT}

Fields:
  - Σ::MT - covariance matrix

Constructor

    GaussianCopula(Σ)

The [Gaussian Copula](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Gaussian_copula) is the
copula of a [Multivariate normal distribution](http://en.wikipedia.org/wiki/Multivariate_normal_distribution). It is constructed as :

```math
C(\\mathbf{x}; \\boldsymbol{\\Sigma}) = F_{\\Sigma}(F_{\\Sigma,i}^{-1}(x_i),i\\in 1,...d)
```
where ``F_{\\Sigma}`` is a cdf of a gaussina random vector and `F_{\\Sigma,i}` is the ith marignal cdf, while ```\\Sigma`` is the covariance matrix.

It can be constructed in Julia via:
```julia
C = GaussianCopula(Σ)
```

The random number generation works as expected:
```julia
rand(C,1000)
# or
Random.rand!(C,u)
```

And yo can fit the distribution via :
```julia
fit(GaussianCopula,data)
```
"""
struct GaussianCopula{d,MT} <: EllipticalCopula{d,MT}
    Σ::MT
    function GaussianCopula(Σ)
        make_cor!(Σ)
        return new{size(Σ,1),typeof(Σ)}(Σ)
    end
end
U(::Type{T}) where T<: GaussianCopula = Distributions.Normal
N(::Type{T}) where T<: GaussianCopula = Distributions.MvNormal
function Distributions.fit(::Type{CT},u) where {CT<:GaussianCopula}
    dd = Distributions.fit(N(CT), quantile.(U(CT)(),u))
    Σ = Matrix(dd.Σ)
    return GaussianCopula(Σ)
end


text:
function sdpblockpreconditoner( H::HermitianBlockMatrix )
    model = JuMP.Model( COSMO.Optimizer )
    JuMP.set_silent( model )

    bs = blocksize( H )

    if( ( bs[1] != bs[2] ) || ( size( H, 1 ) != size( H, 2 ) ) )
        throw( DimensionMismatch( "Hessian H must be square" ) )
    end

    N  = bs[1]
    nu = size( H[Block(1, 1)], 1 )

    # The cholesky of the Hessian is needed for the constraints
    C = cholesky( Matrix(H) )

    E = zeros( GenericAffExpr{Float64,VariableRef}, N*nu, N*nu )

    for i in 1:N
        # Create the diagonal indices
        ind = ( ( nu*(i-1)+1 ):(nu*i) )

        # Create the blocks that go on the diagonal of the constraint matrix
        E[ind, ind] = @variable( model, [j=1:nu, k=1:nu], base_name="e$i" )
    end

    # Initialize the objective
    @variable( model, t )
    @objective( model, Min, t )

    # E should be PSD
    @constraint( model, E in PSDCone() )

    # The coefficient of E is the smallest eigenvalue of the preconditoned matrix
    @SDconstraint( model, Matrix(H) - E >= 0 )
    @SDconstraint( model, [E           Matrix(C.L);
                           Matrix(C.U) t*I(N*nu)] >= 0 )

    JuMP.optimize!( model )

    if( termination_status( model ) == MOI.OPTIMAL )
        sol = value.( E )
        obj = objective_value( model )

    elseif( termination_status( model ) == MOI.TIME_LIMIT && has_values( model ) )
        sol = value.( E )
        obj = objective_value( model )

    else
        # If we can't compute it using a pure diagonal structure, fall back on a full E matrix
        @warn "Unable to compute a preconditioning matrix using pure diagonal structure. Trying a full matrix."

        fullModel = JuMP.Model( COSMO.Optimizer )
        JuMP.set_silent( fullModel )

        # Initialize the objective
        @variable( fullModel, t )
        @objective( fullModel, Min, t )

        # E should be PSD
        @variable( fullModel, E[1:N*nu, 1:N*nu] )
        @constraint( fullModel, E in PSDCone() )

        # The coefficient of E is the smallest eigenvalue of the preconditoned matrix
        @SDconstraint( fullModel, Matrix(H) - E >= 0 )
        @SDconstraint( fullModel, [E           Matrix(C.L);
                                   Matrix(C.U) t*I(N*nu)] >= 0 )

        JuMP.optimize!( fullModel )

        if( termination_status( fullModel ) == MOI.OPTIMAL )
            sol = value.( E )
            obj = objective_value( fullModel )

        elseif( termination_status( fullModel ) == MOI.TIME_LIMIT && has_values( fullModel ) )
            sol = value.( E )
            obj = objective_value( fullModel )
        else
            error( "Unable to compute preconditioner using full matrix" )
        end
    end

    blockSol = BlockMatrix( sol, [nu for i=1:N], [nu for i=1:N] )
    precond  = BlockMatrix( zeros( N*nu, N*nu ), [nu for i=1:N], [nu for i=1:N] )

    for i = 1:N
        eig = eigen( blockSol[Block(i,i )] )

        # We need the square root of the inverse eigenvalues
        v = sqrt.( 1 ./ eig.values )

        # Reassemble the block
        precond[Block(i, i)] = eig.vectors * diagm( v ) * eig.vectors'
    end

    return precond
end

"""
    circulantblockpreconditioner( clqr::ConstrainedTimeInvariantLQR; variant::Symbol = :Strang, form::Symbol = :Symmetric, mat::Symbol = :Full )

Compute the circulant block preconditioner of a constrained LTI LQR problem. This preconditioner is similar
to the one computed using the SDP approach in [`sdpblockpreconditoner`](@ref), but is simpler to compute and
only works for LTI problems.

This will return the diagonal block used to form the preconditioner, not the entire preconditioning matrix.
By default, the block for the symmetric preconditioner ``L^{-1} H L^{-T}`` is returned, but the left/right
preconditioning block can be requested by setting `precondType` to either `:Left` or `:Right`.

"""
function circulantblockpreconditioner( clqr::ConstrainedTimeInvariantLQR; variant::Symbol = :Strang, form::Symbol = :Symmetric, mat::Symbol = :Full )
    sys = getsystem( clqr )

    # We need the solution to the DARE for the prestabilized system (using the original Q matrix)
    P, = ared( sys.A, sys.B, clqr.R, clqr.Q)

    # Form the Hessian block
    H = sys.B' * P * sys.B + clqr.R

    if variant == :Strang
        M = H
    elseif variant == :Chan
        E = dftmatrix( size( H, 1 ), unitary = true )
        D = Diagonal( E * H * E' )
        M = real( E' * D * E )
    else
        error( "Unknown variant specified. Valid types are :Strang and :Chan" )
    end

    # If a symmetric preconditioner is requested, use the lower-diagonal term of the Cholesky decomposition
    L = ( form == :Symmetric ) ? cholesky( Hermitian( M ) ).L : M

    if( mat == :Full)
        bs = ( ninputs(sys) for i=1:clqr.N )
        return BlockMatrix( kron( I(clqr.N), L ), collect( bs ), collect( bs ) )
    else
        return L
    end
end


text:
using ITensorGaussianMPS
using LinearAlgebra
using ITensors

# Half filling
N = 50
Nf = N ÷ 2

@show N, Nf

# Maximum MPS link dimension
_maxlinkdim = 100

@show _maxlinkdim

# DMRG cutoff
_cutoff = 1e-12

# Hopping
t = 1.0

# Electron-electron on-site interaction
U = 1.0

@show t, U

# Free fermion Hamiltonian
ampo = AutoMPO()
for n in 1:(N - 1)
  ampo .+= -t, "Cdag", n, "C", n + 1
  ampo .+= -t, "Cdag", n + 1, "C", n
end

# Hopping Hamiltonian with N spinless fermions
h = hopping_hamiltonian(ampo)

# Get the Slater determinant
Φ = slater_determinant_matrix(h, Nf)

# Create an mps for the free fermion ground state
s = siteinds("Fermion", N; conserve_qns=true)
println("Making free fermion starting MPS")
@time ψ0 = slater_determinant_to_mps(
  s, Φ; eigval_cutoff=1e-4, cutoff=_cutoff, maxdim=_maxlinkdim
)
@show maxlinkdim(ψ0)

# Make an interacting Hamiltonian
for n in 1:(N - 1)
  ampo .+= U, "N", n, "N", n + 1
end
H = MPO(ampo, s)

# Random starting state
ψr = randomMPS(s, n -> n ≤ Nf ? "1" : "0")

println("\nRandom state starting energy")
@show flux(ψr)
@show inner(ψr, H, ψr)

println("\nFree fermion starting energy")
@show flux(ψ0)
@show inner(ψ0, H, ψ0)

println("\nRun dmrg with random starting state")
sweeps = Sweeps(20)
maxdim!(sweeps, 10, 20, 40, _maxlinkdim)
cutoff!(sweeps, _cutoff)
@time dmrg(H, ψr, sweeps)

println("\nRun dmrg with free fermion starting state")
sweeps = Sweeps(4)
maxdim!(sweeps, _maxlinkdim)
cutoff!(sweeps, _cutoff)
@time dmrg(H, ψ0, sweeps)

nothing