text:
#=
    Sorting networks that process 1,2,..16 values are given.
    Sorting networks that process 17..25 values are given in "swapsortr_17to25.jl".

    These sorting networks should minimize conditional swaps.
    The first sixteen are known to minimize conditional swaps.

    Functions are written with groups of internally parallelizable
    statements given together and sequential steps are separated.

    `swapsortr(x1::T, .., xn::T)::NTuple{n,T}`
    `swapsortr(tup::NTuple{n,T})::NTuple{n,T}`
    `swapsortr(vec::AbstractVector{T})::NTuple{N,T}`

    networks were selected using software by John Gamble
        http://pages.ripco.net/~jgamble/nw.html
        http://search.cpan.org/dist/Algorithm-Networksort/
=#

const VALS = ([Val{i} for i=1:ITEMS_MAX]...,)

@inline function swapsortr(vec::AbstractVector{T}) where {T}
    n = length(vec)
    n > ITEMS_MAX && throw(ErrorException("swapsortr(vector) not defined where length(vector) > $ITEMS_MAX"))
    return swapsortr(VALS[length(vec)], vec)
end

#    sort 1 value with 0 minmaxs in 1 stage

function swapsortr(a::T) where {T}
    return (a,)
end

@inline swapsortr(x::NTuple{1,T}) where {T} = swapsortr(x[1])
@inline swapsortr(::Type{Val{1}}, x::AbstractVector{T}) where {T} = swapsortr(x[1])

#    sort 2 values with 1 minmaxs in 1 stage

function swapsortr(a::T, b::T) where {T}

    a, b = max_min(a, b)

    return a,b
end

@inline swapsortr(x::NTuple{2,T}) where {T} = swapsortr(x[1], x[2])
@inline swapsortr(::Type{Val{2}}, x::AbstractVector{T}) where {T} = swapsortr(x[1], x[2])

#    sort 3 values with 3 minmaxs in 3 parallel stages

function swapsortr(a::T, b::T, c::T) where {T}

    b, c = max_min(b, c)

    a, c = max_min(a, c)

    a, b = max_min(a, b)

    return a, b, c
end

@inline swapsortr(x::NTuple{3,T}) where {T} = swapsortr(x[1], x[2], x[3])
@inline swapsortr(::Type{Val{3}}, x::AbstractVector{T}) where {T} = swapsortr(x[1], x[2], x[3])

#    sort 4 values with 5 minmaxs in 3 parallel stages

function swapsortr(a::T, b::T, c::T, d::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)

    a, c = max_min(a, c)
    b, d = max_min(b, d)

    b, c = max_min(b, c)

    return a, b, c, d
end

@inline swapsortr(x::NTuple{4,T}) where {T} = swapsortr(x[1], x[2], x[3], x[4])
@inline swapsortr(::Type{Val{4}}, x::AbstractVector{T}) where {T} = swapsortr(x[1], x[2], x[3], x[4])

#    sort 5 values with 9 minmaxs in 5 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T) where {T}

    b, c = max_min(b, c)
    d, e = max_min(d, e)

    b, d = max_min(b, d)
    a, c = max_min(a, c)

    c, e = max_min(c, e)
    a, d = max_min(a, d)

    a, b = max_min(a, b)
    c, d = max_min(c, d)

    b, c = max_min(b, c)

    return a, b, c, d, e
end

@inline swapsortr(x::NTuple{5,T}) where {T} = swapsortr(x[1], x[2], x[3], x[4], x[5])
@inline swapsortr(::Type{Val{5}}, x::AbstractVector{T}) where {T} = swapsortr(x[1], x[2], x[3], x[4], x[5])

#    sort 6 values with 12 minmaxs in 6 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T) where {T}

    b, c = max_min(b, c)
    e, f = max_min(e, f)

    a, c = max_min(a, c)
    d, f = max_min(d, f)

    a, b = max_min(a, b)
    d, e = max_min(d, e)
    c, f = max_min(c, f)

    a, d = max_min(a, d)
    b, e = max_min(b, e)

    c, e = max_min(c, e)
    b, d = max_min(b, d)

    c, d = max_min(c, d)

    return a, b, c, d, e, f
end

@inline swapsortr(x::NTuple{6,T}) where {T} = swapsortr(x[1], x[2], x[3], x[4], x[5], x[6])
@inline swapsortr(::Type{Val{6}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6])

#    sort 7 values with 16 minmaxs in 6 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T) where {T}

    a, e = max_min(a, e)
    b, f = max_min(b, f)
    c, g = max_min(c, g)

    a, c = max_min(a, c)
    b, d = max_min(b, d)
    e, g = max_min(e, g)

    c, e = max_min(c, e)
    d, f = max_min(d, f)
    a, b = max_min(a, b)

    c, d = max_min(c, d)
    e, f = max_min(e, f)

    b, e = max_min(b, e)
    d, g = max_min(d, g)

    b, c = max_min(b, c)
    d, e = max_min(d, e)
    f, g = max_min(f, g)

    return a, b, c, d, e, f, g
end

@inline swapsortr(x::NTuple{7,T}) where {T} = swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7])
@inline swapsortr(::Type{Val{7}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7])

#    sort 8 values with 19 minmaxs in 7 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)

    a, c = max_min(a, c)
    b, d = max_min(b, d)
    e, g = max_min(e, g)
    f, h = max_min(f, h)

    b, c = max_min(b, c)
    f, g = max_min(f, g)
    a, e = max_min(a, e)
    d, h = max_min(d, h)

    b, f = max_min(b, f)
    c, g = max_min(c, g)

    b, e = max_min(b, e)
    d, g = max_min(d, g)

    c, e = max_min(c, e)
    d, f = max_min(d, f)

    d, e = max_min(d, e)

    return a, b, c, d, e, f, g, h
end

@inline swapsortr(x::NTuple{8,T}) where {T} = swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])
@inline swapsortr(::Type{Val{8}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])

#    sort 9 values with 25 minmaxs in 9 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T) where {T}

    a, b = max_min(a, b)
    d, e = max_min(d, e)
    g, h = max_min(g, h)

    b, c = max_min(b, c)
    e, f = max_min(e, f)
    h, i = max_min(h, i)

    a, b = max_min(a, b)
    d, e = max_min(d, e)
    g, h = max_min(g, h)
    c, f = max_min(c, f)

    a, d = max_min(a, d)
    b, e = max_min(b, e)
    f, i = max_min(f, i)

    d, g = max_min(d, g)
    e, h = max_min(e, h)
    c, f = max_min(c, f)

    a, d = max_min(a, d)
    b, e = max_min(b, e)
    f, h = max_min(f, h)
    c, g = max_min(c, g)

    b, d = max_min(b, d)
    e, g = max_min(e, g)

    c, e = max_min(c, e)
    f, g = max_min(f, g)

    c, d = max_min(c, d)

    return a, b, c, d, e, f, g, h, i
end

@inline swapsortr(x::NTuple{9,T}) where {T} = swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])
@inline swapsortr(::Type{Val{9}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])

#    sort 10 values with 29 minmaxs in 9 parallel stages

function swapsortr(a::T,b::T,c::T,d::T,e::T,
              f::T,g::T,h::T,i::T,j::T) where {T}

    e, j = max_min(e, j)
    d, i = max_min(d, i)
    c, h = max_min(c, h)
    b, g = max_min(b, g)
    a, f = max_min(a, f)

    b, e = max_min(b, e)
    g, j = max_min(g, j)
    a, d = max_min(a, d)
    f, i = max_min(f, i)

    a, c = max_min(a, c)
    d, g = max_min(d, g)
    h, j = max_min(h, j)

    a, b = max_min(a, b)
    c, e = max_min(c, e)
    f, h = max_min(f, h)
    i, j = max_min(i, j)

    b, c = max_min(b, c)
    e, g = max_min(e, g)
    h, i = max_min(h, i)
    d, f = max_min(d, f)

    c, f = max_min(c, f)
    g, i = max_min(g, i)
    b, d = max_min(b, d)
    e, h = max_min(e, h)

    c, d = max_min(c, d)
    g, h = max_min(g, h)

    d, e = max_min(d, e)
    f, g = max_min(f, g)

    e, f = max_min(e, f)

    return a,b,c,d,e,f,g,h,i,j
end

@inline swapsortr(x::NTuple{10,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
@inline swapsortr(::Type{Val{10}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])

#    sort 11 values with 35 minmaxs in 9 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)

    b, d = max_min(b, d)
    f, h = max_min(f, h)
    a, c = max_min(a, c)
    e, g = max_min(e, g)
    i, k = max_min(i, k)

    b, c = max_min(b, c)
    f, g = max_min(f, g)
    j, k = max_min(j, k)
    a, e = max_min(a, e)
    d, h = max_min(d, h)

    b, f = max_min(b, f)
    g, k = max_min(g, k)
    e, i = max_min(e, i)

    f, j = max_min(f, j)
    c, g = max_min(c, g)
    a, e = max_min(a, e)
    d, i = max_min(d, i)

    b, f = max_min(b, f)
    g, k = max_min(g, k)
    c, d = max_min(c, d)
    i, j = max_min(i, j)

    b, e = max_min(b, e)
    h, k = max_min(h, k)
    d, f = max_min(d, f)
    g, i = max_min(g, i)

    c, e = max_min(c, e)
    h, j = max_min(h, j)
    f, g = max_min(f, g)

    d, e = max_min(d, e)
    h, i = max_min(h, i)

    return a,b,c,d,e,f,g,h,i,j,k
end

@inline swapsortr(x::NTuple{11,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11])
@inline swapsortr(::Type{Val{11}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11])

#    sort 12 values with 39 minmaxs in 9 parallel stages

function swapsortr(a::T,b::T,c::T,d::T,e::T,
              f::T,g::T,h::T,i::T,j::T,
              k::T,l::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)

    b, d = max_min(b, d)
    f, h = max_min(f, h)
    j, l = max_min(j, l)
    a, c = max_min(a, c)
    e, g = max_min(e, g)
    i, k = max_min(i, k)

    b, c = max_min(b, c)
    f, g = max_min(f, g)
    j, k = max_min(j, k)
    a, e = max_min(a, e)
    h, l = max_min(h, l)

    b, f = max_min(b, f)
    g, k = max_min(g, k)
    d, h = max_min(d, h)
    e, i = max_min(e, i)

    f, j = max_min(f, j)
    c, g = max_min(c, g)
    a, e = max_min(a, e)
    h, l = max_min(h, l)
    d, i = max_min(d, i)

    b, f = max_min(b, f)
    g, k = max_min(g, k)
    c, d = max_min(c, d)
    i, j = max_min(i, j)

    b, e = max_min(b, e)
    h, k = max_min(h, k)
    d, f = max_min(d, f)
    g, i = max_min(g, i)

    c, e = max_min(c, e)
    h, j = max_min(h, j)
    f, g = max_min(f, g)

    d, e = max_min(d, e)
    h, i = max_min(h, i)

    return a,b,c,d,e,f,g,h,i,j,k,l
end

@inline swapsortr(x::NTuple{12,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12])
@inline swapsortr(::Type{Val{12}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12])

#    sort 13 values with 45 minmaxs in 10 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T) where {T}

    b, h = max_min(b, h)
    j, l = max_min(j, l)
    d, e = max_min(d, e)
    f, i = max_min(f, i)
    a, m = max_min(a, m)
    c, g = max_min(c, g)

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, g = max_min(e, g)
    i, l = max_min(i, l)
    h, m = max_min(h, m)
    f, j = max_min(f, j)

    a, c = max_min(a, c)
    d, h = max_min(d, h)
    k, l = max_min(k, l)
    b, e = max_min(b, e)
    g, m = max_min(g, m)

    h, i = max_min(h, i)
    l, m = max_min(l, m)
    e, j = max_min(e, j)
    g, k = max_min(g, k)

    d, e = max_min(d, e)
    f, g = max_min(f, g)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    b, h = max_min(b, h)

    c, g = max_min(c, g)
    j, l = max_min(j, l)
    b, d = max_min(b, d)
    e, h = max_min(e, h)
    i, k = max_min(i, k)
    a, f = max_min(a, f)

    c, f = max_min(c, f)
    g, i = max_min(g, i)
    j, k = max_min(j, k)

    b, c = max_min(b, c)
    d, f = max_min(d, f)
    h, i = max_min(h, i)
    e, g = max_min(e, g)

    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)

    d, e = max_min(d, e)
    f, g = max_min(f, g)

    return a,b,c,d,e,f,g,h,i,j,k,l,m
end

@inline swapsortr(x::NTuple{13,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13])
@inline swapsortr(::Type{Val{13}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13])

#    sort 14 values with 51 minmaxs in 10 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)

    a, c = max_min(a, c)
    e, g = max_min(e, g)
    i, k = max_min(i, k)
    b, d = max_min(b, d)
    f, h = max_min(f, h)
    j, l = max_min(j, l)

    a, e = max_min(a, e)
    i, m = max_min(i, m)
    b, f = max_min(b, f)
    j, n = max_min(j, n)
    c, g = max_min(c, g)
    d, h = max_min(d, h)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)

    f, k = max_min(f, k)
    g, j = max_min(g, j)
    d, m = max_min(d, m)
    h, l = max_min(h, l)
    b, c = max_min(b, c)
    e, i = max_min(e, i)

    b, e = max_min(b, e)
    h, n = max_min(h, n)
    c, i = max_min(c, i)
    f, g = max_min(f, g)
    j, k = max_min(j, k)

    c, e = max_min(c, e)
    l, n = max_min(l, n)
    d, i = max_min(d, i)
    h, m = max_min(h, m)

    g, i = max_min(g, i)
    k, m = max_min(k, m)
    d, f = max_min(d, f)
    h, j = max_min(h, j)

    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)

    g, h = max_min(g, h)
    i, j = max_min(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n
end

@inline swapsortr(x::NTuple{14,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14])
@inline swapsortr(::Type{Val{14}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14])

#    sort 15 values with 56 minmaxs in 10 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)

    a, c = max_min(a, c)
    e, g = max_min(e, g)
    i, k = max_min(i, k)
    m, o = max_min(m, o)
    b, d = max_min(b, d)
    f, h = max_min(f, h)
    j, l = max_min(j, l)

    a, e = max_min(a, e)
    i, m = max_min(i, m)
    b, f = max_min(b, f)
    j, n = max_min(j, n)
    c, g = max_min(c, g)
    k, o = max_min(k, o)
    d, h = max_min(d, h)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)

    f, k = max_min(f, k)
    g, j = max_min(g, j)
    d, m = max_min(d, m)
    n, o = max_min(n, o)
    h, l = max_min(h, l)
    b, c = max_min(b, c)
    e, i = max_min(e, i)

    b, e = max_min(b, e)
    h, n = max_min(h, n)
    c, i = max_min(c, i)
    l, o = max_min(l, o)
    f, g = max_min(f, g)
    j, k = max_min(j, k)

    c, e = max_min(c, e)
    l, n = max_min(l, n)
    d, i = max_min(d, i)
    h, m = max_min(h, m)

    g, i = max_min(g, i)
    k, m = max_min(k, m)
    d, f = max_min(d, f)
    h, j = max_min(h, j)

    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)

    g, h = max_min(g, h)
    i, j = max_min(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o
end

@inline swapsortr(x::NTuple{15,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15])
@inline swapsortr(::Type{Val{15}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15])

#    sort 16 values with 60 minmaxs in 10 parallel stages

function swapsortr(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T) where {T}

    a, b = max_min(a, b)
    c, d = max_min(c, d)
    e, f = max_min(e, f)
    g, h = max_min(g, h)
    i, j = max_min(i, j)
    k, l = max_min(k, l)
    m, n = max_min(m, n)
    o, p = max_min(o, p)

    a, c = max_min(a, c)
    e, g = max_min(e, g)
    i, k = max_min(i, k)
    m, o = max_min(m, o)
    b, d = max_min(b, d)
    f, h = max_min(f, h)
    j, l = max_min(j, l)
    n, p = max_min(n, p)

    a, e = max_min(a, e)
    i, m = max_min(i, m)
    b, f = max_min(b, f)
    j, n = max_min(j, n)
    c, g = max_min(c, g)
    k, o = max_min(k, o)
    d, h = max_min(d, h)
    l, p = max_min(l, p)

    a, i = max_min(a, i)
    b, j = max_min(b, j)
    c, k = max_min(c, k)
    d, l = max_min(d, l)
    e, m = max_min(e, m)
    f, n = max_min(f, n)
    g, o = max_min(g, o)
    h, p = max_min(h, p)

    f, k = max_min(f, k)
    g, j = max_min(g, j)
    d, m = max_min(d, m)
    n, o = max_min(n, o)
    h, l = max_min(h, l)
    b, c = max_min(b, c)
    e, i = max_min(e, i)

    b, e = max_min(b, e)
    h, n = max_min(h, n)
    c, i = max_min(c, i)
    l, o = max_min(l, o)
    f, g = max_min(f, g)
    j, k = max_min(j, k)

    c, e = max_min(c, e)
    l, n = max_min(l, n)
    d, i = max_min(d, i)
    h, m = max_min(h, m)

    g, i = max_min(g, i)
    k, m = max_min(k, m)
    d, f = max_min(d, f)
    h, j = max_min(h, j)

    d, e = max_min(d, e)
    f, g = max_min(f, g)
    h, i = max_min(h, i)
    j, k = max_min(j, k)
    l, m = max_min(l, m)

    g, h = max_min(g, h)
    i, j = max_min(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p
end

@inline swapsortr(x::NTuple{16,T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15], x[16])
@inline swapsortr(::Type{Val{16}}, x::AbstractVector{T}) where {T} =
    swapsortr(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15], x[16])


text:
using Documenter, HoTRG_lemon

makedocs()

deploydocs(
		   deps = Deps.pip("mkdocs", "python-markdown-math"),
		   repo = "github.com/yoyoyoju/HoTRG_lemon.jl.git",
		   julia = "0.5",
		   osname = "osx"
)

ext: jl


text:
load("debug.jl")

module TestCollect
import Base.*
import Debug.*

code = :(function f(n)
    x = 0
    for k=1:n
        x += k
        println("k = ", k)
    end
    x = x*x
end)

println(getdefs(code), '\n')

end


text:
include(joinpath(depsdir, "elkflags.jl"))

struct ElkFlagInfo{T}
    name::Symbol
    default::Union{T, Nothing}  #check again v0.7 Some
    description::String
end
ElkFlagInfo() = ElkFlagInfo{Nothing}(:error, "")
Base.eltype(x::ElkFlagInfo{T}) where T = T

struct ElkControlBlockInfo <: AbstractBlockInfo
    name::Symbol
    flags::Vector{<:ElkFlagInfo}
    description::String
end

const ELK_CONTROLBLOCKS = _ELK_CONTROLBLOCKS()

const ELK_EXECS = ["elk", "elk-omp"]

function elk_flaginfo(flag::Symbol)
	for b in ELK_CONTROLBLOCKS
		for f in b.flags
			if f.name == flag
				return f
			end
		end
	end
end

elk_block_info(name::Symbol) = getfirst(x->x.name == name, ELK_CONTROLBLOCKS)

function elk_block_variable(flag_name::Symbol)
	for b in ELK_CONTROLBLOCKS
		for f in b.flags
			if f.name == flag_name
				return b
			end
		end
	end
end

flagtype(::DFInput{Elk}, flag::Symbol) = eltype(elk_flaginfo(flag))


text:
function TabulatedProvider(arg0::Frame, arg1::List, arg2::jint, arg3::AngularDerivativesFilter)
    return TabulatedProvider((Frame, List, jint, AngularDerivativesFilter), arg0, arg1, arg2, arg3)
end

function get_attitude(obj::TabulatedProvider, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::TabulatedProvider, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
# NB. attributes for builtin symbols are mostly set in protected_symbols.jl

# Check if we are trying to add to a symbol bound to itself. If
# so, we warn and return the unevaluated expression. This avoids
# an infinite evaluation loop. This is what Mma does.
macro checkunbound(mx,x,sv)
    esc( quote
           $sv = symval($x)   # could replace with isbound($x)
           if $sv == $x
             @warn("The symbol " * string($x) * " does not have a value, so it's value cannot be changed")
             setfixed($mx)
             return $mx
          end
    end)
end

protectstr(s) = "Symbol '" * string(s) * "' is proctected. Did you mean `==` ?"

function checkprotect(s::Qsym)
    get_attribute(s,:Protected) &&
    symerror(protectstr(s))
end

function checkprotect(s::SJSym)
    isProtected(s) && symerror(protectstr(symname(s)))
end

checkprotect(mx::Mxpr) = checkprotect(mhead(mx))

function warncheckprotect(s::SJSym)
    if isProtected(s)
        @symwarn(protectstr(symname(s)))
        return false
    else
        return true
    end
end
warncheckprotect(mx::Mxpr) = warncheckprotect(mhead(mx))

### Attributes

@mkapprule Attributes :nargs => 1  """
    Attributes(s)

returns attributes associated with symbol `s`. Builtin symbols have
the attribute `Protected`, and may have others, including `HoldFirst`, `SequenceHold`,
`HoldRest`, `HoldAll`, `ReadProtected`, `Constant`, `Locked`, `Flat`, `Listable`,
`NumericFunction`, `OneIdentity`.
"""

@doap Attributes(s::SymString) = get_attributesList(Symbol(s))

get_attributesList(sj::SJSymbol) = tolistfixed(get_attributes(sj))

### SetAttributes

@mkapprule SetAttributes :nargs => 2  """
    SetAttributes(sym,attr)

add `attr` to the list of attributes for `sym`.

    SetAttributes(list,attr)

add `attr` to the list of attributes for each symbol in `list`.

    SetAttributes(sym,list)

add each attribute in `list`  to the list of attributes for `sym`.
"""

function check_set_attributes(sym::SJSymbol, attr::SJSym)
    checkprotect(sym)
    set_attribute(sym,attr)
    Null
end

set_attributes(sym::SJSymbol, attr::SJSym) = check_set_attributes(sym,attr)

function set_attributes(sym::SJSymbol, attrs::T) where T<:Array
    foreach(a -> check_set_attributes(sym,a), attrs)
end

function set_attributes(syms::T, attr::SJSym) where T<:Array
    foreach(sym -> check_set_attributes(sym,attr), syms)
end

@doap SetAttributes(sym::SJSymbol, attr::SJSym) = set_attributes(sym,attr)
@doap SetAttributes(syms::Mxpr{:List}, attr::SJSym) = set_attributes(margs(syms),attr)
@doap SetAttributes(sym::SJSymbol, attrs::Mxpr{:List}) = set_attributes(sym,margs(attrs))

### ClearAttributes

@sjdoc ClearAttributes """
    ClearAttributes(sym,attr)

removes `attr` from the list of attributes of symbol `sym`.
"""
@mkapprule ClearAttributes :nargs => 2

@doap ClearAttributes(sym,attr) = (unset_attribute(sym,attr); Null)

### Unprotect

###  Wolfram does this: Protect[a,b,c,....] returns a list of symbol names (strings) of those symbols which were not
###  protected before, but now are. Non-symbols are silently ignored. Symbols which were already protected are
###  silently ignored.

@sjdoc Unprotect """
    Unprotect(z1,z2,...)

removes the `Protected` attribute from the symbols `z1, z2, ...`.
"""

function apprules(mx::Mxpr{:Unprotect})
    nargs = newargs()
    for i in 1:length(mx)
        ret = do_unprotect(mx,mx[i])
        if ret != false
            push!(nargs, ret)
        end
    end
    mxprcf(:List, nargs)
end

function do_unprotect(mx,a::SJSym)
    if isProtected(a)
        unprotect(a)
        return string(a)
    end
    return false
end

do_unprotect(mx,a) = false

### Protect

@sjdoc Protect """
    Protect(z1,z2,...)

add the `Protected` attribute to the lists of attributes for the symbols `z1, z2, ...`.
"""

# FIXME: Following broke after updating to v0.6.0-rc2 from an older v0.6.0
#@mkapprule Protect :nargs => 1:Inf  :nodefault => true
@mkapprule Protect :nodefault => true

@doap function Protect(args...)
    nargs = newargs()
    for x in args
        ret = do_protect(x)
        if ret != false
            push!(nargs, ret)
        end
    end
    mxprcf(:List, nargs)
end

function do_protect(a::SJSym)
    isProtected(a) && return false
    protect(a)
    string(a)
end

### Set and SetDelayed

@sjdoc Set """
    Set(a,b), a = b

Sets the value of `a` to `b`. `b` is evaluated only once, when `a=b` is evaluated.

    obj[i,j,...] = val

set a part of `obj` to `val`. `obj` can be a Symata expression
or a Julia object, such as an `Array` or `Dict`.
"""

@sjseealso_group(Set, SetDelayed, UpSet, DownValues, UpValues, OwnValues)
@sjexamp( Set,
         ("Clear(a,b,c)",""),
         ("b = a", "a"),
         ("a = 1", "1"),
         ("c = a", "1"),
         ("a = 2", "2"),
         ("b", "2"),
         ("c", "1"))

@sjdoc SetDelayed """
    SetDelayed(a,b), a := b

Whenever `a` is evaluated, `b` is evaluated and the result is assigned to `a`.

Note that `a` is not set to the value of `b` at the time `a := b` is evaluated, but
rather to the current value of `b` every time `a` is evaluated.
"""

# Set SJSym value.
# Set has HoldFirst, SetDelayed has HoldAll.

@mkapprule Set nargs => 1:Inf

# This is what Mma does.
# When is this used ?
function do_Set(mx::Mxpr{:Set}, lhs::SJSym)
    checkprotect(lhs)
    rhs = mxprcf(:Sequence)
    setsymval(lhs, rhs)
    setdefinition(lhs, mx) # This is correct. Definition includes `Set`.
    rhs
end

## Destructured. *Can* do [x,y] = [y,x]
@doap function Set(lhs::Mxpr{:List},rhs::Mxpr{:List})
    length(lhs) != length(rhs) && return mx
    for i in 1:length(lhs)
        doeval(mxpr(:Set,lhs[i],rhs[i]))
    end
    rhs
end

@mkapprule SetDelayed nargs => 1:Inf  nodefault => true

@doap SetDelayed(lhs,rhs) = setdelayed(mx,lhs,rhs)

@doap SetDelayed(lhs, rhs...) = mxpr(:SetDelayed, lhs, mxpr(:Sequence,rhs...))

@doap SetDelayed() = mx

# getsym(symname(lhs)) is because a copy of symbol is being made somewhere
# so we look up the original in the table
# SetDelayed is not correct. rhs is also evaluated because
# lhs is evaluated twice in order to find fixed point.
# we have treat this specially somehow, not ordinary evaluation.
#
# Note added (Apr 2016): This is working AFAICT
# For evaluation, we can use setsymval for both set and setdelayed
# The only difference is whether we return the rhs.
# But, to save the definitions to a file, we need to record whether we had set or setdelayed,
# so we use setdelayedval
function setdelayed(mx::Mxpr{:SetDelayed},lhs::SJSymbol, rhs)
    checkprotect(lhs)
    setsymval(lhs,rhs)   # Using this works just as well. But, for printing, we don't whether Set or SetDelayed
    setdefinition(lhs, mx)
    add_completion_symbols(string(lhs))
    Null
end

#function do_Set(mx::Mxpr{:Set}, lhs::SJSymbol, rhs)
@doap function Set(lhs::SJSymbol, rhs)
    checkprotect(lhs)
    setsymval(lhs,rhs)
    setdefinition(lhs, mx)
    add_completion_symbols(string(lhs))
    rhs
end

# Create DownValue. "function" definition
# eg f(x_) := x  defines a DownValue for the SJSym f
function setdelayed(mx,lhs::Mxpr, rhs)
    checkprotect(lhs)
    rule = mxpr(:RuleDelayed,mxpr(:HoldPattern,lhs),rhs)
    h = mhead(lhs)
    set_downvalue(mx,h,rule) # push DownValue
    if isa(h,Symbol) add_completion_symbols(string(h)) end
    rule
    Null
end

function do_Set(mx::Mxpr{:Set},lhs::Mxpr{:Part}, rhs::Mxpr{:Module})
    symerror("$mx is not implemented")
end

# Mma is not clear but seems to evaluate the first arg to the lhs (the expression
# whose part we want) exactly once. We should document what we do.
# We check is_Number several times, because we may have a Dict.
@doap function Set(lhs::Mxpr{:Part}, rhs)
    ex0 = meval(expr(lhs))  # evaluate once, eg, to get expr from symbol.
    tinds = indices(lhs)
    if isa(tinds[1],Mxpr{:Sequence}) ## is this the correct place to unwrap the Sequence ?
        tinds[1] = margs(tinds[1])
    end
    inds = map(doeval,tinds)
    val = doeval(rhs)
    ex = ex0
    exlast = ex0
    for j in 1:length(inds)-1
        exlast = ex
        ex = get_part_one_ind(ex,inds[j])
    end
    ind = inds[end] # doeval(tinds[end])
    ind = posnegi(ex,ind)
    if isa(ind,Number) && ind == 0
        if length(inds) == 1
            newexpr = mxprnewhead(ex,val)
            setsymval(lhs[1], newexpr)
            ex0 = newexpr
        else
            exlast[tinds[length(inds)-1]] = mxprnewhead(ex,val)
        end
    else
        ex[ind] = val
    end
    unsetfixed(ex0)
    val
end

# Optimize a bit. Localize variables once, not every time pattern is evaluated
setdelayed(mx,lhs::Mxpr, rhs::Mxpr{:Module}) = setdelayed(mx,lhs,localize_module!(rhs))

@doap Set(lhs::Mxpr, rhs::Mxpr{:Module}) = do_Set(mx,lhs,localize_module!(rhs))

### Increment

@sjdoc Increment """
    Increment(n)

increments the value of `n` by `1` and returns the old value.
"""

@mkapprule Increment :nargs => 1

@doap function Increment(x::SJSym)
    @checkunbound(mx,x,xval)
    do_increment1(mx,x,xval)
end

function do_increment1(mx,x,xval::T) where T<:Number
    setsymval(x,mplus(xval,1))  # maybe + is ok here.
    return xval
end

function do_increment1(mx,x,val)
    setsymval(x,doeval(mxpr(:Plus,val,1)))
    return val
end

### Decrement

@sjdoc Decrement """
    Decrement(n)

decrements the value of `n` by `1` and returns the old value.
"""

@mkapprule Decrement :nargs => 1

#function do_Decrement(mx, x::SJSym)
@doap function Decrement(x::SJSym)
    @checkunbound(mx,x,xval)
    do_decrement1(mx,x,xval)
end

function do_decrement1(mx,x,xval::Number)
    setsymval(x,mplus(xval,-1))  # maybe + is ok here.
    return xval
end

function do_decrement1(mx,x,val)
    setsymval(x,doeval(mxpr(:Plus,val,-1)))
    return val
end

### TimesBy

@sjdoc TimesBy """
    TimesBy(a,b), or a *= b

set `a` to `a * b` and returns the new value. This is currently
faster than `a = a * b` for numbers.
"""

@mkapprule TimesBy :nargs => 2

function do_TimesBy(mx::Mxpr{:TimesBy}, x::SJSym,val)
    @checkunbound(mx,x,xval)
    do_TimesBy1(mx,x,xval,val)
end

function do_TimesBy1(mx,x,xval::T, val::V) where {T<:Number,V<:Number}
    r = mmul(xval,val)
    setsymval(x,r)
    return r
end
function do_TimesBy1(mx,x,xval,val)
    setsymval(x,doeval(mxpr(:Set,x, mxpr(:Times,xval,val))))
    return symval(x)
end

#### AddTo

@sjdoc AddTo """
    AddTo(a,b), or a += b,

sets `a` to `a + b` and returns the new value. This is currently
faster than `a = a + b` for numbers.
"""

@mkapprule AddTo :nargs => 2

function do_AddTo(mx::Mxpr{:AddTo},x::SJSym,val)
    @checkunbound(mx,x,xval)
    do_AddTo1(mx,x,xval,val)
end

function do_AddTo1(mx,x,xval::T, val::V) where {T<:Number,V<:Number}
    r = mplus(xval,val)
    setsymval(x,r)
    return r
end
function do_AddTo1(mx,x,xval,val)
    setsymval(x,doeval(mxpr(:Set,x, mxpr(:Plus,xval,val))))
    return symval(x)
end

#### Dump and DumpHold

@sjdoc Dump """
    Dump(expr)

print an internal representation of `expr`.

`Dump` is similar to Julia `dump`.
"""

@sjdoc DumpHold """
    DumpHold(expr)

print an internal representation of `expr`.

`DumpHold` is similar to Julia `dump`. In contrast to `Dump`, `expr` is not evaluated before its internal
representation is printed.
"""

@sjseealso_group(Dump,DumpHold)

# DumpHold does not evaluate args before dumping
apprules(mx::T) where {T<:Union{Mxpr{:Dump},Mxpr{:DumpHold}}} = for a in margs(mx) is_SJSym(a) ? dump(getssym(a)) : dump(a) end

## TODO: completely redesign this
### Contexts

@mkapprule Contexts  :nargs => 0:1

@sjdoc Contexts """
    Contexts()

return a `List` of all contexts.
"""

@doap Contexts() = tolist(getcontexts())
#     mxpr(:List, getcontexts()...)
# end

@mkapprule ContextSymbols  :nargs => 1

@sjdoc ContextsSymbols """
    Contexts(context)

return a `List` of all symbols in `context`.
"""

@doap ContextSymbols(s) = tolist(getsymbolsincontext(s))
#     mxpr(:List, getsymbolsincontext(s)...)
# end

### UpSet

@sjdoc UpSet """
    UpSet(a(g(x_)),b),  a(g(x_)) ^= b

associate the transformation rule with `g`.
"""

@mkapprule UpSet nargs => 1:Inf
@doap UpSet(lhs,rhs) = upset(mx,lhs,rhs)

## Note, mx is needed here because the entire expression is recorded for the definition of lhs
function upset(mx,lhs::Mxpr, rhs)
    _upset(mx,lhs,rhs)
    return rhs
end

function _upset_one(mx,m::SJSym,rule)
    warncheckprotect(m)
    set_upvalue(mx, m,rule)
end

_upset_one(mx,m::Mxpr,rule) = _upset_one(mx,mhead(m),rule)

_upset_one(args...) = nothing

function _upset(mx,lhs::Mxpr, rhs)
    rule = mxpr(:RuleDelayed,mxpr(:HoldPattern,lhs),rhs)
    foreach( m -> _upset_one(mx,m,rule), lhs)
end

@mkapprule UpSetDelayed nargs => 1:Inf
@doap UpSetDelayed(lhs,rhs) = upsetdelayed(mx,lhs,rhs)

## I think the only difference with UpSet again, is we don't return the rhs
function upsetdelayed(mx,lhs::Mxpr, rhs)
    _upset(mx,lhs,rhs)
    return Null
end

# unsetfixed allows f(3) to return 9, rather than 3^2 with f(x_) = x^2
# Not sure if this is the correct solution
function do_Set(mx::Mxpr{:Set}, lhs::Mxpr, rhs)
    checkprotect(lhs)
    unsetfixed(rhs)
    # Next line. We don't want to evaluate lhs, it might have a value
    # Say g(4) == 3
    # But if  f = g, n = 4, f(n) = 1 assigns 1 to g(4)
    newlhs = mxpra(doeval(mhead(lhs)), mapmargs(doeval, margs(lhs)))
    rule = mxpr(:RuleDelayed, mxpr(:HoldPattern, newlhs), rhs)
    set_downvalue(mx, mhead(newlhs), rule) # push DownValue
    return rhs # Set always returns the rhs. This is checked
end

### Symbol

@sjdoc Symbol """
    Symbol(str)

converts the string `str` to a symbol. For example if `a` is `1`,
then Symbol("a") returns `1`.
"""

@mkapprule Symbol :nargs => 1

# function apprules(mx::Mxpr{:Symbol})
#     dosymbol(mx,mx[1])
# end

@doap Symbol(s::String) = getsym(Symbol(s))
@doap Symbol(x) = (@symwarn("Symbol: expected a string"); mx)

#### Clear

@sjdoc Clear """
    Clear(x,y,z,...)

remove the values associated with `x,y,z,...`. It does not remove
their `DownValues`.

    Clear(Out)

delete all of the saved Output lines. It actually replaces them with
the value `Null`.
"""

@sjseealso_group(Clear, ClearAll)

# 'Clear' a value. ie. set symbol's value to its name
## NOTE: Mma 1) allows 0 or more args with no message.
##           2) complains about non-symbols, but does not abort, ie. will clear remaining symbols
##           3) the symbol can be given as a string
##           4) always returns Null

@mkapprule Clear  nodefault => true

@doap Clear(args...) = (foreach(_clear, args); Null)

_clear(x) = Null

function _clear(ins::SymString)
    s = Symbol(ins)
    if s == :Out
        clear_all_output()
        return
    end
    checkprotect(s)
    clear_downvalues(s)
    clear_upvalues(s)
    clear_ownvalue_definition(s)
    setsymval(s,s)        # May 2016. This works ?
end

### ClearAll

@sjdoc ClearAll """
    ClearAll(x,y,z,...)

remove all values and `DownValues` associated with `x,y,z`.

The symbols are removed from the symbol table and will not appear in the list returned
by `UserSyms()`.

`Apply(ClearAll, UserSyms())` clears all user symbols.
"""

# FIXME. remove SymPy properties from symbol
# Remove all values associate with SJSym. values and DownValues
function apprules(mx::Mxpr{:ClearAll})  # already threaded
    for a in margs(mx)
        b = Symbol(a)
        checkprotect(b)
        delete_sym(b)
        remove_completion_symbols(b)
    end
end

### Unique

@sjdoc Unique """
    Unique()

creates a unique symbol.
"""

@mkapprule Unique :nargs => 0:1

@doap Unique() = (s = gensym(); setsymval(s,s); s)

### SymbolName

@sjdoc SymbolName """
    SymbolName(symbol)

returns the name of `symbol` as a string.
"""
@mkapprule SymbolName :nargs => 1
@doap SymbolName(x::SJSym) = string(x)

### OwnValues

@sjdoc OwnValues """
    OwnValues(s)

Returns a list of rules corresponding values bound to `s`.
"""
@mkapprule OwnValues :nargs => 1
@doap function OwnValues(x::SJSym)
    value = getsymata(x)
    value == x && return List()
    return List(mxprcf(:RuleDelayed, mxprcf(:HoldPattern, x), value))
end

### DownValues

@sjdoc DownValues """
    DownValues(s)

return a `List` of `DownValues` associated with symbol `s`. These are values
that are typically set with the declarative "function definition".

For example `f(x_) := x` sets a `DownValue` for `f`.
"""

@sjexamp( DownValues,
         ("ClearAll(f)",""),
         ("f(x_) := x^2",""),
         ("DownValues(f)", "[HoldPattern(f(x_))->(x^2)]"))

## TODO. implement options
@mkapprule DownValues :nargs => 1
@doap DownValues(s::SJSymbol) = mxpr(:List,downvalues(s)...)

### UpValues

@sjdoc UpValues """
    UpValues(s)

returns a List of UpValues associated with symbol `s`. These are values
that are typically set with `UpSet`.
"""
@mkapprule UpValues
@doap UpValues(s::SJSymbol) = mxpr(:List,upvalues(s)...)

#### HAge, FixedQ and UnFix

@sjdoc HAge """
    HAge(s)

return the timestamp for the expression or symbol `s`.
Using this timestamp to avoid unnecessary evaluation is a partially
implemented feature.
"""

@sjseealso_group(HAge, Age, FixedQ, Syms, DirtyQ, Unfix)
# Get the last-altered timestamp of an expression or symbol
apprules(mx::Mxpr{:HAge}) = hdo_getage(mx,mx[1])
hdo_getage(mx,s::SJSym) = Int(symage(s))
hdo_getage(mx,s::Mxpr) = Int(getage(s))
hdo_getage(mx,x) = mx

apprules(mx::Mxpr{:Age}) = do_getage(mx,mx[1])
do_getage(mx,s::SJSym) = Int(symage(s))
do_getage(mx,s::Mxpr) = Int(getage(s))
do_getage(mx,x) = mx

@sjdoc FixedQ """
    FixedQ(expr)

returns the status of the fixed point bit, which tells whether `expr`
is expected to evaluate to itself in the current environment. This is partially
implemented.
"""

# Get fixed-point bit. Idea is to set it if expr evaluates to itself.
apprules(mx::Mxpr{:FixedQ}) = is_fixed(mx[1])

@sjdoc Unfix """
    Unfix(expr)

unsets the fixed flag on `expr`, causing it to be evaluated.
This is a workaround for bugs that cause an expression to be marked fixed
before it is completely evaluated.
"""
@mkapprule Unfix  nodefault => true  options => Dict( :Deep => false )

function do_Unfix(mx,expr::Mxpr; Deep=false)
#    (deep,val) = kws[1]
    Deep ? deepunsetfixed(expr) : unsetfixed(expr)
    unsetfixed(expr)
    expr
end

#do_Unfix(mx,args...) = mx
do_Unfix(mx,args...;kws...) = mx

# function apprules(mx::Mxpr{:Unfix})
#     unsetfixed(mx[1])
#     mx[1]
# end

### Syms

@sjdoc Syms """
    Syms(m)

return a `List` of the symbols that the expression `m` "depends" on. The
list is wrapped in `HoldForm` in order to prevent evaluation of the symbols.
"""

function apprules(mx::Mxpr{:Syms})
    mxpr(:HoldForm, do_syms(mx[1]))
end

do_syms(mx::Mxpr) = mxpr(:List, listsyms(mx)...)
do_syms(s) = mxpr(:List,)

### BuiltIns

@sjdoc BuiltIns """
    BuiltIns()

return a `List` of all "builtin" symbols. These are in fact all symbols that
have the `Protected` attribute.
"""

## This does not cause infinite evaluation loop.
#apprules(mx::Mxpr{:BuiltIns}) = protectedsymbols()

function apprules(mx::Mxpr{:BuiltIns})
    syms = get_builtins()
    nargs = newargs(length(syms))
    copyto!(nargs, syms)
    setfixed(MListA(nargs))
end

function get_builtins()
    syms = get_system_symbols();
    symidx = something(findfirst(isequal(:ans), syms), 0)
    deleteat!(syms, symidx)  ## These two lines avoid infinite evaluation loop
    deleteat!(syms,findall((x -> match(r"^O+$",string(x)) !== nothing),syms))
#    deleteat!(syms,find((x -> match(r"^O+$",string(x)) !== nothing),syms))
    return syms
end

@sjdoc UserSyms """
    UserSyms()

return a `List` of symbols that have not been imported from the `System` namespace.
This is all user defined symbols (unless you have imported symbols from elsewhere).
"""

@mkapprule UserSyms  :nargs => 0

@doap UserSyms() = tolistfixed(usersymbols())

@sjdoc CurrentContext """
    CurrentContext()

return the name of the current context.
"""

@mkapprule CurrentContext :nargs => 0

# This does not return a context or module type, because we need to
# keep types out of the language as much as possible. Everything is
# an expression! In fact, we should probably return a string.

@doap CurrentContext() = string(get_current_context_name())


text:
function check_file_existence(url::String)

    if (isfile(url) == false)
        throw(ArgumentError("Oppps! opening file $(url): No such file or directory"))
    end
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
using Multigraphs, Graphs, SparseArrays
using Test
try
    m2 = spzeros(Int, 2, 3)
    dg = DiMultigraph(m2)
catch e
    @test e !== nothing
end
try
    m2 = spzeros(Int, 2, 2)
    m2[1, 2] = 2
    dg = DiMultigraph(m2)
catch e
    @test e === nothing
end
try
    m2 = spzeros(Int, 2, 2)
    m2[1, 2] = -1
    m2[2, 1] = -1
    dg = DiMultigraph(m2)
catch e
    @test e !== nothing
end

m = spzeros(Int, 4, 4)
m
m[1,2] = 2
m[2,1] = 1
m[2,3] = 2
m[3,2] = 1
m[3,4] += 1
m[3,4] = 0
m[4,3] += 1
m[4,3] = 0
g = DiMultigraph(m)
g = DiMultigraph(Matrix(m))
g0 = DiMultigraph(2)
@test !add_edge!(g0, 2, 3) && !rem_edge!(g0, 1, 2)
g1 = DiMultigraph(path_digraph(3))
@test adjacency_matrix(g) == m

@test is_directed(g)
@test edgetype(g) == MultipleEdge{Int, Int}
@test size(adjacency_matrix(g), 1) == 4

@test nv(g) == 4 && ne(g, count_mul = true) == 6 && ne(g) == 4

add_edge!(g1, 1, 1, 2)
@test ne(g1, count_mul = true) == 4 && ne(g1) == 3

add_vertices!(g, 3)
@test nv(g) == 7

@test has_edge(g, 1, 2, 2)
@test rem_vertices!(g, [7, 5, 4, 6])
add_edge!(g, [2, 3, 2])
rem_edge!(g, [2, 3, 2])
add_edge!(g, 2, 3)
rem_edge!(g, 2, 3)
add_edge!(g, 2, 3, 2)
rem_edge!(g, 2, 3, 1)

@test has_edge(g, 2, 3) && has_edge(g, [2, 3])
@test has_edge(g, 2, 3, 2) && has_edge(g, (2, 3, 2))
@test !has_edge(g, 2, 2) && !has_edge(g, 2, 5)
@test has_vertex(g, 1) && !has_vertex(g, 5)

collect(edges(g))
@test inneighbors(g, 2) == outneighbors(g, 2)
@test degree(g, 2) != indegree(g, 2)
add_vertex!(g)
@test indegree(g) != outdegree(g)

dmg0 = DiMultigraph(0)
@test nv(dmg0) == ne(dmg0) == 0


text:
# -------------------------------------------------------
#           Green's function indexing stuff
# -------------------------------------------------------
"""
Get linear site index of cartesian coordinates (x,y) respecting PBC.
"""
@inline siteidx(mc,sql,x,y) = siteidx(mc.p.L, mc.l.sites, sql, x, y)
@inline siteidx(L,N,sql,x,y) = begin
  xpbc = mod1(x, L)
  ypbc = mod1(y, L)
  sql[ypbc, xpbc] # to linear idx
end

"""
Get column/row idx of particular flv ∈ (xu, yd, xd, yu) and site ∈ 1:N in Green's function.
"""
@inline greensidx(N::Int, flv, site) = (flv-1)*N + site
@inline greensidx(mc::AbstractDQMC, flv, site) = greensidx(mc.l.sites, flv, site)

"""
Access full (4*N, 4*N) Green's function for any OPDIM efficiently.
"""
@inline function G(mc, i, j, greens=mc.s.greens)
  gt = eltype(greens)
  N = mc.l.sites
  opdim = mc.p.opdim
  half = 2*N

  @inbounds if opdim == 3
    return greens[i,j]

  else
    # Virtually expand Green's function
    if i>half && j<=half # lower left block
      return zero(gt)

    elseif i<=half && j>half # upper right block
      return zero(gt)

    elseif i>half && j>half # lower right block
      return conj(greens[i - half, j - half])

    else # upper left block
      return greens[i, j]
    end
  end
end

"""
Access full (4*N, 4*N) `Gtilde = I - G` where G is the Green's function for any OPDIM efficiently.
"""
@inline function Gtilde(mc, i, j, greens=mc.s.greens)
  if i == j # on diagonal
    return 1 - G(mc, i, j, greens)
  else
    return -G(mc, i, j, greens)
  end
end

"""
Construct full (4*N, 4*N) Green's function for any OPDIM efficiently.
"""
function fullG(mc, greens=mc.s.greens)
  gt = eltype(greens)
  N = mc.l.sites
  g = zeros(gt, 4*N, 4*N)

  @inbounds for j in 1:4*N
    for i in 1:4*N
      g[i,j] = G(mc, i, j, greens)
    end
  end
  return g
end

"""
Construct full (4*N, 4*N) `Gtilde = I - G` where G is the Green's function for any OPDIM efficiently.
"""
function fullGtilde(mc, greens=mc.s.greens)
  gt = eltype(greens)
  N = mc.l.sites
  g = zeros(gt, 4*N, 4*N)

  @inbounds for j in 1:4*N
    for i in 1:4*N
      g[i,j] = Gtilde(mc, i, j, greens)
    end
  end
  return g
end

# test: fullGtilde(mc) == (I - fullG(mc))

# -------------------------------------------------------
#         Postprocessing/Analysis
# -------------------------------------------------------
# Go from xup, ydown, xdown, yup -> xup, yup, xdown, ydown
function permute_greens(greens::AbstractMatrix)
  perm = [1,4,3,2] # flv*spin: xup, ydown, xdown, yup -> xup, yup, xdown, ydown
  N = Int(sqrt(length(greens))/4)
  return reshape(reshape(greens, (N,4,N,4))[:,perm,:,perm], (4*N,4*N)); # rfs1, rfs2
end

"""
    fft_greens(greens; fftshift=false, flv=4, L=sqrt(size(greens, 1)/flv))

Fourier transform real-space Green's function (`flv*N, flv*N`) to momentum space.

If `fftshift == true` momenta are shifted such that k=0 (Γ point) is centered.
"""
function fft_greens(greens::AbstractMatrix; fftshift=false, flv::Integer=4, L::Integer=Int(sqrt(size(greens, 1)/flv)))
  N = L^2
  g = reshape(greens, (L,L,flv,L,L,flv)) # y1,x1,f1, y2,x2,f2
  gk = ifft( fft(g, (1,2)), (4,5))
  if fftshift
    gk = FFTW.fftshift(gk, (1,2,4,5))
  end
  gk = reshape(gk, (flv*N,flv*N))
end

fft_greens(mc::AbstractDQMC, greens::AbstractMatrix; fftshift=false) = fft_greens(greens; fftshift=fftshift, flv=mc.p.flv, L=mc.p.L)

"""
    ifft_greens(gk; ifftshift=false, flv=4, L=sqrt(size(greens, 1)/flv))

Inverse Fourier transform momentum-space Green's function (`flv*N, flv*N`) to real space.
Assumes fftshifted momenta as it applies ifftshift before the transformation.

Inverse of `fft_greens`.
"""
function ifft_greens(gk::AbstractMatrix; ifftshift=false, flv::Integer=4, L::Integer=Int(sqrt(size(gk, 1)/flv)))
  N = L^2
  g = reshape(gk, (L,L,flv,L,L,flv)) # y1,x1,f1, y2,x2,f2
  if ifftshift
    g = FFTW.ifftshift(g, (1,2,4,5))
  end
  greens = fft( ifft(g, (1,2)), (4,5))
  greens = reshape(greens, (flv*N,flv*N))
end

ifft_greens(mc::AbstractDQMC, gk::AbstractMatrix; ifftshift=false) = ifft_greens(gk; ifftshift=ifftshift, flv=mc.p.flv, L=mc.p.L)

"""
    fftmomenta(L; fftshift=false) -> qs (= qxs = qys)

Returns momenta of a finite size system (Brillouin zone discretization).

If `fftshift==true` the momenta are shifted such that k=0 (Γ point) is centered.
In this case, the output is identical to `collect(range(-pi, pi, length=7))[1:end-1]`.
"""
function fftmomenta(L::Int; fftshift=false)
  qs = 2*pi*fftfreq(L) # 2*pi because fftfreq returns linear momenta
  if fftshift
    qs = FFTW.fftshift(qs)
  end

  # ωs = 2*pi*fftfreq(M, 1/Δτ)
  return qs
end

fftmomenta(mc::AbstractDQMC; kwargs...) = fftmomenta(mc.p.L; kwargs...)

"""
Checks whether a correlation function, like etpc_minus, is reflection symmetric.
"""
function is_reflection_symmetric(x; tol=1e-12)
  @views inner = x[2:end, 2:end]

  tmp = inner - reverse(inner, dims=1)
  maximum(real(tmp)) < tol || return false
  maximum(imag(tmp)) < tol || return false

  tmp = inner - reverse(inner, dims=2)
  maximum(real(tmp)) < tol || return false
  maximum(imag(tmp)) < tol || return false

  @views border_top = x[1,2:end]
  tmp = border_top - reverse(border_top)
  maximum(real(tmp)) < tol || return false
  maximum(imag(tmp)) < tol || return false

  @views border_left = x[2:end, 1]
  tmp = border_left - reverse(border_left)
  maximum(real(tmp)) < tol || return false
  maximum(imag(tmp)) < tol || return false

  return true
end

# -------------------------------------------------------
#             Equal times Green's function
# -------------------------------------------------------
function measure_greens_and_logdet(mc::AbstractDQMC, safe_mult::Int=mc.p.safe_mult, slice::Int=1)
  greens, greens_logdet = calc_greens_and_logdet(mc, slice, safe_mult)
  effective_greens2greens!(mc, greens)
  return greens, greens_logdet
end

function measure_greens(mc::AbstractDQMC, safe_mult::Int=mc.p.safe_mult, slice::Int=1)
  greens = calc_greens(mc, slice, safe_mult)
  effective_greens2greens!(mc, greens)
  return greens
end

# -------------------------------------------------------
#   Equal times Green's function (all of them, i.e. ∀τ)
# -------------------------------------------------------
"""
Calculates (true i.e. not effective) equal-time Green's functions
at every time slice, i.e. G(tau) for all tau in 1:M where M is # slices.

This method is using the dqmc stack and its logic.
"""
function measure_all_greens(mc::AbstractDQMC)
  etgfs = calc_all_greens(mc)
  effective_greens2greens!.(Ref(mc), etgfs)
  return etgfs
end

"""
Calculates (effective) equal-time Green's functions at every
time slice, i.e. G(tau) for all tau in 1:M where M is # slices.

This method is using the dqmc stack and its logic.
"""
function calc_all_greens(mc::AbstractDQMC)
  build_stack(mc)
  M = mc.p.slices
  safe_mult = mc.p.safe_mult

  etgfs = Vector{Matrix{geltype(mc)}}(undef, M)

  for tau in M:-1:1
    propagate(mc)
    # @assert tau == mc.s.current_slice
    @inbounds etgfs[tau] = copy(mc.s.greens)
  end

  return etgfs
end

"""
Calculates (effective) equal-time Green's functions at every
time slice, i.e. G(tau) for all tau in 1:M where M is # slices.
"""
function calc_all_greens_explicit(mc::AbstractDQMC)
  # TODO: remove calculation redundancy by storing intermediate Bchains
  M = mc.p.slices
  safe_mult = mc.p.safe_mult

  etgfs = Vector{Matrix{geltype(mc)}}(undef, M)

  etgfs[1] = calc_greens(mc, 1)

  for tau in 2:M
    if mod(tau-1, safe_mult) == 0
      etgfs[tau] = calc_greens(mc, tau)
    else
      etgfs[tau] = wrap_greens(mc, etgfs[tau-1], tau-1, 1)
    end
  end

  return etgfs
end

# -------------------------------------------------------
#         Occupation
# -------------------------------------------------------
"""
Overall occupation, averaged over sites and fermion flavors (both spin and band).

Example: half-filling of spin 1/2 fermions on a lattice corresponds to `n=0.5`, not `n=1`.
"""
function occupation(mc::AbstractDQMC, greens::AbstractMatrix=mc.s.greens)
  n = mean(1 .- diag(greens))
  # only per site:
  # N = mc.l.sites
  # M = mc.p.slices
  # n = 1/(N*M) * sum(1 .- diag(greens))
  real(n)
end

"""
Flavor resolved occupations (only averaged over sites).
"""
function occupations_flv(mc::AbstractDQMC, greens::AbstractMatrix=mc.s.greens)
  N = mc.l.sites
  ns = mean.(Iterators.partition(1 .- diag(greens), N))
  # only per site:
  # N = mc.l.sites
  # M = mc.p.slices
  # n = 1/(N*M) * sum(1 .- diag(greens))
  real(ns)
end

# -------------------------------------------------------
#           Equal-time Pairing Correlations
# -------------------------------------------------------
function allocate_etpc!(mc)
  L = mc.p.L
  meas = mc.s.meas

  meas.etpc_minus = zeros(Float64,L,L)
  meas.etpc_plus = zeros(Float64,L,L)

  println("Allocated memory for ETPC measurements.")
  nothing
end

"""
Calculate equal time pairing susceptibilities (Eq. 4 in notes)

  Pη(y,x) = << Δη^†(r)Δη(0) >> = etpc_plus/minus(y,x)

where η=± (s and d-wave), r=(y,x), and <<·>> indicates a fermion average for fixed ϕ.

x,y ∈ [0,L-1]

Details:
- we mean over r_0, that is "0".
- we do not mean over τ_0=0 (τ not shown above).
"""
function measure_etpc!(mc::AbstractDQMC, greens::AbstractMatrix)
  L = mc.p.L
  N = mc.l.sites
  Pm = mc.s.meas.etpc_minus # d-wave
  Pp = mc.s.meas.etpc_plus # s-wave

  fill!(Pm, 0.0)
  fill!(Pp, 0.0)

  # we only need four combinations of j1,j2,j3,j4 for etpc_ev
  xu, yd, xd, yu = 1,2,3,4
  js = ((xd,xu,xu,xd), (xd,xu,yu,yd), (yd,yu,xu,xd), (yd,yu,yu,yd))

  sql = reshape(1:N,L,L)

  @inbounds for x in 0:(L-1), y in 0:(L-1) # r (displacement)
    for x0 in 1:L, y0 in 1:L # r0 (origin)
      r1 = r2 = siteidx(mc, sql, x0+x, y0+y)
      r3 = r4 = siteidx(mc, sql, x0, y0)

      ev1 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 1)
      ev2 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 2)
      ev3 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 3)
      ev4 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 4)

      Pm[y+1,x+1] += ev1 - ev2 - ev3 + ev4
      Pp[y+1,x+1] += ev1 + ev2 + ev3 + ev4
    end
  end

  # r0 mean
  Pm ./= N
  Pp ./= N

  nothing
end

"""
Pairing correlation expectation value

  etpc_ev = << c_j1(r1) c_j2(r2) c_j3(r3)^† c_j4(r4)^† >>

(Wick's theorem.)
"""
@inline function _pc_ev(mc, greens, r1, r2, r3, r4, js, j)
  N = mc.l.sites

  i1 = greensidx(N, js[j][1], r1)
  i2 = greensidx(N, js[j][2], r2)
  i3 = greensidx(N, js[j][3], r3)
  i4 = greensidx(N, js[j][4], r4)

  # Wick's theorem. Result is real because greens is symm. in r and has anti-unitary symm.
  result = G(mc, i1, i4, greens)*G(mc, i2, i3, greens) -
              G(mc, i1, i3, greens)*G(mc, i2, i4, greens)
  # @show maximum(imag(result))
  # @assert isreal(result)
  return real(result)
end

# """
# FFT of ETPC, i.e. P(ky, kx)
# """
# measure_etpc_k!(args...) = begin measure_etpc!(args...); rfft(mc.s.meas.etpc, (1,2)); end

# """
# P = Σ_r P(r), with P(r) from `etpc`.
# """
# measure_etpc_uniform!(mc, η, greens) = begin measure_etpc!(mc, η, greens); sum(mc.s.meas.etpc); end

# """
# P = P(q=0), with P(q)=fft(P(r)) from `etpc_k`. Slower than `etpc_uniform!`.
# """
# measure_etpc_uniform_alternative!(mc, η, greens) = real(measure_etpc_k!(mc, η, greens)[1,1])

# -------------------------------------------------------
#      Zero-Frequency (static) Pairing Correlations
# -------------------------------------------------------
function allocate_zfpc!(mc)
  L = mc.p.L
  meas = mc.s.meas

  meas.zfpc_minus = zeros(Float64,L,L)
  meas.zfpc_plus = zeros(Float64,L,L)

  println("Allocated memory for ZFPC measurements.")
  nothing
end

"""
Calculate zero-frequency pairing susceptibilities (Eq. 3 in notes)

  Pη(y,x) = int_tau << Δη^†(r, tau)Δη(0, 0) >> = zfpc_plus/minus(y,x)

where η=± (s and d-wave), r=(y,x), and <<·>> indicates a fermion average for fixed ϕ.

x,y ∈ [0,L-1]

Details:
- we mean over r_0, that is "0".
- we do not mean over τ_0=0.
"""
function measure_zfpc!(mc::AbstractDQMC, Gt0s::AbstractVector{T}) where T <: AbstractMatrix
  L = mc.p.L
  N = mc.l.sites
  M = mc.p.slices
  Pm = mc.s.meas.zfpc_minus # d-wave
  Pp = mc.s.meas.zfpc_plus # s-wave

  fill!(Pm, 0.0)
  fill!(Pp, 0.0)

  # we only need four combinations of j1,j2,j3,j4 for etpc_ev
  xu, yd, xd, yu = 1,2,3,4
  js = ((xd,xu,xu,xd), (xd,xu,yu,yd), (yd,yu,xu,xd), (yd,yu,yu,yd))

  sql = reshape(1:N,L,L)

  @inbounds for tau in 1:M
    greens = Gt0s[tau]

    # Identical to etpc now
    @inbounds for x in 0:(L-1), y in 0:(L-1) # r (displacement)
      for x0 in 1:L, y0 in 1:L # r0 (origin)
        r1 = r2 = siteidx(mc, sql, x0+x, y0+y)
        r3 = r4 = siteidx(mc, sql, x0, y0)

        ev1 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 1)
        ev2 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 2)
        ev3 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 3)
        ev4 = _pc_ev(mc, greens, r1, r2, r3, r4, js, 4)

        Pm[y+1,x+1] += ev1 - ev2 - ev3 + ev4
        Pp[y+1,x+1] += ev1 + ev2 + ev3 + ev4
      end
    end
  end

  # r0 mean
  Pm ./= N
  Pp ./= N

  nothing
end

# -------------------------------------------------------
#         Equal-time Charge Density Correlations
# -------------------------------------------------------
function allocate_etcdc!(mc)
  L = mc.p.L
  meas = mc.s.meas

  meas.etcdc_minus = zeros(Float64,L,L)
  meas.etcdc_plus = zeros(Float64,L,L)

  println("Allocated memory for ETCDC measurements.")
  nothing
end

"""
Calculate equal time charge density susceptibilities

  Cη(y,x) = << Δη^†(r)Δη(0) >> = etcdc_plus/minus(y,x)

where η=± (s and d-wave), r=(y,x), and <<·>> indicates a fermion average for fixed ϕ.

x,y ∈ [0,L-1]

Details:
- we mean over r_0, that is "0".
- we do not mean over τ_0=0 (τ not shown above).
"""
function measure_etcdc!(mc::AbstractDQMC, greens::AbstractMatrix)
  L = mc.p.L
  N = mc.l.sites
  Cm = mc.s.meas.etcdc_minus # d-wave
  Cp = mc.s.meas.etcdc_plus # s-wave

  fill!(Cm, 0.0)
  fill!(Cp, 0.0)

  # xu, yd, xd, yu = 1,2,3,4
  X, Y = 1, 2
  UP, DOWN = 0, 2
  flv = (band, spin) -> band + spin

  sql = reshape(1:N,L,L)

  @inbounds for x0 in 1:L, y0 in 1:L # r0 (origin)
    for s1 in (UP, DOWN), s2 in (UP, DOWN)
      for x in 0:(L-1), y in 0:(L-1) # r (displacement)
        ri = siteidx(mc, sql, x0+x, y0+y)
        r0 = siteidx(mc, sql, x0, y0)

        ev1 = _etcdc_ev(mc, greens, ri, ri, r0, r0, flv(X, s1), flv(X, s1), flv(X, s2), flv(X, s2))
        ev2 = _etcdc_ev(mc, greens, ri, ri, r0, r0, flv(X, s1), flv(X, s1), flv(Y, s2), flv(Y, s2))
        ev3 = _etcdc_ev(mc, greens, ri, ri, r0, r0, flv(Y, s1), flv(Y, s1), flv(X, s2), flv(X, s2))
        ev4 = _etcdc_ev(mc, greens, ri, ri, r0, r0, flv(Y, s1), flv(Y, s1), flv(Y, s2), flv(Y, s2))

        Cm[y+1,x+1] += ev1 - ev2 - ev3 + ev4
        Cp[y+1,x+1] += ev1 + ev2 + ev3 + ev4
      end
    end
  end

  # r0 mean
  Cm ./= N
  Cp ./= N

  nothing
end

"""
Charge density correlation expectation value

  etcdc_ev = << c_j1^†(r1) c_j2(r2) c_j3^†(r3) c_j4(r4) >>

(Wick's theorem.)
"""
@inline function _etcdc_ev(mc, greens, r1, r2, r3, r4, j1, j2, j3, j4)
  N = mc.l.sites
  i1 = greensidx(N, j1, r1)
  i2 = greensidx(N, j2, r2)
  i3 = greensidx(N, j3, r3)
  i4 = greensidx(N, j4, r4)

  uncorr = (kd(i1, i2) - G(mc, i2, i1, greens)) * (kd(i3, i4) - G(mc, i4, i3, greens)) # optimize wrt speed?
  corr1 = (kd(i1, i4) - G(mc, i4, i1, greens))
  corr2 = G(mc, i2, i3, greens)

  return uncorr + corr1 * corr2
end

# -------------------------------------------------------
#  Zero-Frequency (static) Charge Density Correlations
# -------------------------------------------------------
function allocate_zfcdc!(mc)
  L = mc.p.L
  meas = mc.s.meas

  meas.zfcdc_minus = zeros(Float64,L,L)
  meas.zfcdc_plus = zeros(Float64,L,L)

  println("Allocated memory for ZFCDC measurements.")
  nothing
end

"""
Calculate zero-frequency charge density susceptibilities

  Cη(y,x) = int_tau << Δη^†(r, tau) Δη(0, 0) >> = zfcdc_plus/minus(y,x)

where η=± (s and d-wave), r=(y,x), and <<·>> indicates a fermion average for fixed ϕ.

x,y ∈ [0,L-1]

The input `greens` can either be a single ETGF or all ETGFs, i.e. G(tau).

Details:
- we mean over r_0, that is "0".
- we do not mean over τ_0=0 (τ not shown above).
"""
function measure_zfcdc!(mc::AbstractDQMC, greens::Union{V, W}, Gt0s::AbstractVector{T}, G0ts::AbstractVector{T}) where {T <: AbstractMatrix, V <: AbstractVector, W <: AbstractMatrix}
  L = mc.p.L
  N = mc.l.sites
  M = mc.p.slices
  Cm = mc.s.meas.zfcdc_minus # d-wave
  Cp = mc.s.meas.zfcdc_plus # s-wave

  fill!(Cm, 0.0)
  fill!(Cp, 0.0)

  # xu, yd, xd, yu = 1,2,3,4
  X, Y = 1, 2
  UP, DOWN = 0, 2
  flv = (band, spin) -> band + spin

  sql = reshape(1:N,L,L)

  if typeof(greens) <: AbstractMatrix
    get_Gtau = (tau) -> greens
    get_G0 = (tau) -> greens
  else
    get_Gtau = (tau) -> greens[tau]
    get_G0 = (tau) -> greens[1]
  end

  @inbounds for tau in 1:M
    Gtau = get_Gtau(tau)
    G0 = get_G0(tau)
    Gt0 = Gt0s[tau]
    G0t = G0ts[tau]

    for x0 in 1:L, y0 in 1:L # r0 (we average over r0)
      for s1 in (UP, DOWN), s2 in (UP, DOWN)
        for x in 0:(L-1), y in 0:(L-1) # r (displacement)
          ri = siteidx(mc, sql, x0+x, y0+y)
          r0 = siteidx(mc, sql, x0, y0)

          ev1 = _zfcdc_ev(mc, tau, Gtau, G0, Gt0, G0t, ri, ri, r0, r0, flv(X, s1), flv(X, s1), flv(X, s2), flv(X, s2))
          ev2 = _zfcdc_ev(mc, tau, Gtau, G0, Gt0, G0t, ri, ri, r0, r0, flv(X, s1), flv(X, s1), flv(Y, s2), flv(Y, s2))
          ev3 = _zfcdc_ev(mc, tau, Gtau, G0, Gt0, G0t, ri, ri, r0, r0, flv(Y, s1), flv(Y, s1), flv(X, s2), flv(X, s2))
          ev4 = _zfcdc_ev(mc, tau, Gtau, G0, Gt0, G0t, ri, ri, r0, r0, flv(Y, s1), flv(Y, s1), flv(Y, s2), flv(Y, s2))

          Cm[y+1,x+1] += ev1 - ev2 - ev3 + ev4
          Cp[y+1,x+1] += ev1 + ev2 + ev3 + ev4
        end
      end
    end
  end

  # r0 mean
  Cm ./= N
  Cp ./= N

  nothing
end

"""
Zero-frequency charge density correlation expectation value

  zfcdc_ev = << c_j1^†(r1, tau) c_j2(r2, tau) c_j3^†(r3, 0) c_j4(r4, 0) >>

(Wick's theorem.)
"""
@inline function _zfcdc_ev(mc, tau, Gtau, G0, Gt0, G0t, r1, r2, r3, r4, j1, j2, j3, j4)
  N = mc.l.sites
  i1 = greensidx(N, j1, r1)
  i2 = greensidx(N, j2, r2)
  i3 = greensidx(N, j3, r3)
  i4 = greensidx(N, j4, r4)

  uncorr = (kd(i1, i2) - G(mc, i2, i1, Gtau)) * (kd(i3, i4) - G(mc, i4, i3, G0)) # optimize wrt speed?
  # case τ=1: corr1 = I - G(0,0) = I - (G(0, τ=0) + I) = - G(0, τ=0)
  corr1 = -G(mc, i4, i1, G0t)
  corr2 = G(mc, i2, i3, Gt0)

  return uncorr + corr1 * corr2
end

# -------------------------------------------------------
#  Zero-Frequency (static) Current-Current Correlations
# -------------------------------------------------------
function allocate_zfccc!(mc)
  L = mc.p.L
  meas = mc.s.meas

  meas.zfccc_xx = zeros(Float64,L,L)
  meas.zfccc_yy = zeros(Float64,L,L)

  println("Allocated memory for ZFCCC measurements.")
  nothing
end

"""
Calculate zero-frequency current-current susceptibilities

  Λxx(y,x) = int_tau << j_x(r, tau) j_x(0, 0) >> = zfccc(y,x)

where r=(y,x), and <<·>> indicates a fermion average for fixed ϕ.

x,y ∈ [0,L-1]

The input `greens` can either be a single ETGF or all ETGFs, i.e. G(tau).

Details:
- we mean over r_0, that is "0".
- we do not mean over τ_0=0 (τ not shown above).
"""
function measure_zfccc!(mc::AbstractDQMC, greens::Union{V, W}, Gt0s::AbstractVector{S}, G0ts::AbstractVector{S}) where {S <: AbstractMatrix, V <: AbstractVector, W <: AbstractMatrix}
  L = mc.p.L
  N = mc.l.sites
  M = mc.p.slices
  T = mc.l.hopping_matrix

  # xu, yd, xd, yu = 1,2,3,4
  X, Y = 1, 2
  UP, DOWN = 0, 2
  flv = (band, spin) -> band + spin

  sql = reshape(1:N,L,L)

  if typeof(greens) <: AbstractMatrix
    get_Gtau = (tau) -> greens
    get_G0 = (tau) -> greens
  else
    get_Gtau = (tau) -> greens[tau]
    get_G0 = (tau) -> greens[1]
  end

  for dir in (:x, :y)
    Lambda = dir == :x ? mc.s.meas.zfccc_xx : mc.s.meas.zfccc_yy
    fill!(Lambda, 0.0)

    @inbounds for tau in 1:M
      Gtau = get_Gtau(tau)
      G0 = get_G0(tau)
      Gt0 = Gt0s[tau]
      G0t = G0ts[tau]

      # TODO: use symmetries to reduce a1,a2,s1,s2 sum?

      for x0 in 1:L, y0 in 1:L # r0 (we average over r0)
        for a1 in (X, Y), a2 in (X, Y)
          for s1 in (UP, DOWN), s2 in (UP, DOWN)
            for x in 0:(L-1), y in 0:(L-1) # r (displacement)
              ri = siteidx(mc, sql, x0+x, y0+y)
              rj = dir == :x ? siteidx(mc, sql, x0+x+1, y0+y) : siteidx(mc, sql, x0+x, y0+y+1) # ri + x̂ / ri + ŷ
              r0 = siteidx(mc, sql, x0, y0)
              r0p = dir == :x ? siteidx(mc, sql, x0+1, y0) : siteidx(mc, sql, x0, y0+1) # r0 + x̂ / r0 + ŷ

              ais = greensidx(N, flv(a1, s1), ri)
              ajs = greensidx(N, flv(a1, s1), rj)
              ap0sp = greensidx(N, flv(a2, s2), r0)
              ap0psp = greensidx(N, flv(a2, s2), r0p)

              # hoppings
              tij = G(mc, ais, ajs, T)
              tji = G(mc, ajs, ais, T)
              t00p = G(mc, ap0sp, ap0psp, T)
              t0p0 = G(mc, ap0psp, ap0sp, T)

              # Uncorrelated part
              Lambda[y+1, x+1] += (tij * G(mc, ajs, ais, Gtau) - tji * G(mc, ais, ajs, Gtau)) *
                                  (t00p * G(mc, ap0psp, ap0sp, G0) - t0p0 * G(mc, ap0sp, ap0psp, G0))

              # Correlated part
              Lambda[y+1, x+1] += - tij * t00p * G(mc, ap0psp, ais, G0t) * G(mc, ajs, ap0sp, Gt0) +
                                  + tij * t0p0 * G(mc, ap0sp, ais, G0t) * G(mc, ajs, ap0psp, Gt0) +
                                  + tji * t00p * G(mc, ap0psp, ajs, G0t) * G(mc, ais, ap0sp, Gt0) +
                                  - tji * t0p0 * G(mc, ap0sp, ajs, G0t) * G(mc, ais, ap0psp, Gt0)
            end
          end
        end
      end
    end

    # r0 mean + overall minus sign due to i^2
    Lambda ./= -N
  end

  nothing
end

# -------------------------------------------------------
#  Superfluid density
# -------------------------------------------------------

"""
    measure_sfdensity(mc, zfccc_xx=mc.s.meas.zfccc_xx, zfccc_yy=mc.s.meas.zfccc_yy)

Calculate the superfluid density from zero-frequency
current-current correlations `zfccc_xx`, `zfccc_yy`.
"""
function measure_sfdensity(mc, zfccc_xx=mc.s.meas.zfccc_xx, zfccc_yy=mc.s.meas.zfccc_yy)
  Λxxq = rfft(real(zfccc_xx))
  Λyyq = rfft(real(zfccc_yy))

  # Our model is symmetric under the combination of a pi/2 rotation, a PH trafo and the exchange of the bands
  # Hence, Λxx ≈ Λyy' (transpose == rotate coordiantes by pi/2)
  ρs = 1/4 * (real(Λxxq[1,2]) - real(Λxxq[2,1]))
  # Otherwise
  # ρs = 1/8 * (real(Λxxq[1,2]) - real(Λxxq[2,1]) + real(Λyyq[1,2]) - real(Λyyq[2,1]))
end

# -------------------------------------------------------
#         Equal times Green's function (effective)
# -------------------------------------------------------
"""
QR DECOMPOSITION: Calculate effective(!) Green's function (direct, i.e. without stack)
"""
function calc_Bchain(mc::AbstractDQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  # Calculate Ul, Dl, Tl =B(stop) ... B(start)
  flv = mc.p.flv
  N = mc.l.sites
  G = geltype(mc)

  @assert 0 < start <= mc.p.slices
  @assert 0 < stop <= mc.p.slices
  @assert start <= stop

  U = Matrix{G}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  T = Matrix{G}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in start:stop
    if mod(k,safe_mult) == 0 || k == stop # always decompose in the end
      multiply_B_left!(mc,k,U)
      rmul!(U, Diagonal(D))
      U, Tnew = decompose_udt!(U, D)
      mul!(mc.s.tmp, Tnew, T)
      T .=  mc.s.tmp
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_B_left!(mc,k,U)
    end
  end
  return (U,D,T,svs)
end

# Calculate Ul, Dl, Tl = [B(stop) ... B(start)]^(-1) = B(start)^(-1) ... B(stop)^(-1)
function calc_Bchain_inv(mc::AbstractDQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  flv = mc.p.flv
  N = mc.l.sites
  G = geltype(mc)

  @assert 0 < start <= mc.p.slices
  @assert 0 < stop <= mc.p.slices
  @assert start <= stop

  U = Matrix{G}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  T = Matrix{G}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in reverse(start:stop)
    if mod(k,safe_mult) == 0 || k == start # always decompose in the end
      multiply_B_inv_left!(mc,k,U)
      rmul!(U, Diagonal(D))
      U, Tnew = decompose_udt!(U, D)
      mul!(mc.s.tmp, Tnew, T)
      T .=  mc.s.tmp
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_B_inv_left!(mc,k,U)
    end
  end
  return (U,D,T,svs)
end

# Calculate (Ur, Dr, Tr)' = B(stop) ... B(start)  => Ur,Dr, Tr = B(start)' ... B(stop)'
function calc_Bchain_dagger(mc::AbstractDQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  flv = mc.p.flv
  N = mc.l.sites
  G = geltype(mc)

  @assert 0 < start <= mc.p.slices
  @assert 0 < stop <= mc.p.slices
  @assert start <= stop

  U = Matrix{G}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  T = Matrix{G}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in reverse(start:stop)
    if mod(k,safe_mult) == 0 || k == start # always decompose in the end
      multiply_daggered_B_left!(mc,k,U)
      rmul!(U, Diagonal(D))
      U, Tnew = decompose_udt!(U, D)
      mul!(mc.s.tmp, Tnew, T)
      T .=  mc.s.tmp
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_daggered_B_left!(mc,k,U)
    end
  end
  return (U,D,T,svs)
end

# Calculate G(slice) = [1+B(slice-1)...B(1)B(M) ... B(slice)]^(-1) and its singular values in a stable manner
function calc_greens(mc::AbstractDQMC, slice::Int=mc.s.current_slice, safe_mult::Int=mc.p.safe_mult)
  s = mc.s
  _calc_greens_helper(mc, slice, safe_mult)

  rmul!(s.U, Diagonal(s.d))
  return s.U * s.T
end
function calc_greens_and_logdet(mc::AbstractDQMC, slice::Int=mc.s.current_slice, safe_mult::Int=mc.p.safe_mult)
  s = mc.s
  _calc_greens_helper(mc, slice, safe_mult)

  ldet = real(log(complex(det(s.U))) + sum(log.(s.d)) + log(complex(det(s.T))))

  rmul!(s.U, Diagonal(s.d))
    return s.U * s.T, ldet
end

# Result in s.U, s.d, and s.T
function _calc_greens_helper(mc::AbstractDQMC, slice::Int, safe_mult::Int)
  flv = mc.p.flv
  N = mc.l.sites
  s = mc.s
  G = geltype(mc)

  # Calculate Ur,Dr,Tr=B(slice)' ... B(M)'
  Ur, Dr, Tr = calc_Bchain_dagger(mc,slice,mc.p.slices, safe_mult)

  # Calculate Ul,Dl,Tl=B(slice-1) ... B(1)
  if slice-1 >= 1
    Ul, Dl, Tl = calc_Bchain(mc,1,slice-1, safe_mult)
  else
    Ul = Matrix{G}(I, flv*N, flv*N)
    Dl = ones(Float64, flv*N)
    Tl = Matrix{G}(I, flv*N, flv*N)
  end

  inv_one_plus_two_udts!(mc, s.U,s.d,s.T, Ul, Dl, Tl, Ur,Dr,Tr)

  nothing
end

# -------------------------------------------------------
#    Effective Green's function -> Green's function
# -------------------------------------------------------
function effective_greens2greens!(mc::DQMC_CBTrue, greens::AbstractMatrix)
  chkr_hop_half_minus = mc.l.chkr_hop_half
  chkr_hop_half_plus = mc.l.chkr_hop_half_inv
  n_groups = mc.l.n_groups
  tmp = mc.s.tmp

  @inbounds @views begin
      for i in reverse(1:n_groups)
        mul!(tmp, greens, chkr_hop_half_minus[i])
        greens .= tmp
      end
      for i in reverse(1:n_groups)
        mul!(tmp, chkr_hop_half_plus[i], greens)
        greens .= tmp
      end
  end
  nothing
end
function effective_greens2greens!(mc::DQMC_CBTrue, U::AbstractMatrix, T::AbstractMatrix)
  chkr_hop_half_minus = mc.l.chkr_hop_half
  chkr_hop_half_plus = mc.l.chkr_hop_half_inv
  n_groups = mc.l.n_groups
  tmp = mc.s.tmp

  @inbounds @views begin
      for i in reverse(1:n_groups)
        mul!(tmp, T, chkr_hop_half_minus[i])
        T .= tmp
      end
      for i in reverse(1:n_groups)
        mul!(tmp, chkr_hop_half_plus[i], U)
        U .= tmp
      end
  end
    nothing
end
function greens2effective_greens!(mc::DQMC_CBTrue, greens::AbstractMatrix)
  chkr_hop_half_minus = mc.l.chkr_hop_half
  chkr_hop_half_plus = mc.l.chkr_hop_half_inv
  n_groups = mc.l.n_groups

  @inbounds @views begin
      for i in 1:n_groups
        mul!(tmp, greens, chkr_hop_half_plus[i])
        greens .= tmp
      end
      for i in 1:n_groups
        mul!(tmp, chkr_hop_half_minus[i], greens)
        greens .= tmp
      end
  end
  nothing
end
function effective_greens2greens!(mc::DQMC_CBFalse, greens::AbstractMatrix)
  eTminus = mc.l.hopping_matrix_exp
  eTplus = mc.l.hopping_matrix_exp_inv
  tmp = mc.s.tmp

  mul!(tmp, greens, eTminus)
  mul!(greens, eTplus, tmp)
  nothing
end
function effective_greens2greens!(mc::DQMC_CBFalse, U::AbstractMatrix, T::AbstractMatrix)
  eTminus = mc.l.hopping_matrix_exp
  eTplus = mc.l.hopping_matrix_exp_inv

  T .= T * eTminus
  U .= eTplus * U
  nothing
end
function greens2effective_greens!(mc::DQMC_CBFalse, greens::AbstractMatrix)
  eTminus = mc.l.hopping_matrix_exp
  eTplus = mc.l.hopping_matrix_exp_inv
  tmp = mc.s.tmp

  mul!(tmp, greens, eTplus)
  mul!(greens, eTminus, tmp)
  nothing
end
function effective_greens2greens(mc::AbstractDQMC, greens::AbstractMatrix)
  g = copy(greens)
  effective_greens2greens!(mc, g)
  return g
end

# -------------------------------------------------------
#         Time-displaced Green's function
# -------------------------------------------------------
@enum Direction begin
  LEFT
  RIGHT
end

function allocate_tdgfs!(mc)
  @stackshortcuts
  M = mc.p.slices
  Nflv = N*flv
  meas = mc.s.meas

  nranges = length(mc.s.ranges)

  meas.Gt0 = Matrix{G}[zeros(G, Nflv, Nflv) for _ in 1:M]
  meas.G0t = Matrix{G}[zeros(G, Nflv, Nflv) for _ in 1:M]

  mc.s.meas.BT0Inv_u_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BT0Inv_d_stack = Vector{Float64}[zeros(Float64, flv*N) for _ in 1:nranges]
  mc.s.meas.BT0Inv_t_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BBetaT_u_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BBetaT_d_stack = Vector{Float64}[zeros(Float64, flv*N) for _ in 1:nranges]
  mc.s.meas.BBetaT_t_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BT0_u_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BT0_d_stack = Vector{Float64}[zeros(Float64, flv*N) for _ in 1:nranges]
  mc.s.meas.BT0_t_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BBetaTInv_u_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]
  mc.s.meas.BBetaTInv_d_stack = Vector{Float64}[zeros(Float64, flv*N) for _ in 1:nranges]
  mc.s.meas.BBetaTInv_t_stack = Matrix{G}[zeros(G, flv*N, flv*N) for _ in 1:nranges]

  println("Allocated memory for TDGF measurement.")
  nothing
end

function deallocate_tdgfs_stacks!(mc)

  mc.s.meas.BT0Inv_u_stack = Matrix{G}[]
  mc.s.meas.BT0Inv_d_stack = Vector{Float64}[]
  mc.s.meas.BT0Inv_t_stack = Matrix{G}[]
  mc.s.meas.BBetaT_u_stack = Matrix{G}[]
  mc.s.meas.BBetaT_d_stack = Vector{Float64}[]
  mc.s.meas.BBetaT_t_stack = Matrix{G}[]
  mc.s.meas.BT0_u_stack = Matrix{G}[]
  mc.s.meas.BT0_d_stack = Vector{Float64}[]
  mc.s.meas.BT0_t_stack = Matrix{G}[]
  mc.s.meas.BBetaTInv_u_stack = Matrix{G}[]
  mc.s.meas.BBetaTInv_d_stack = Vector{Float64}[]
  mc.s.meas.BBetaTInv_t_stack = Matrix{G}[]

  println("Deallocated UDT stacks memory of TDGF measurement.")
  nothing
end

"""
Returns the actual memory usage of fields related to `measure_tdgfs` (in MB).
"""
function memory_usage_tdgfs(mc)
  s = x -> Base.summarysize(x) / 1024 / 1024 # size in MB
  m = mc.s.meas

  mem = 0.

  mem += s(m.Gt0)
  mem += s(m.G0t)

  mem += s(m.BT0Inv_u_stack)
  mem += s(m.BT0Inv_d_stack)
  mem += s(m.BT0Inv_t_stack)
  mem += s(m.BBetaT_u_stack)
  mem += s(m.BBetaT_d_stack)
  mem += s(m.BBetaT_t_stack)
  mem += s(m.BT0_u_stack)
  mem += s(m.BT0_d_stack)
  mem += s(m.BT0_t_stack)
  mem += s(m.BBetaTInv_u_stack)
  mem += s(m.BBetaTInv_d_stack)
  mem += s(m.BBetaTInv_t_stack)

  return mem
end

"""
    estimate_memory_usage_tdgfs(L, beta; flv=4, safe_mult=10, delta_tau=0.1)

Estimates the expected memory usage of fields related to `measure_tdgfs` (in MB).
"""
function estimate_memory_usage_tdgfs(L, beta; flv=4, safe_mult=10, delta_tau=0.1)
  N = L^2
  M = Int(beta / delta_tau)

  gfmem = (N * flv)^2 * 128 / 8 # gfmem in bytes

  mem = 2 * gfmem * M # Gt0, G0t
  mem += 8 * gfmem * M / safe_mult # u and t stacks
  mem += 4 * (N * 64 / 8) # d_stack

  return round(mem / 1024 / 1024, digits=1)
end

"""
Estimates the expected memory usage of fields related to `measure_tdgfs` (in MB).
"""
estimate_memory_usage_tdgfs(mc) = estimate_memory_usage_tdgfs(mc.p.L, mc.p.beta;
                                                                flv=mc.p.flv,
                                                                safe_mult=mc.p.safe_mult,
                                                                delta_tau=mc.p.delta_tau)

# TODO: Comment!
function measure_tdgfs!(mc)
  G = geltype(mc)
  M = mc.p.slices
  safe_mult = mc.p.safe_mult

  Gt0 = mc.s.meas.Gt0
  G0t = mc.s.meas.G0t

  BT0Inv_u_stack = mc.s.meas.BT0Inv_u_stack
  BT0Inv_d_stack = mc.s.meas.BT0Inv_d_stack
  BT0Inv_t_stack = mc.s.meas.BT0Inv_t_stack
  BBetaT_u_stack = mc.s.meas.BBetaT_u_stack
  BBetaT_d_stack = mc.s.meas.BBetaT_d_stack
  BBetaT_t_stack = mc.s.meas.BBetaT_t_stack
  BT0_u_stack = mc.s.meas.BT0_u_stack
  BT0_d_stack = mc.s.meas.BT0_d_stack
  BT0_t_stack = mc.s.meas.BT0_t_stack
  BBetaTInv_u_stack = mc.s.meas.BBetaTInv_u_stack
  BBetaTInv_d_stack = mc.s.meas.BBetaTInv_d_stack
  BBetaTInv_t_stack = mc.s.meas.BBetaTInv_t_stack

  # ---- first, calculate Gt0 and G0t only at safe_mult slices
  # right mult (Gt0)
  calc_Bchain_udts!(mc, BT0Inv_u_stack, BT0Inv_d_stack, BT0Inv_t_stack, invert=true, dir=LEFT);
  calc_Bchain_udts!(mc, BBetaT_u_stack, BBetaT_d_stack, BBetaT_t_stack, invert=false, dir=RIGHT);

  # left mult (G0t)
  calc_Bchain_udts!(mc, BT0_u_stack, BT0_d_stack, BT0_t_stack, invert=false, dir=LEFT);
  calc_Bchain_udts!(mc, BBetaTInv_u_stack, BBetaTInv_d_stack, BBetaTInv_t_stack, invert=true, dir=RIGHT);

  safe_mult_taus = 1:safe_mult:mc.p.slices
  @inbounds for i in 1:length(safe_mult_taus) # i = ith safe mult time slice
    tau = safe_mult_taus[i] # tau = tauth (overall) time slice
    if i != 1
      # Gt0
      inv_sum_udts_scalettar!(mc, Gt0[tau], BT0Inv_u_stack[i-1], BT0Inv_d_stack[i-1], BT0Inv_t_stack[i-1],
                   BBetaT_u_stack[i], BBetaT_d_stack[i], BBetaT_t_stack[i]) # G(i,0) = G(mc.s.ranges[i][1], 0), i.e. G(21, 1) for i = 3
      effective_greens2greens!(mc, Gt0[tau])

      # G0t
      inv_sum_udts_scalettar!(mc, G0t[tau], BT0_u_stack[i-1], BT0_d_stack[i-1], BT0_t_stack[i-1],
                   BBetaTInv_u_stack[i], BBetaTInv_d_stack[i], BBetaTInv_t_stack[i]) # G(i,0) = G(mc.s.ranges[i][1], 0), i.e. G(21, 1) for i = 3
      effective_greens2greens!(mc, G0t[tau])
    else
      # Gt0
      inv_one_plus_udt_scalettar!(mc, Gt0[tau], BBetaT_u_stack[1], BBetaT_d_stack[1], BBetaT_t_stack[1])
      effective_greens2greens!(mc, Gt0[tau])

      # G0t
      inv_one_plus_udt_scalettar!(mc, G0t[tau], BBetaTInv_u_stack[1], BBetaTInv_d_stack[1], BBetaTInv_t_stack[1])
      effective_greens2greens!(mc, G0t[tau]) # TODO: check analytically that we can still do this
    end
  end

  # ---- fill time slices between safe_mult slices
  fill_tdgf!(mc, Gt0, G0t)

  @inbounds for i in 1:M
    G0t[i] .*= -1 # minus sign
  end

  nothing
end

"""
Calculate UDTs at safe_mult time slices of
dir = LEFT:
inv=false:  B(tau, 1) = B(tau) * B(tau-1) * ... * B(1)                    # mult left, 1:tau
inv=true:   [B(tau, 1)]^-1 = B(1)^-1 * B(2)^-1 * ... B(tau)^-1            # mult inv right, 1:tau

udv[i] = from 1 to mc.s.ranges[i][end]

dir = RIGHT:
inv=false:  B(beta, tau) = B(beta) * B(beta-1) * ... * B(tau)             # mult right, beta:tau
inv=true:   [B(beta, tau)]^-1 = B(tau)^-1 * B(tau+1)^-1 * ... B(beta)^-1  # mult inv left, beta:tau

udv[i] = from mc.s.ranges[i][1] to mc.p.slices (beta)
"""
function calc_Bchain_udts!(mc::AbstractDQMC, u_stack, d_stack, t_stack; invert::Bool=false, dir::Direction=LEFT)
  G = geltype(mc)
  flv = mc.p.flv
  N = mc.l.sites
  curr_U_or_T = mc.s.curr_U
  ranges = mc.s.ranges

  rightmult = false
  ((dir == RIGHT && !invert) || (dir == LEFT && invert)) && (rightmult = true)

  # @show rightmult

  rng = 1:length(ranges)
  dir == RIGHT && (rng = reverse(rng))

  # Calculate udv[i], given udv[i-1]
  @inbounds for (i, rngidx) in enumerate(rng)

    if i == 1
      copyto!(curr_U_or_T, I)
    else
      if !rightmult
        copyto!(curr_U_or_T, u_stack[i-1])
      else
        copyto!(curr_U_or_T, t_stack[i-1])
      end
    end

    slice_range = dir == RIGHT ? reverse(ranges[rngidx]) : ranges[rngidx]

    for slice in slice_range
      if invert == false
        if dir == LEFT
          multiply_B_left!(mc, slice, curr_U_or_T)
        else
          # rightmult
          multiply_B_right!(mc, slice, curr_U_or_T)
        end
      else
        if dir == LEFT
          # rightmult
          multiply_B_inv_right!(mc, slice, curr_U_or_T)
        else
          multiply_B_inv_left!(mc, slice, curr_U_or_T)
        end
      end
    end

    if i != 1
      if !rightmult
        rmul!(curr_U_or_T, Diagonal(d_stack[i-1]))
      else
        lmul!(Diagonal(d_stack[i-1]), curr_U_or_T)
      end
    end

    if !rightmult
      u_stack[i], T = decompose_udt!(curr_U_or_T, d_stack[i])
    else
      U, t_stack[i] = decompose_udt!(curr_U_or_T, d_stack[i])
    end

    if i == 1
      if !rightmult
        mul!(t_stack[i], T, I)
      else
        mul!(u_stack[i], I, U)
      end
    else
      if !rightmult
        mul!(t_stack[i], T, t_stack[i-1])
      else
        mul!(u_stack[i], u_stack[i-1], U)
      end
    end
  end

  if dir == RIGHT
    reverse!(u_stack), reverse!(d_stack), reverse!(t_stack)
  end

  nothing
end

# Given Gt0 and G0t at safe mult slices (mc.s.ranges[i][1])
# propagate to all other slices.
function fill_tdgf!(mc, Gt0, G0t)
  safe_mult = mc.p.safe_mult
  M = mc.p.slices

  safe_mult_taus = 1:safe_mult:M
  Mhalf = Int(M/2)+1
  @inbounds for tau in Mhalf:M
    (tau in safe_mult_taus) && continue # skip safe mult taus

    Gt0[tau] .= Gt0[tau-1] # copy
    multiply_B_left!(mc, tau, Gt0[tau])

    G0t[tau] .= G0t[tau-1] # copy
    multiply_B_inv_right!(mc, tau, G0t[tau])
  end

  @inbounds for tau in (Mhalf-1):-1:1
    (tau in safe_mult_taus) && continue # skip safe mult taus

    Gt0[tau] .= Gt0[tau+1] # copy
    multiply_B_inv_left!(mc, tau+1, Gt0[tau])
    # sec. arg: tau seems to sometimes give better results?!?!

    G0t[tau] .= G0t[tau+1] # copy
    multiply_B_right!(mc, tau+1, G0t[tau])
  end

  nothing
end

# Calculate "G(tau, 0)", i.e. G(slice,1) as G(slice,1) = [B(slice, 1)^-1 + B(beta, slice)]^-1 which is equal to B(slice,1)G(1)
function calc_tdgf_direct(mc::AbstractDQMC, slice::Int, safe_mult::Int=mc.p.safe_mult; scalettar=true)
  if slice != 1
    Ul, Dl, Tl = calc_Bchain_inv(mc, 1, slice-1, safe_mult)
  else
    Ul, Dl, Tl = mc.s.eye_full, mc.s.ones_vec, mc.s.eye_full
  end

  if slice != mc.p.slices
    Ur, Dr, Tr = calc_Bchain(mc, slice, mc.p.slices, safe_mult)
  else
    Ur, Dr, Tr = mc.s.eye_full, mc.s.ones_vec, mc.s.eye_full
  end

  # time displace
  if !scalettar
    U, D, T = inv_sum_udts(Ul, Dl, Tl, Ur, Dr, Tr)
  else
    U, D, T = inv_sum_udts_scalettar(Ul, Dl, Tl, Ur, Dr, Tr)
  end
  effective_greens2greens!(mc, U, T)

  rmul!(U, Diagonal(D))
  return U*T
end

# function test_stacks()
#   nr = length(mc.s.ranges)

#   check_unitarity(BT0_u_stack)
#   check_unitarity(BBetaTInv_u_stack)
#   check_unitarity(BT0Inv_u_stack)
#   check_unitarity(BBetaT_u_stack)

#   # test left multiplications
#   B2 = BT0_u_stack[3] * Diagonal(BT0_d_stack[3]) * BT0_t_stack[3];
#   U, D, T = calc_Bchain(mc, 1, mc.s.ranges[3][end]); B1 = U*Diagonal(D)*T;
#   compare(B1, B2) # this is exactly the same

#   B2 = BBetaTInv_u_stack[3] * Diagonal(BBetaTInv_d_stack[3]) * BBetaTInv_t_stack[3];
#   U, D, T = calc_Bchain_inv(mc, mc.s.ranges[3][1], mc.p.slices); B1 = U*Diagonal(D)*T;
#   compare(B1, B2) # why is there a difference here at all?

#   # test right multiplications
#   B2 = BT0Inv_u_stack[3] * Diagonal(BT0Inv_d_stack[3]) * BT0Inv_t_stack[3];
#   U, D, T = calc_Bchain_inv(mc, 1, mc.s.ranges[3][end]); B1 = U*Diagonal(D)*T;
#   compare(B1, B2)

#   B2 = BBetaT_u_stack[3] * Diagonal(BBetaT_d_stack[3]) * BBetaT_t_stack[3];
#   U, D, T = calc_Bchain(mc, mc.s.ranges[3][1], mc.p.slices); B1 = U*Diagonal(D)*T;
#   compare(B1, B2)

#   # compare B(beta,1) from BT0 and BBetaT
#   BT0_full = BT0_u_stack[end] * Diagonal(BT0_d_stack[end]) * BT0_t_stack[end];
#   BBetaT_full = BBetaT_u_stack[1] * Diagonal(BBetaT_d_stack[1]) * BBetaT_t_stack[1];
#   U, D, T = calc_Bchain(mc, 1, mc.s.ranges[end][end]); BBeta0 = U*Diagonal(D)*T;
#   compare(BT0_full, BBeta0)
#   compare(BBetaT_full, BBeta0) # we have (large) abs errors here. maybe it's still ok

#   # compare resulting greens
#   gT0_full = inv_one_plus_udt(BT0_u_stack[end], BT0_d_stack[end], BT0_t_stack[end])
#   gBetaT_full = inv_one_plus_udt(BBetaT_u_stack[1], BBetaT_d_stack[1], BBetaT_t_stack[1])
#   gBeta0 = calc_greens(mc, 1)
#   compare(gT0_full, gBeta0) # 1e-16
#   compare(gBetaT_full, gBeta0) # 1e-16

#   # compare B(beta, 1), build by combining BBetaT and BT0
#   i = floor(Int, nr/2)

#   U = BBetaT_t_stack[i+1] * BT0_u_stack[i]
#   rmul!(U, Diagonal(BT0_d_stack[i]))
#   lmul!(Diagonal(BBetaT_d_stack[i+1]), U)
#   u,d,t = decompose_udt(U)
#   u = BBetaT_u_stack[i+1] * u
#   t = t * BT0_t_stack[i]
#   Bfull = u * Diagonal(d) * t
#   Bfull_d = copy(d)
#   U, D, T = calc_Bchain(mc, 1, mc.p.slices); Bfull2 = U * Diagonal(D) * T;
#   Bfull2_d = copy(D)
#   compare(Bfull, Bfull2) # max absdiff: 7.7e+03, max reldiff: 3.4e+01
#   compare(Bfull_d, Bfull2_d) # max absdiff: 9.2e+03, max reldiff: 2.3e-13

#   # compare resulting greens
#   g1 = inv_one_plus_udt(u, d, t)
#   g2 = calc_greens(mc, 1)
#   compare(g1, g2) # 1e-15
# end

# function check_unitarity(u_stack)
#   for i in 1:length(u_stack)
#     U = u_stack[i]
#     !isapprox(U * adjoint(U), I) && (return false) # I was eye(U)
#   end
#   return true
# end

# SVD STUFF ----------------------------------------------------
"""
SVD DECOMPOSITION: Calculate effective(!) Green's function (direct, i.e. without stack)
"""
function calc_Bchain_udv(mc::AbstractDQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  # Returns: tuple of results (U, D, and V) and log singular values of the intermediate products
  flv = mc.p.flv
  N = mc.l.sites
  G = geltype(mc)

  @assert 0 < start <= mc.p.slices
  @assert 0 < stop <= mc.p.slices
  @assert start <= stop

  U = Matrix{G}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  Vt = Matrix{G}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in start:stop
    if mod(k,safe_mult) == 0 || k == stop # always decompose in the end
      multiply_B_left!(mc,k,U)
      rmul!(U, Diagonal(D))
      U, D, Vtnew = decompose_udv(U)
      mul!(mc.s.tmp, Vtnew, Vt)
      Vt .=  mc.s.tmp
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_B_left!(mc,k,U)
    end
  end
  return (U,D,Vt,svs)
end

function calc_Bchain_inv_udv(mc::AbstractDQMC, start::Int, stop::Int, safe_mult::Int=mc.p.safe_mult)
  # Calculate Ul, Dl, Vtl = [B(stop) ... B(start)]^(-1) = B(start)^(-1) ... B(stop)^(-1)
  flv = mc.p.flv
  slices = mc.p.slices
  N = mc.l.sites
  G = geltype(mc)

  @assert 0 < start <= slices
  @assert 0 < stop <= slices
  @assert start <= stop

  U = Matrix{G}(I, flv*N, flv*N)
  D = ones(Float64, flv*N)
  Vt = Matrix{G}(I, flv*N, flv*N)
  Vtnew = Matrix{G}(I, flv*N, flv*N)

  svs = zeros(flv*N,length(start:stop))
  svc = 1
  for k in reverse(start:stop)
    if mod(k,safe_mult) == 0
      multiply_B_inv_left!(mc,k,U)
      U *= Diagonal(D)
      U, D, Vtnew = decompose_udv!(U)
      # not yet in-place
      Vt =  Vtnew * Vt
      svs[:,svc] = log.(D)
      svc += 1
    else
      multiply_B_left!(mc,k,U)
    end
  end
  return (U,D,Vt,svs)
end

# Calculate G(slice) = [1+B(slice-1)...B(1)B(M) ... B(slice)]^(-1) and its logdet in a stable manner
function calc_greens_and_logdet_udv(mc::AbstractDQMC, slice::Int, safe_mult::Int=mc.p.safe_mult)
  flv = mc.p.flv
  N = mc.l.sites
  G = geltype(mc)

  # Calculate Ur,Dr,Vtr=B(M) ... B(slice)
  Ur, Dr, Vtr = calc_Bchain_udv(mc,slice,mc.p.slices,safe_mult)

  # Calculate Ul,Dl,Vtl=B(slice-1) ... B(1)
  if slice-1 >= 1
    Ul, Dl, Vtl = calc_Bchain_udv(mc,1,slice-1,safe_mult)
  else
    Ul = Matrix{G}(I, flv*N, flv*N)
    Dl = ones(Float64, flv*N)
    Vtl = Matrix{G}(I, flv*N, flv*N)
  end

  # Calculate Greens function
  tmp = Vtl * Ur
  inner = adjoint(Vtr * Ul) + Diagonal(Dl) * tmp * Diagonal(Dr)
  I = decompose_udv!(inner)
  U = adjoint(I[3] * Vtr)
  D = Diagonal(1 ./ I[2])
  Vt = adjoint(Ul * I[1])
  return U*D*Vt, sum(log.(diag(D)))
end

"""
Calculate ETGF (slice=1) as inv_one_plus_udv(calc_Bchain_udv)
"""
function calc_greens_udv_simple(mc::AbstractDQMC)
  Budv = calc_Bchain_udv(mc, 1, mc.p.slices);
  gudv = inv_one_plus_udv_scalettar(Budv[1], Budv[2], Budv[3])
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using Test

using Revise

import CureMIDI

 @testset "tick_to_ms test" begin
    @test CureMIDI.tick_to_ms(UInt(100), Int16(100), 100) == 600

end

@testset "ms_to_tick test" begin
    @test CureMIDI.ms_to_tick(600, Int16(100), 100) == 100

end


text:
using VibronicToolkit
using Test

test_S = 3
test_M = 2
test_energy = zeros(test_S, test_S)
test_energy_diag = zeros(test_S, test_S)
test_freq = zeros(test_M, test_S)
# The JSON format doesn't support different frequencies between surfaces.
test_freq_flat = zeros(test_M, test_S)
test_lin = zeros(test_M, test_S, test_S)
test_lin_diag = zeros(test_M, test_S, test_S)
test_quad = zeros(test_M, test_M, test_S, test_S)
test_quad_diag = zeros(test_M, test_M, test_S, test_S)
for s1 in 1:test_S
    for m1 in 1:test_M
        test_freq[m1, s1] = s1+m1
        test_freq_flat[m1, s1] = m1
    end
    test_energy_diag[s1, s1] = 0.01*(s1+s1) - 1.0
    for m1 in 1:test_M
        test_lin_diag[m1, s1, s1] = 0.001*(s1+s1+m1)
        for m2 in 1:test_M
            test_quad_diag[m2, m1, s1, s1] = 0.001*(s1+s1+m1+m2)
        end
    end
    for s2 in 1:test_S
        test_energy[s2, s1] = 0.01*(s1+s2) - 1.0
        for m1 in 1:test_M
            test_lin[m1, s2, s1] = 0.001*(s1+s2+m1)
            for m2 in 1:test_M
                test_quad[m2, m1, s2, s1] = 0.001*(s1+s2+m1+m2)
            end
        end
    end
end
test_params = (test_energy, test_freq, test_lin, test_quad)
test_params_diag = (test_energy_diag, test_freq, test_lin_diag, test_quad_diag)
test_params_flat = (test_energy, test_freq_flat, test_lin, test_quad)
test_params_flat_diag = (test_energy_diag, test_freq_flat, test_lin_diag, test_quad_diag)
test_params_trivial = (ones(test_S, test_S),
                       ones(test_M, test_S),
                       zeros(test_M, test_S, test_S),
                       zeros(test_M, test_M, test_S, test_S))

test_S_sampling = 5
test_energy_sampling = zeros(test_S_sampling, test_S_sampling)
test_freq_sampling = zeros(test_M, test_S_sampling)
test_lin_sampling = zeros(test_M, test_S_sampling, test_S_sampling)
test_quad_sampling = zeros(test_M, test_M, test_S_sampling, test_S_sampling)
for s1 in 1:test_S_sampling
    test_energy_sampling[s1, s1] = 0.01*(s1+s1) - 1.0
    for m1 in 1:test_M
        test_freq_sampling[m1, s1] = s1+m1
        test_lin_sampling[m1, s1, s1] = 0.001*(s1+s1+m1)
    end
end
test_params_sampling = (test_energy_sampling, test_freq_sampling, test_lin_sampling, test_quad_sampling)

test_energy_sym = zeros(test_S, test_S)
test_lin_sym = zeros(test_M, test_S, test_S)
test_lin_sym3 = zeros(test_M, test_S, test_S)
test_quad_sym = zeros(test_M, test_M, test_S, test_S)
for m1 in 1:test_M
    for s1 in 1:test_S
        s2 = m1 % 2 == 0 ? s1 : test_S+1-s1
        test_lin_sym[m1, s2, s1] = 1.0 + (s1 == 2 ? 1.0 : 0.0)
        test_lin_sym3[m1, s2, s1] = 1.0
    end
end
test_params_sym = (test_energy_sym, test_freq_flat, test_lin_sym, test_quad_sym)
test_params_sym3 = (test_energy_sym, test_freq_flat, test_lin_sym3, test_quad_sym)

@testset "VibronicToolkit" begin
    @testset "utilities.jl" begin include("utilities.jl") end

    @testset "system.jl" begin include("system.jl") end
    @testset "operators.jl" begin include("operators.jl") end

    @testset "jackknife.jl" begin include("jackknife.jl") end

    @testset "analytical.jl" begin include("analytical.jl") end
    @testset "sos.jl" begin include("sos.jl") end
    @testset "trotter.jl" begin include("trotter.jl") end

    @testset "sampling.jl" begin include("sampling.jl") end
    @testset "sampling_finite_difference.jl" begin include("sampling_finite_difference.jl") end
    @testset "sampling_primitive_thermodynamic.jl" begin include("sampling_primitive_thermodynamic.jl") end

    @testset "pes.jl" begin include("pes.jl") end

    @testset "iterative_decomposition.jl" begin include("iterative_decomposition.jl") end
end