text:
"Banking holidays for Brazil (federal holidays plus Carnival)."
struct BRSettlement <: HolidayCalendar end
const Brazil = BRSettlement

"B3 Exchange holidays (https://www.b3.com.br)."
struct BrazilExchange <: HolidayCalendar end
const BrazilBMF = BrazilExchange
const BrazilB3 = BrazilExchange

# Brazilian Banking Holidays
function isholiday(::BRSettlement, dt::Dates.Date)

    yy = Dates.year(dt)
    mm = Dates.month(dt)
    dd = Dates.day(dt)

    # Bisection
    if mm >= 8
        # Fixed holidays
        if (
                # Independencia do Brasil
                ((mm == 9) && (dd == 7))
                ||
                # Nossa Senhora Aparecida
                ((mm == 10) && (dd == 12))
                ||
                # Finados
                ((mm == 11) && (dd == 2))
                ||
                # Proclamacao da Republica
                ((mm == 11) && (dd == 15))
                ||
                # Natal
                ((mm == 12) && (dd == 25))
            )
            return true
        end
    else
        # mm < 8
        # Fixed holidays
        if (
                # Confraternizacao Universal
                ((mm == 1) && (dd == 1))
                ||
                # Tiradentes
                ((mm == 4) && (dd == 21))
                ||
                # Dia do Trabalho
                ((mm == 5) && (dd == 1))
            )
            return true
        end

        # Easter occurs up to April, so Corpus Christi will be up to July in the worst case, which is before August (mm < 8). See `test/easter-min-max.jl`.
        # Holidays based on easter date.
        dt_rata::Int = Dates.days(dt)
        e_rata::Int = easter_rata(Dates.Year(yy))

        if (
                # Segunda de Carnaval
                ( dt_rata == ( e_rata - 48 ) )
                ||
                # Terca de Carnaval
                ( dt_rata == ( e_rata - 47 ) )
                ||
                # Sexta-feira Santa
                ( dt_rata == ( e_rata -  2 ) )
                ||
                # Corpus Christi
                ( dt_rata == ( e_rata + 60 ) )
            )
            return true
        end
    end

    return false
end

function isholiday(::BrazilExchange, dt::Dates.Date)
    yy = Dates.year(dt)
    mm = Dates.month(dt)
    dd = Dates.day(dt)

    if (
        # Aniversário de São Paulo
        ( mm == 1 && dd == 25 && yy < 2022 )
        ||
        # Revolucão
        ( mm == 7 && dd == 9 && yy != 2020 && yy < 2022 )
        ||
        # Consciência Negra (since 2007)
        ( yy >= 2007 && mm == 11 && dd == 20 && yy != 2020 && yy < 2022 )
        # Christmas Eve
        ||
        ( mm == 12 && dd == 24)
        ||
        # Último dia útil do ano
        ( mm == 12 && (dd == 31 || (dd>=29 && Dates.dayofweek(dt) == Dates.Friday) ))
        ||
        # National holidays
        isholiday(Brazil(), dt)
       )
        return true
    end

    return false
end


text:
using ADCME

function test_thread_pool(a)
    test_thread_pool_ = load_op_and_grad("./build/libTestThreadPool","test_thread_pool")
    a = convert_to_tensor(Any[a], [Float64]); a = a[1]
    test_thread_pool_(a)
end

# TODO: specify your input parameters
a = 1.0
u = test_thread_pool(a)
config = tf.ConfigProto(inter_op_parallelism_threads=2, intra_op_parallelism_threads=3, device_count=Dict("CPU"=>1))
sess = Session(config = config); init(sess)
@show run(sess, u)

ext: jl


text:
immutable DBN
    layers::Vector{(String, RBM)}
end

function Base.show(io::IO, dbn::DBN)
    layer_names = [getname(dbn, i) for i=1:length(dbn.layers)]
    layer_str = join(layer_names, ",")
    print(io, "DBN($layer_str)")
end

# DBN fields may change in the future, so it's worth to work through accessors
getname(dbn::DBN, k::Int) = dbn.layers[k][1]
getmodel(dbn::DBN, k::Int) = dbn.layers[k][2]
function getmodel(dbn::DBN, name::String)
    k = findfirst(p -> p[1] == name, dbn.layers)
    return (k != 0) ? dbn.layers[k][2] : error("No layer named '$name'")
end

# DBN may be thought of as a list of named RBMs, thus, we allow to access these RBMs
# with a short syntax like dbn[k] or dbn["layer_name"]
getindex(dbn::DBN, k::Int) = getmodel(dbn, k)
getindex(dbn::DBN, name::String) = getmodel(dbn, name)
Base.length(dbn::DBN) = length(dbn.layers)
Base.endof(dbn::DBN) = length(dbn)

function mh_at_layer(dbn::DBN, batch::Array{Float64, 2}, layer::Int)
    hiddens = Array(Array{Float64, 2}, layer)
    hiddens[1] = mean_hiddens(dbn[1], batch)
    for k=2:layer
        hiddens[k] = mean_hiddens(dbn[k], hiddens[k-1])
    end
    hiddens[end]
end

function transform(dbn::DBN, X::Mat{Float64})
    return mh_at_layer(dbn, X, length(dbn))
end

function fit(dbn::DBN, X::Mat{Float64}; lr=0.1, n_iter=10, batch_size=100, n_gibbs=1)
    @assert minimum(X) >= 0 && maximum(X) <= 1
    n_samples = size(X,2)
    n_batches = int(ceil(n_samples / batch_size))
    for k = 1:length(dbn.layers)
        w_buf = zeros(size(dbn[k].W))
        for itr=1:n_iter
            info("Layer $(k), iteration $(itr)")
            for i=1:n_batches
                batch = X[:, ((i-1)*batch_size + 1):min(i*batch_size, end)]
                input = k == 1 ? batch : mh_at_layer(dbn, batch, k-1)
                fit_batch!(dbn[k], input, buf=w_buf, n_gibbs=n_gibbs)
            end
            # There are some serious memory issues in the batches loop
            # cause we never free any memory between loops
            gc()
        end
    end
end


text:
using uFDTD
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

function smoke_test()
    include("../src/main.jl")
end

smoke_test()


text:
# indicator of the L2 norm ball with given radius

export IndBallL2

"""
**Indicator of a Euclidean ball**

    IndBallL2(r=1.0)

Returns the indicator function of the set
```math
S = \\{ x : \\|x\\| \\leq r \\},
```
where ``\\|\\cdot\\|`` is the ``L_2`` (Euclidean) norm. Parameter `r` must be positive.
"""
struct IndBallL2{R <: Real} <: ProximableFunction
    r::R
    function IndBallL2{R}(r::R) where {R <: Real}
        if r <= 0
            error("parameter r must be positive")
        else
            new(r)
        end
    end
end

is_convex(f::IndBallL2) = true
is_set(f::IndBallL2) = true

IndBallL2(r::R=1.0) where {R <: Real} = IndBallL2{R}(r)

function (f::IndBallL2)(x::AbstractArray{T}) where {R <: Real, T <: RealOrComplex{R}}
    if isapprox_le(norm(x), f.r, atol=eps(R), rtol=sqrt(eps(R)))
        return R(0)
    end
    return R(Inf)
end

function prox!(y::AbstractArray{T}, f::IndBallL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    scal = f.r/norm(x)
    if scal > 1
        y .= x
        return R(0)
    end
    for k in eachindex(x)
        y[k] = scal*x[k]
    end
    return R(0)
end

fun_name(f::IndBallL2) = "indicator of an L2 norm ball"
fun_dom(f::IndBallL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::IndBallL2) = "x ↦ 0 if ||x|| ⩽ r, +∞ otherwise"
fun_params(f::IndBallL2) = "r = $(f.r)"

function prox_naive(f::IndBallL2, x::AbstractArray{T}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    normx = norm(x)
    if normx > f.r
        y = (f.r/normx)*x
    else
        y = x
    end
    return y, R(0)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

if !isdefined(@__MODULE__, :DefsFileId)
    const DefsFileId = String
else
    @warn("Skipping redefinition of DefsFileId to String")
end


text:
function getTableauVSPARK(s, σ, o, tsym, g, h, lq, lp, ω, d=Nothing)
    # α_q_1 and α_p_2 need to be conjugate symplectic
    # α_q_1 and α_p_3 need to be conjugate symplectic
    # α_q_2 and α_p_2 need to be conjugate symplectic
    # α_q_2 and α_p_3 need to be conjugate symplectic
    # -> α_q_1 and α_q_2 need to be identical
    # -> α_p_2 and α_p_3 need to be identical

    α_q_1 = lq.a
    α_q_2 = lq.a
    α_p_2 = lp.a
    α_p_3 = lp.a

    b_q_1 = lq.b
    b_q_2 = lq.b
    b_p_2 = lp.b
    b_p_3 = lp.b

    # α_p_1 is free, but determines a_q_1 and a_q_2
    α_p_1 = h.a

    # a_p_1, a_p_2 and a_p_3 are free (they are not used in VPSARKpSecondary integrators anymore)
    a_p_1 = g.a
    b_p_1 = g.b

    β_q_1 = lq.b
    β_q_2 = lq.b

    β_p_1 = lp.b
    β_p_2 = lp.b
    β_p_3 = lp.b

    # a_q_1_ij = b_q_1_j / b_p_1_i * (b_p_1_i - α_p_1_ji)
    a_q_1 = zeros(s,σ)
    for i in 1:s
        for j in 1:σ
            a_q_1[i,j] = b_q_1[j] / b_p_1[i] * (b_p_1[i] - α_p_1[j,i])
        end
    end

    # a_q_2_ij = b_q_2_j / b_p_1_i * (b_p_1_i - α_p_1_ji)
    a_q_2 = zeros(s,σ)
    for i in 1:s
        for j in 1:σ
            a_q_2[i,j] = b_q_2[j] / b_p_1[i] * (b_p_1[i] - α_p_1[j,i])
        end
    end

    a_p_2 = a_q_1
    a_p_3 = a_q_2

    a_q = (a_q_1, a_q_2)
    b_q = (b_q_1, b_q_2)
    c_q = g.c

    a_p = (a_p_1, a_p_2, a_p_3)
    b_p = (b_p_1, b_p_2, b_p_3)
    c_p = g.c

    α_q = (α_q_1, α_q_2)
    β_q = (β_q_1, β_q_2)
    γ_q = lq.c

    α_p = (α_p_1, α_p_2, α_p_3)
    β_p = (β_p_1, β_p_2, β_p_3)
    γ_p = lp.c

    coeff_q = CoefficientsSPARK(tsym, o, s, σ, a_q, b_q, c_q)
    coeff_p = CoefficientsSPARK(tsym, o, s, σ, a_p, b_p, c_p)
    coeff_q̃ = CoefficientsSPARK(tsym, o, σ, s, α_q, β_q, γ_q)
    coeff_p̃ = CoefficientsSPARK(tsym, o, σ, s, α_p, β_p, γ_p)

    TableauVSPARKsecondary(tsym, o, s, σ, coeff_q, coeff_p, coeff_q̃, coeff_p̃, ω, d)
end

function getTableauVSPARKLobIII(s, lq, lp)
    o = 2s-2
    tsym = Symbol("VSPARKLobIII")
    getTableauVSPARK(s, s, o, tsym, lq, lp, lq, lp, get_lobatto_ω_matrix(s), get_lobatto_d_vector(s))
end

function getTableauVSPARKLobIIIAB(s)
    lq = getCoefficientsLobIIIA(s)
    lp = getCoefficientsLobIIIB(s)
    getTableauVSPARKLobIII(s, lq, lp)
end

function getTableauVSPARKLobIIIC(s)
    l = getCoefficientsLobIII(s)
    getTableauVSPARKLobIII(s, l, l)
end

function getTableauVSPARKLobIIID(s)
    l = getCoefficientsLobIIID(s)
    getTableauVSPARKLobIII(s, l, l)
end

function getTableauVSPARKLobIIIE(s)
    l = getCoefficientsLobIIIE(s)
    getTableauVSPARKLobIII(s, l, l)
end

function getTableauVSPARKGLRKLobIII(s, σ, lq, lp)
    o = 2s
    g = getCoefficientsGLRK(s)
    tsym = Symbol("VSPARKGLRKLobIII")
    getTableauVSPARK(s, σ, o, tsym, g, get_lobatto_interstage_coefficients(s, σ), lq, lp, get_GLRK_ω_matrix(σ), get_lobatto_d_vector(σ))
end

function getTableauVSPARKGLRKLobIIIAB(s, σ=s+1)
    getTableauVSPARKGLRKLobIII(s, σ, getCoefficientsLobIIIA(σ), getCoefficientsLobIIIB(σ))
end

function getTableauVSPARKGLRKLobIIIC(s, σ=s+1)
    getTableauVSPARKGLRKLobIII(s, σ, getCoefficientsLobIII(σ), getCoefficientsLobIIIC(σ))
end

function getTableauVSPARKGLRKLobIIID(s, σ=s+1)
    getTableauVSPARKGLRKLobIII(s, σ, getCoefficientsLobIIID(σ), getCoefficientsLobIIID(σ))
end

function getTableauVSPARKGLRKLobIIIE(s, σ=s+1)
    getTableauVSPARKGLRKLobIII(s, σ, getCoefficientsLobIIIE(σ), getCoefficientsLobIIIE(σ))
end


text:
using LinearAlgebra
using SparseArrays

"""
Compute the thermal state `exp(-βH) / Z`, with `Z` the partition function.
"""
function thermal_state(βH::AbstractMatrix)
    σ = exp(-Matrix(βH))
    return σ / tr(σ)
end

max_stars_repo_name: keefehuang/vqt

max_issues_repo_name: keefehuang/vqt

max_forks_repo_name: keefehuang/vqt


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Base.Threads
using Core: ConcurrencyViolationError
import Base: copy

mutable struct ARefxy{T}
    @atomic x::T
    y::T
    ARefxy(x::T, y::T) where {T} = new{T}(x, y)
    ARefxy{T}(x, y) where {T} = new{T}(x, y)
    ARefxy{T}() where {T} = new{T}()
end

mutable struct Refxy{T}
    x::T
    y::T
    Refxy(x::T, y::T) where {T} = new{T}(x, y)
    Refxy{T}(x, y) where {T} = new{T}(x, y)
    Refxy{T}() where {T} = new() # unused, but sets ninitialized to 0
end

@test_throws ErrorException("invalid redefinition of constant ARefxy") @eval mutable struct ARefxy{T}
    @atomic x::T
    @atomic y::T
end
@test_throws ErrorException("invalid redefinition of constant ARefxy") @eval mutable struct ARefxy{T}
    x::T
    y::T
end
@test_throws ErrorException("invalid redefinition of constant ARefxy") @eval mutable struct ARefxy{T}
    x::T
    @atomic y::T
end
@test_throws ErrorException("invalid redefinition of constant Refxy") @eval mutable struct Refxy{T}
    x::T
    @atomic y::T
end

copy(r::Union{Refxy,ARefxy}) = typeof(r)(r.x, r.y)
function add(x::T, y)::T where {T}; x + y; end
swap(x, y) = y

let T1 = Refxy{NTuple{3,UInt8}},
    T2 = ARefxy{NTuple{3,UInt8}}
    @test sizeof(T1) == 6
    @test sizeof(T2) == 8
    @test fieldoffset(T1, 1) == 0
    @test fieldoffset(T2, 1) == 0
    @test fieldoffset(T1, 2) == 3
    @test fieldoffset(T2, 2) == 4
    @test !Base.datatype_haspadding(T1)
    @test Base.datatype_haspadding(T2)
    @test Base.datatype_alignment(T1) == 1
    @test Base.datatype_alignment(T2) == 4
end

# check that very large types are getting locks
let (x, y) = (Complex{Int128}(10, 30), Complex{Int128}(20, 40))
    ar = ARefxy(x, y)
    r = Refxy(x, y)
    @test 64 == sizeof(r) < sizeof(ar)
    @test sizeof(r) == sizeof(ar) - Int(fieldoffset(typeof(ar), 1))
end

struct PadIntA <: Number # internal padding
    a::Int8
    b::Int16
    PadIntA(x) = new(82, x)
end
struct PadIntB <: Number # external padding
    a::UInt8
    b::UInt8
    c::UInt8
    PadIntB(x) = new(x & 0xff, (x >> 8) & 0xff, (x >> 16) & 0xff)
end
primitive type Int24 <: Signed 24 end # integral padding
Int24(x::Int) = Core.Intrinsics.trunc_int(Int24, x)
Base.Int(x::PadIntB) = x.a + (Int(x.b) << 8) + (Int(x.c) << 16)
Base.:(+)(x::PadIntA, b::Int) = PadIntA(x.b + b)
Base.:(+)(x::PadIntB, b::Int) = PadIntB(Int(x) + b)
Base.:(+)(x::Int24, b::Int) = Core.Intrinsics.add_int(x, Int24(b))
Base.show(io::IO, x::PadIntA) = print(io, "PadIntA(", x.b, ")")
Base.show(io::IO, x::PadIntB) = print(io, "PadIntB(", Int(x), ")")
Base.show(io::IO, x::Int24) = print(io, "Int24(", Core.Intrinsics.zext_int(Int, x), ")")

@noinline function _test_field_operators(r)
    r = r[]
    T = typeof(getfield(r, :x))
    @test getfield(r, :x, :sequentially_consistent) === T(123_10)
    @test setfield!(r, :x, T(123_1), :sequentially_consistent) === T(123_1)
    @test getfield(r, :x, :sequentially_consistent) === T(123_1)
    @test replacefield!(r, :x, 123_1 % UInt, T(123_30), :sequentially_consistent, :sequentially_consistent) === (T(123_1), false)
    @test replacefield!(r, :x, T(123_1), T(123_30), :sequentially_consistent, :sequentially_consistent) === (T(123_1), true)
    @test getfield(r, :x, :sequentially_consistent) === T(123_30)
    @test replacefield!(r, :x, T(123_1), T(123_1), :sequentially_consistent, :sequentially_consistent) === (T(123_30), false)
    @test getfield(r, :x, :sequentially_consistent) === T(123_30)
    @test modifyfield!(r, :x, add, 1, :sequentially_consistent) === (T(123_30), T(123_31))
    @test modifyfield!(r, :x, add, 1, :sequentially_consistent) === (T(123_31), T(123_32))
    @test getfield(r, :x, :sequentially_consistent) === T(123_32)
    @test swapfield!(r, :x, T(123_1), :sequentially_consistent) === T(123_32)
    @test getfield(r, :x, :sequentially_consistent) === T(123_1)
    nothing
end
@noinline function test_field_operators(r)
    _test_field_operators(Ref(copy(r)))
    _test_field_operators(Ref{Any}(copy(r)))
    nothing
end
test_field_operators(ARefxy{Int}(123_10, 123_20))
test_field_operators(ARefxy{Any}(123_10, 123_20))
test_field_operators(ARefxy{Union{Nothing,Int}}(123_10, nothing))
test_field_operators(ARefxy{Complex{Int32}}(123_10, 123_20))
test_field_operators(ARefxy{Complex{Int128}}(123_10, 123_20))
test_field_operators(ARefxy{PadIntA}(123_10, 123_20))
test_field_operators(ARefxy{PadIntB}(123_10, 123_20))
#FIXME: test_field_operators(ARefxy{Int24}(123_10, 123_20))
test_field_operators(ARefxy{Float64}(123_10, 123_20))

@noinline function _test_field_orderings(r, x, y)
    @nospecialize x y
    r = r[]

    @test getfield(r, :x) === x
    @test_throws ConcurrencyViolationError("invalid atomic ordering") getfield(r, :x, :u)
    @test_throws ConcurrencyViolationError("getfield: atomic field cannot be accessed non-atomically") getfield(r, :x, :not_atomic)
    @test getfield(r, :x, :unordered) === x
    @test getfield(r, :x, :monotonic) === x
    @test getfield(r, :x, :acquire) === x
    @test_throws ConcurrencyViolationError("invalid atomic ordering") getfield(r, :x, :release) === x
    @test_throws ConcurrencyViolationError("invalid atomic ordering") getfield(r, :x, :acquire_release) === x
    @test getfield(r, :x, :sequentially_consistent) === x
    @test isdefined(r, :x)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") isdefined(r, :x, :u)
    @test_throws ConcurrencyViolationError("isdefined: atomic field cannot be accessed non-atomically") isdefined(r, :x, :not_atomic)
    @test isdefined(r, :x, :unordered)
    @test isdefined(r, :x, :monotonic)
    @test isdefined(r, :x, :acquire)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") isdefined(r, :x, :release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") isdefined(r, :x, :acquire_release)
    @test isdefined(r, :x, :sequentially_consistent)

    @test getfield(r, :y) === y
    @test_throws ConcurrencyViolationError("invalid atomic ordering") getfield(r, :y, :u)
    @test getfield(r, :y, :not_atomic) === y
    @test_throws ConcurrencyViolationError("getfield: non-atomic field cannot be accessed atomically") getfield(r, :y, :unordered)
    @test_throws ConcurrencyViolationError("getfield: non-atomic field cannot be accessed atomically") getfield(r, :y, :monotonic)
    @test_throws ConcurrencyViolationError("getfield: non-atomic field cannot be accessed atomically") getfield(r, :y, :acquire)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") getfield(r, :y, :release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") getfield(r, :y, :acquire_release)
    @test_throws ConcurrencyViolationError("getfield: non-atomic field cannot be accessed atomically") getfield(r, :y, :sequentially_consistent)
    @test isdefined(r, :y)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") isdefined(r, :y, :u)
    @test isdefined(r, :y, :not_atomic)
    @test_throws ConcurrencyViolationError("isdefined: non-atomic field cannot be accessed atomically") isdefined(r, :y, :unordered)
    @test_throws ConcurrencyViolationError("isdefined: non-atomic field cannot be accessed atomically") isdefined(r, :y, :monotonic)
    @test_throws ConcurrencyViolationError("isdefined: non-atomic field cannot be accessed atomically") isdefined(r, :y, :acquire)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") isdefined(r, :y, :release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") isdefined(r, :y, :acquire_release)
    @test_throws ConcurrencyViolationError("isdefined: non-atomic field cannot be accessed atomically") isdefined(r, :y, :sequentially_consistent)

    @test_throws ConcurrencyViolationError("invalid atomic ordering") setfield!(r, :x, y, :u)
    @test_throws ConcurrencyViolationError("setfield!: atomic field cannot be written non-atomically") setfield!(r, :x, y)
    @test_throws ConcurrencyViolationError("setfield!: atomic field cannot be written non-atomically") setfield!(r, :x, y, :not_atomic)
    @test getfield(r, :x) === x
    @test setfield!(r, :x, y, :unordered) === y
    @test setfield!(r, :x, y, :monotonic) === y
    @test_throws ConcurrencyViolationError("invalid atomic ordering") setfield!(r, :x, y, :acquire) === y
    @test setfield!(r, :x, y, :release) === y
    @test_throws ConcurrencyViolationError("invalid atomic ordering") setfield!(r, :x, y, :acquire_release) === y
    @test setfield!(r, :x, y, :sequentially_consistent) === y
    @test getfield(r, :x) === y

    @test_throws ConcurrencyViolationError("invalid atomic ordering") setfield!(r, :y, x, :u)
    @test_throws ConcurrencyViolationError("setfield!: non-atomic field cannot be written atomically") setfield!(r, :y, x, :unordered)
    @test_throws ConcurrencyViolationError("setfield!: non-atomic field cannot be written atomically") setfield!(r, :y, x, :monotonic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") setfield!(r, :y, x, :acquire)
    @test_throws ConcurrencyViolationError("setfield!: non-atomic field cannot be written atomically") setfield!(r, :y, x, :release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") setfield!(r, :y, x, :acquire_release)
    @test_throws ConcurrencyViolationError("setfield!: non-atomic field cannot be written atomically") setfield!(r, :y, x, :sequentially_consistent)
    @test getfield(r, :y) === y
    @test setfield!(r, :y, x) === x
    @test setfield!(r, :y, x, :not_atomic) === x
    @test getfield(r, :y) === x

    @test_throws ConcurrencyViolationError("invalid atomic ordering") swapfield!(r, :y, y, :u)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") swapfield!(r, :y, y, :unordered)
    @test_throws ConcurrencyViolationError("swapfield!: non-atomic field cannot be written atomically") swapfield!(r, :y, y, :monotonic)
    @test_throws ConcurrencyViolationError("swapfield!: non-atomic field cannot be written atomically") swapfield!(r, :y, y, :acquire)
    @test_throws ConcurrencyViolationError("swapfield!: non-atomic field cannot be written atomically") swapfield!(r, :y, y, :release)
    @test_throws ConcurrencyViolationError("swapfield!: non-atomic field cannot be written atomically") swapfield!(r, :y, y, :acquire_release)
    @test_throws ConcurrencyViolationError("swapfield!: non-atomic field cannot be written atomically") swapfield!(r, :y, y, :sequentially_consistent)
    @test swapfield!(r, :y, y, :not_atomic) === x

    @test_throws ConcurrencyViolationError("invalid atomic ordering") modifyfield!(r, :y, swap, y, :u)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") modifyfield!(r, :y, swap, y, :unordered)
    @test_throws ConcurrencyViolationError("modifyfield!: non-atomic field cannot be written atomically") modifyfield!(r, :y, swap, y, :monotonic)
    @test_throws ConcurrencyViolationError("modifyfield!: non-atomic field cannot be written atomically") modifyfield!(r, :y, swap, y, :acquire)
    @test_throws ConcurrencyViolationError("modifyfield!: non-atomic field cannot be written atomically") modifyfield!(r, :y, swap, y, :release)
    @test_throws ConcurrencyViolationError("modifyfield!: non-atomic field cannot be written atomically") modifyfield!(r, :y, swap, y, :acquire_release)
    @test_throws ConcurrencyViolationError("modifyfield!: non-atomic field cannot be written atomically") modifyfield!(r, :y, swap, y, :sequentially_consistent)
    @test modifyfield!(r, :y, swap, x, :not_atomic) === (y, x)

    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :u, :not_atomic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :unordered, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: non-atomic field cannot be written atomically") replacefield!(r, :y, y, y, :monotonic, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: non-atomic field cannot be written atomically") replacefield!(r, :y, y, y, :acquire, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: non-atomic field cannot be written atomically") replacefield!(r, :y, y, y, :release, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: non-atomic field cannot be written atomically") replacefield!(r, :y, y, y, :acquire_release, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: non-atomic field cannot be written atomically") replacefield!(r, :y, y, y, :sequentially_consistent, :not_atomic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :u)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :unordered)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :monotonic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :acquire)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :acquire_release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :y, y, y, :not_atomic, :sequentially_consistent)
    @test replacefield!(r, :y, x, y, :not_atomic, :not_atomic) === (x, true)
    @test replacefield!(r, :y, x, y, :not_atomic, :not_atomic) === (y, x === y)
    @test replacefield!(r, :y, y, y, :not_atomic) === (y, true)
    @test replacefield!(r, :y, y, y) === (y, true)

    @test_throws ConcurrencyViolationError("invalid atomic ordering") swapfield!(r, :x, x, :u)
    @test_throws ConcurrencyViolationError("swapfield!: atomic field cannot be written non-atomically") swapfield!(r, :x, x, :not_atomic)
    @test_throws ConcurrencyViolationError("swapfield!: atomic field cannot be written non-atomically") swapfield!(r, :x, x)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") swapfield!(r, :x, x, :unordered) === y
    @test swapfield!(r, :x, x, :monotonic) === y
    @test swapfield!(r, :x, x, :acquire) === x
    @test swapfield!(r, :x, x, :release) === x
    @test swapfield!(r, :x, x, :acquire_release) === x
    @test swapfield!(r, :x, x, :sequentially_consistent) === x

    @test_throws ConcurrencyViolationError("invalid atomic ordering") modifyfield!(r, :x, swap, x, :u)
    @test_throws ConcurrencyViolationError("modifyfield!: atomic field cannot be written non-atomically") modifyfield!(r, :x, swap, x, :not_atomic)
    @test_throws ConcurrencyViolationError("modifyfield!: atomic field cannot be written non-atomically") modifyfield!(r, :x, swap, x)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") modifyfield!(r, :x, swap, x, :unordered)
    @test modifyfield!(r, :x, swap, x, :monotonic) === (x, x)
    @test modifyfield!(r, :x, swap, x, :acquire) === (x, x)
    @test modifyfield!(r, :x, swap, x, :release) === (x, x)
    @test modifyfield!(r, :x, swap, x, :acquire_release) === (x, x)
    @test modifyfield!(r, :x, swap, x, :sequentially_consistent) === (x, x)

    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :u, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be written non-atomically") replacefield!(r, :x, x, x)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be written non-atomically") replacefield!(r, :x, y, x, :not_atomic, :not_atomic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :unordered, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be accessed non-atomically") replacefield!(r, :x, x, x, :monotonic, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be accessed non-atomically") replacefield!(r, :x, x, x, :acquire, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be accessed non-atomically") replacefield!(r, :x, x, x, :release, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be accessed non-atomically") replacefield!(r, :x, x, x, :acquire_release, :not_atomic)
    @test_throws ConcurrencyViolationError("replacefield!: atomic field cannot be accessed non-atomically") replacefield!(r, :x, x, x, :sequentially_consistent, :not_atomic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :u)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :unordered)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :monotonic)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :acquire)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :acquire_release)
    @test_throws ConcurrencyViolationError("invalid atomic ordering") replacefield!(r, :x, x, x, :not_atomic, :sequentially_consistent)
    @test replacefield!(r, :x, x, y, :sequentially_consistent, :sequentially_consistent) === (x, true)
    @test replacefield!(r, :x, x, y, :sequentially_consistent, :sequentially_consistent) === (y, x === y)
    @test replacefield!(r, :x, y, x, :sequentially_consistent) === (y, true)
    nothing
end
@noinline function test_field_orderings(r, x, y)
    _test_field_orderings(Ref(copy(r)), x, y)
    _test_field_orderings(Ref{Any}(copy(r)), x, y)
    nothing
end
@noinline test_field_orderings(x, y) = (@nospecialize; test_field_orderings(ARefxy(x, y), x, y))
test_field_orderings(10, 20)
test_field_orderings(true, false)
test_field_orderings("hi", "bye")
test_field_orderings(:hi, :bye)
test_field_orderings(nothing, nothing)
test_field_orderings(ARefxy{Any}(123_10, 123_20), 123_10, 123_20)
test_field_orderings(ARefxy{Any}(true, false), true, false)
test_field_orderings(ARefxy{Union{Nothing,Missing}}(nothing, missing), nothing, missing)
test_field_orderings(ARefxy{Union{Nothing,Int}}(nothing, 123_1), nothing, 123_1)
test_field_orderings(Complex{Int128}(10, 30), Complex{Int128}(20, 40))
test_field_orderings(10.0, 20.0)
test_field_orderings(NaN, Inf)

struct UndefComplex{T}
    re::T
    im::T
    UndefComplex{T}() where {T} = new{T}()
end
Base.convert(T::Type{<:UndefComplex}, S) = T()
@noinline function _test_field_undef(r)
    r = r[]
    TT = fieldtype(typeof(r), :x)
    x = convert(TT, 12345_10)
    @test_throws UndefRefError getfield(r, :x)
    @test_throws UndefRefError getfield(r, :x, :sequentially_consistent)
    @test_throws UndefRefError modifyfield!(r, :x, add, 1, :sequentially_consistent)
    @test_throws (TT === Any ? UndefRefError : TypeError) replacefield!(r, :x, 1, 1.0, :sequentially_consistent)
    @test_throws UndefRefError replacefield!(r, :x, 1, x, :sequentially_consistent)
    @test_throws UndefRefError getfield(r, :x, :sequentially_consistent)
    @test_throws UndefRefError swapfield!(r, :x, x, :sequentially_consistent)
    @test getfield(r, :x, :sequentially_consistent) === x === getfield(r, :x)
    nothing
end
@noinline function test_field_undef(TT)
    _test_field_undef(Ref(TT()))
    _test_field_undef(Ref{Any}(TT()))
    nothing
end
test_field_undef(ARefxy{BigInt})
test_field_undef(ARefxy{Any})
test_field_undef(ARefxy{Union{Nothing,Integer}})
test_field_undef(ARefxy{UndefComplex{Any}})
test_field_undef(ARefxy{UndefComplex{UndefComplex{Any}}})

@test_throws ErrorException @macroexpand @atomic foo()
@test_throws ErrorException @macroexpand @atomic foo += bar
@test_throws ErrorException @macroexpand @atomic foo += bar
@test_throws ErrorException @macroexpand @atomic foo = bar
@test_throws ErrorException @macroexpand @atomic foo()
@test_throws ErrorException @macroexpand @atomic foo(bar)
@test_throws ErrorException @macroexpand @atomic foo(bar, baz)
@test_throws ErrorException @macroexpand @atomic foo(bar, baz, bax)
@test_throws ErrorException @macroexpand @atomicreplace foo bar

# test macroexpansions
let a = ARefxy(1, -1)
    @test 1 === @atomic a.x
    @test 2 === @atomic :sequentially_consistent a.x = 2
    @test 3 === @atomic :monotonic a.x = 3
    local four = 4
    @test 4 === @atomic :monotonic a.x = four
    @test 3 === @atomic :monotonic a.x = four - 1
    @test_throws ConcurrencyViolationError @atomic :not_atomic a.x = 2
    @test_throws ConcurrencyViolationError @atomic :not_atomic a.x
    @test_throws ConcurrencyViolationError @atomic :not_atomic a.x += 1

    @test 3 === @atomic :monotonic a.x
    @test 5 === @atomic a.x += 2
    @test 4 === @atomic :monotonic a.x -= 1
    @test 12 === @atomic :monotonic a.x *= 3

    @test 12 === @atomic a.x
    @test (12, 13) === @atomic a.x + 1
    @test (13, 15) === @atomic :monotonic a.x + 2
    @test (15, 19) === @atomic a.x max 19
    @test (19, 20) === @atomic :monotonic a.x max 20
    @test_throws ConcurrencyViolationError @atomic :not_atomic a.x + 1
    @test_throws ConcurrencyViolationError @atomic :not_atomic a.x max 30

    @test 20 === @atomic a.x
    @test 20 === @atomicswap a.x = 1
    @test 1 === @atomicswap :monotonic a.x = 2
    @test_throws ConcurrencyViolationError @atomicswap :not_atomic a.x = 1

    @test 2 === @atomic a.x
    @test (2, true) === @atomicreplace a.x 2 => 1
    @test (1, false) === @atomicreplace :monotonic a.x 2 => 1
    @test (1, false) === @atomicreplace :monotonic :monotonic a.x 2 => 1
    @test_throws ConcurrencyViolationError @atomicreplace :not_atomic a.x 1 => 2
    @test_throws ConcurrencyViolationError @atomicreplace :monotonic :acquire a.x 1 => 2

    @test 1 === @atomic a.x
    xchg = 1 => 2
    @test (1, true) === @atomicreplace a.x xchg
    @test (2, false) === @atomicreplace :monotonic a.x xchg
    @test (2, false) === @atomicreplace :acquire_release :monotonic a.x xchg
    @test_throws ConcurrencyViolationError @atomicreplace :not_atomic a.x xchg
    @test_throws ConcurrencyViolationError @atomicreplace :monotonic :acquire a.x xchg
end

max_stars_repo_name: amontoison/julia

max_issues_repo_name: amontoison/julia

max_forks_repo_name: amontoison/julia


text:
include("api.jl")
include("util.jl")


text:
module Module1

f = 5

end

max_stars_repo_name: SebastianCallh/julia-baby-steps

max_issues_repo_name: SebastianCallh/julia-baby-steps

max_forks_repo_name: SebastianCallh/julia-baby-steps


text:
# ---
# title: 885. Spiral Matrix III
# id: problem885
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/spiral-matrix-iii/description/>
# hidden: true
# ---
#
# On a 2 dimensional grid with `R` rows and `C` columns, we start at `(r0, c0)`
# facing east.
#
# Here, the north-west corner of the grid is at the first row and column, and
# the south-east corner of the grid is at the last row and column.
#
# Now, we walk in a clockwise spiral shape to visit every position in this grid.
#
# Whenever we would move outside the boundary of the grid, we continue our walk
# outside the grid (but may return to the grid boundary later.)
#
# Eventually, we reach all `R * C` spaces of the grid.
#
# Return a list of coordinates representing the positions of the grid in the
# order they were visited.
#
#
#
# **Example 1:**
#
#
#
#     Input: R = 1, C = 4, r0 = 0, c0 = 0
#     Output: [[0,0],[0,1],[0,2],[0,3]]
#
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/08/24/example_1.png)
#
#
#
#
# **Example 2:**
#
#
#
#     Input: R = 5, C = 6, r0 = 1, c0 = 4
#     Output: [[1,4],[1,5],[2,5],[2,4],[2,3],[1,3],[0,3],[0,4],[0,5],[3,5],[3,4],[3,3],[3,2],[2,2],[1,2],[0,2],[4,5],[4,4],[4,3],[4,2],[4,1],[3,1],[2,1],[1,1],[0,1],[4,0],[3,0],[2,0],[1,0],[0,0]]
#
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/08/24/example_2.png)
#
#
#
#
# **Note:**
#
#   1. `1 <= R <= 100`
#   2. `1 <= C <= 100`
#   3. `0 <= r0 < R`
#   4. `0 <= c0 < C`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
mutable struct UnequalTimeStack{GT<:Number, GMT<:AbstractArray{GT}} <: AbstractDQMCStack
    # B_{n*safe_mult} ... B_1
    forward_idx::Int64
    forward_u_stack::Vector{GMT}
    forward_d_stack::Vector{Vector{Float64}}
    forward_t_stack::Vector{GMT}

    # B_{n*safe_mult+1}^† ... B_M^†
    backward_idx::Int64
    backward_u_stack::Vector{GMT}
    backward_d_stack::Vector{Vector{Float64}}
    backward_t_stack::Vector{GMT}

    # B_{n*safe_mult+1}^-1 .. B_{(n+1)*safe_mult}^-1
    inv_done::Vector{Bool}
    inv_u_stack::Vector{GMT}
    inv_d_stack::Vector{Vector{Float64}}
    inv_t_stack::Vector{GMT}

    # Greens construction
    greens::GMT
    tmp::GMT
    U::GMT
    D::Vector{Float64}
    T::GMT

    # temps for complex Greens
    # TODO rework measurements to work well with StructArrays and remove this
    complex_greens_temp1::Matrix{ComplexF64}
    complex_greens_temp2::Matrix{ComplexF64}
    complex_greens_temp3::Matrix{ComplexF64}

    # To avoid recalculating
    last_update::Int64
    last_slices::Tuple{Int64, Int64}

    function UnequalTimeStack{GT, GMT}() where {GT<:Number, GMT<:AbstractArray{GT}}
        @assert isconcretetype(GT);
        @assert isconcretetype(GMT);
        new{GT, GMT}()
    end
end

# TODO
# function Base.sizeof(::Type{UnequalTimeStack}, mc::DQMC)
#     GreensEltype =  MonteCarlo.geltype(mc)
#     N = length(mc.model.l)
#     flv = mc.model.flv
#     blocks = mc.stack.n_elements

#     sizeof(GreensEltype) * flv*N * flv*N * blocks * 4 +         # forward+backward U, T
#     sizeof(GreensEltype) * flv*N * blocks * 2 +                 # forward+backward D
#     sizeof(GreensEltype) * flv*N * flv*N * (blocks-1) * 2 +     # inv U, T
#     sizeof(GreensEltype) * flv*N * (blocks-1) +                 # inv Dr
#     5*8 + sizeof(Bool) * (blocks-1) +                           # book-keeping/skipping
#     sizeof(GreensEltype) * flv*N * flv*N * 3 +                  # greens, U, Tr
#     sizeof(GreensEltype) * flv*N                                # D
# end

function initialize_stack(mc::DQMC, s::UnequalTimeStack)
    GreensElType  = geltype(mc)
    GreensMatType = gmattype(mc)
    N = length(lattice(mc))
    flv = nflavors(field(mc))
    # mirror stack
    M = mc.stack.n_elements

    # B_{n*safe_mult} ... B_1
    s.forward_idx = 1
    s.forward_u_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:M]
    s.forward_d_stack = [zeros(Float64, flv*N) for _ in 1:M]
    s.forward_t_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:M]

    # B_{n*safe_mult+1}^† ... B_M^†
    s.backward_idx = M-1
    s.backward_u_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:M]
    s.backward_d_stack = [zeros(Float64, flv*N) for _ in 1:M]
    s.backward_t_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:M]

    # B_{n*safe_mult+1}^-1 .. B_{(n+1)*safe_mult}^-1
    s.inv_done = fill(false, M-1)
    s.inv_u_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:M-1]
    s.inv_d_stack = [zeros(Float64, flv*N) for _ in 1:M-1]
    s.inv_t_stack = [GreensMatType(undef, flv*N, flv*N) for _ in 1:M-1]

    # Greens construction
    s.greens = GreensMatType(undef, flv*N, flv*N)
    s.tmp = GreensMatType(undef, flv*N, flv*N)
    s.U = GreensMatType(undef, flv*N, flv*N)
    s.D = zeros(Float64, flv*N)
    s.T = GreensMatType(undef, flv*N, flv*N)

    # TODO rework measurements to work well with StructArrays and remove this
    if (GreensMatType <: BlockDiagonal{ComplexF64, N, CMat64} where N) || (GreensMatType <: CMat64)
        s.complex_greens_temp1 = Matrix(s.greens)
        s.complex_greens_temp2 = Matrix(s.greens)
        s.complex_greens_temp3 = Matrix(s.greens)
    end

    # To avoid recalculating
    s.last_update = -1
    s.last_slices = (-1, -1)

    # maybe skip identities?
    copyto!(s.forward_u_stack[1], I)
    s.forward_d_stack[1] .= 1
    copyto!(s.forward_t_stack[1], I)

    copyto!(s.backward_u_stack[end], I)
    s.backward_d_stack[end] .= 1
    copyto!(s.backward_t_stack[end], I)

    s
end

################################################################################
### Stack building (full & lazy)
################################################################################

@bm function build_stack(mc::DQMC, s::UnequalTimeStack)
    # forward
    @bm "forward build" begin
        @inbounds for idx in 1:length(mc.stack.ranges)
            copyto!(mc.stack.curr_U, s.forward_u_stack[idx])
            for slice in mc.stack.ranges[idx]
                multiply_slice_matrix_left!(mc, mc.model, slice, mc.stack.curr_U)
            end
            vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(s.forward_d_stack[idx]))
            udt_AVX_pivot!(
                s.forward_u_stack[idx+1], s.forward_d_stack[idx+1],
                mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv
            )
            vmul!(s.forward_t_stack[idx+1], mc.stack.tmp1, s.forward_t_stack[idx])
        end
    end

    # backward
    @bm "backward build" begin
        @inbounds for idx in length(mc.stack.ranges):-1:1
            copyto!(mc.stack.curr_U, s.backward_u_stack[idx + 1])
            for slice in reverse(mc.stack.ranges[idx])
                multiply_daggered_slice_matrix_left!(mc, mc.model, slice, mc.stack.curr_U)
            end
            vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(s.backward_d_stack[idx + 1]))
            udt_AVX_pivot!(
                s.backward_u_stack[idx], s.backward_d_stack[idx],
                mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv
            )
            vmul!(s.backward_t_stack[idx], mc.stack.tmp1, s.backward_t_stack[idx+1])
        end
    end

    # inverse
    # TODO should this multiply to U's instead?
    @bm "inverse build" begin
        @inbounds for idx in 1:length(mc.stack.ranges)
            copyto!(s.inv_t_stack[idx], I)
            for slice in reverse(mc.stack.ranges[idx])
                multiply_slice_matrix_inv_left!(mc, mc.model, slice, s.inv_t_stack[idx])
            end
            udt_AVX_pivot!(
                s.inv_u_stack[idx], s.inv_d_stack[idx], s.inv_t_stack[idx],
                mc.stack.pivot, mc.stack.tempv
            )
        end
    end

    s.inv_done .= true
    s.forward_idx = length(mc.stack.ranges)+1
    s.backward_idx = 1
    s.last_update = mc.last_sweep
    s.last_slices = (-1, -1)

    nothing
end

@bm function lazy_build_stack(mc::DQMC, s::UnequalTimeStack)
    # stack = [0, Δτ, 2Δτ, ..., β]
    if mc.last_sweep == s.last_update && all(s.inv_done) &&
        (s.forward_idx == length(mc.stack.ranges)+1) && (s.backward_idx == 1)
        return nothing
    end

    build_stack(mc, s)

    nothing
end

@bm function lazy_build_forward!(mc::DQMC, s::UnequalTimeStack, upto)
    # stack = [0, Δτ, 2Δτ, ..., β]
    if s.last_update != mc.last_sweep
        s.last_update = mc.last_sweep
        s.inv_done .= false
        s.forward_idx = 1
        s.backward_idx = length(mc.stack.ranges)+1
    end

    # @assert 1 <= s.forward_idx
    # @assert upto-1 <= length(mc.stack.ranges)

    # forward
    @bm "forward build" begin
        @inbounds for idx in s.forward_idx:upto-1
            copyto!(mc.stack.curr_U, s.forward_u_stack[idx])
            for slice in mc.stack.ranges[idx]
                multiply_slice_matrix_left!(mc, mc.model, slice, mc.stack.curr_U)
            end
            vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(s.forward_d_stack[idx]))
            udt_AVX_pivot!(
                s.forward_u_stack[idx+1], s.forward_d_stack[idx+1],
                mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv
            )
            vmul!(s.forward_t_stack[idx+1], mc.stack.tmp1, s.forward_t_stack[idx])
        end
    end

    s.forward_idx = max(upto, s.forward_idx)
    nothing
end

@bm function lazy_build_backward!(mc::DQMC, s::UnequalTimeStack, downto)
    if s.last_update != mc.last_sweep
        s.last_update = mc.last_sweep
        s.inv_done .= false
        s.forward_idx = 1
        s.backward_idx = length(mc.stack.ranges)+1
    end

    # @assert s.backward_idx-1 <= length(mc.stack.ranges)
    # @assert 1 <= downto

    # backward
    @bm "backward build" begin
        @inbounds for idx in s.backward_idx-1:-1:downto
            copyto!(mc.stack.curr_U, s.backward_u_stack[idx + 1])
            for slice in reverse(mc.stack.ranges[idx])
                multiply_daggered_slice_matrix_left!(mc, mc.model, slice, mc.stack.curr_U)
            end
            vmul!(mc.stack.tmp1, mc.stack.curr_U, Diagonal(s.backward_d_stack[idx + 1]))
            udt_AVX_pivot!(
                s.backward_u_stack[idx], s.backward_d_stack[idx],
                mc.stack.tmp1, mc.stack.pivot, mc.stack.tempv
            )
            vmul!(s.backward_t_stack[idx], mc.stack.tmp1, s.backward_t_stack[idx+1])
        end
    end

    s.backward_idx = min(downto, s.backward_idx)
end

@bm function lazy_build_inv!(mc::DQMC, s::UnequalTimeStack, from, to)
    if s.last_update != mc.last_sweep
        s.last_update = mc.last_sweep
        s.inv_done .= false
        s.forward_idx = 1
        s.backward_idx = length(mc.stack.ranges)+1
    end

    # @assert 1 <= from
    # @assert to <= length(mc.stack.ranges)

    # inverse
    @bm "inverse build" begin
        @inbounds for idx in from:to
            s.inv_done[idx] && continue
            s.inv_done[idx] = true
            copyto!(s.inv_t_stack[idx], I)
            for slice in reverse(mc.stack.ranges[idx])
                multiply_slice_matrix_inv_left!(mc, mc.model, slice, s.inv_t_stack[idx])
            end
            udt_AVX_pivot!(
                s.inv_u_stack[idx], s.inv_d_stack[idx], s.inv_t_stack[idx],
                mc.stack.pivot, mc.stack.tempv
            )
        end
    end

    nothing
end

########################################
# Calculating greens
########################################

"""
    greens(mc::DQMC, k, l)

Calculates the unequal-time Greens function at slice `k` and `l`, i.e.
`G(k <- l) = G(kΔτ <- lΔτ) = ⟨cᵢ(k⋅Δτ)cⱼ(l⋅Δτ)^†⟩` where `nslices(mc) ≥ k > l ≥ 0`.

Note that `G(0, 0) = G(nslices, nslices) = G(β, β) = G(0, 0)`.

This method requires the `UnequalTimeStack` to be available. The stack
variables `Ul`, `Dl`, `Tl`, `Ur`, `Dr`, `Tr`, `curr_U` and `tmp1` will
be overwritten (but are avaible as `output` and `temp`).
"""
@bm greens(mc::DQMC, slice1::Int64, slice2::Int64) = copy(_greens!(mc, slice1, slice2))
"""
    greens!(mc::DQMC, k, l[; output=mc.stack.greens_temp, temp=mc.stack.tmp1])
"""
@bm function greens!(
        mc::DQMC, slice1::Int64, slice2::Int64;
        output = mc.stack.greens_temp, temp = mc.stack.tmp1
    )
    _greens!(mc, slice1, slice2, output, temp)
end
"""
    _greens!(mc::DQMC, k, l, output=mc.stack.greens_temp, temp=mc.stack.tmp1)
"""
function _greens!(
        mc::DQMC, slice1::Int64, slice2::Int64,
        output = mc.stack.greens_temp, temp = mc.stack.tmp1
    )
    calculate_greens(mc, slice1, slice2)
    G = _greens!(mc, output, mc.ut_stack.greens, temp)
    return GreensMatrix(slice1, slice2, G)
end
@bm function calculate_greens(mc::DQMC, slice1::Int64, slice2::Int64)
    @assert 0 ≤ slice1 ≤ mc.parameters.slices
    @assert 0 ≤ slice2 ≤ mc.parameters.slices
    s = mc.ut_stack
    if (s.last_slices != (slice1, slice2)) || (s.last_update != mc.last_sweep)
        s.last_slices = (slice1, slice2)
        if slice1 ≥ slice2
            calculate_greens_full1!(mc, s, slice1, slice2)
        else
            calculate_greens_full2!(mc, s, slice1, slice2)
        end
    end
    s.greens
end

# Notes on ranges:
# Let's say we are looking for the first full block in a product chain
# i i+1 ... j-1 j. We can find this with _find_range_with_value(mc, i - 1) + 1
# - if i is the first element of a block ... i-1 | i ... then
#   _find_range_with_value(mc, i-1) is the block just before that
# - if i is not the first element, then i-1 is in the same block. Thus
#   _find_range_with_value(mc, i-1) == _find_range_with_value(mc, i) and the
#   first full block is _find_range_with_value(mc, i-1) + 1
# Similarly for last full block we can check _find_range_with_value(mc, i+1) - 1
function _find_range_with_value(mc, val)
    # This returns the index of the range in mc.stack.ranges which contains val.
    # If val is below or above all of the ranges it returns 0 or length(ranges)+1
    ranges = mc.stack.ranges
    if val < 1
        return 0
    elseif val > last(last(ranges))
        return length(ranges) + 1
    else
        # this might actually be exact?
        estimate = clamp(
            trunc(Int, length(mc.stack.ranges) * (val-1) / (last(last(ranges)) - 1) + 1),
            1, length(ranges)
        )
        low = first(ranges[estimate])
        high = last(ranges[estimate])
        while 1 <= estimate <= length(ranges)
            if low <= val <= high
                return estimate
            elseif val < low
                estimate -= 1
                high = low - 1
                low = first(ranges[estimate])
            else
                estimate += 1
                low = high + 1
                high = last(ranges[estimate])
            end
        end
        error("Failed to find slice index in ranges. Oh no.")
    end
end

# low = slice2, high = slice1
"""
    compute_inverse_udt_block!(mc, s::UnequalTimeStack, low, high[, ...])

Computes a inverse UDT block `U D T = B_{low+1}^-1 B_{low+2}^-1 ⋯ B_high^-1`.

Default extra args are: (in order)
* `U = s.U`
* `D = s.D`
* `T = s.T`
* `tmp1 = mc.stack.tmp1`
* `tmp2 = mc.stack.tmp2`
"""
function compute_inverse_udt_block!(
        mc, s, low, high,
        U = s.U, D = s.D, T = s.T, tmp1 = mc.stack.tmp1, tmp2 = mc.stack.tmp2
    )
    # @assert low ≤ high
    # B_{low+1}^-1 B_{low+2}^-1 ⋯ B_{high-1}^-1 B_high^-1

    # Combine pre-computed blocks
    @bm "inverse pre-computed" begin
        # lower = div(low+1 + mc.parameters.safe_mult - 2, mc.parameters.safe_mult) + 1
        lower = _find_range_with_value(mc, low) + 1
        # upper = div(high, mc.parameters.safe_mult)
        upper = _find_range_with_value(mc, high+1) - 1

        lazy_build_inv!(mc, s, lower, upper)
        copyto!(U, I)
        D .= 1
        copyto!(T, I)

        # UDT combining style
        for idx in lower:upper
            # () is operation order, [] marks UDT decomposition
            # U [(D (T stack_U)) stack_D] stack_T
            # (U [U') D (T'] stack_T)
            # U D T
            vmul!(tmp1, T, s.inv_u_stack[idx])
            vmul!(tmp2, Diagonal(D), tmp1)
            vmul!(tmp1, tmp2, Diagonal(s.inv_d_stack[idx]))
            udt_AVX_pivot!(tmp2, D, tmp1, mc.stack.pivot, mc.stack.tempv)
            vmul!(T, tmp1, s.inv_t_stack[idx])
            vmul!(tmp1, U, tmp2)
            copyto!(U, tmp1)
        end
    end

    # remaining multiplications to reach specified bounds
    @bm "inverse fine tuning" begin
        # lower_slice = (lower-1) * mc.parameters.safe_mult + 1
        lower_slice = if lower <= length(mc.stack.ranges)
            first(mc.stack.ranges[lower]) else last(last(mc.stack.ranges)) + 1
        end
        # upper_slice = upper * mc.parameters.safe_mult
        upper_slice = upper > 0 ? last(mc.stack.ranges[upper]) : 0

        for slice in min(lower_slice-1, high) : -1 : low+1
            multiply_slice_matrix_inv_left!(mc, mc.model, slice, U)
        end
        if min(lower_slice-1, high) ≥ low+1
            vmul!(tmp1, U, Diagonal(D))
            udt_AVX_pivot!(U, D, tmp1, mc.stack.pivot, mc.stack.tempv)
            vmul!(tmp2, tmp1, T)
            copyto!(T, tmp2)
        end
        for slice in max(upper_slice+1, min(lower_slice-1, high)+1) : high
            multiply_slice_matrix_inv_right!(mc, mc.model, slice, T)
        end
        # U D T = B_{l+1}^-1 B_{l+2}^-1 ⋯ B_{k-1}^-1 B_k^-1
    end
    return nothing
end

"""
    compute_forward_udt_block!(mc, s::UnequalTimeStack, slice[, ...])

Computes a UDT block `U D T = B_slice B_{slice-1} ⋯ B_1`.

Default extra args are: (in order)
* `U = mc.stack.Ul`
* `D = mc.stack.Dl`
* `T = mc.stack.Tl`
* `tmp = mc.stack.tmp1`
"""
function compute_forward_udt_block!(
        mc, s, slice, apply_pivot = Val(true),
        U = mc.stack.Ul, D = mc.stack.Dl, T = mc.stack.Tl, tmp = mc.stack.tmp1,
    )
    # B(slice, 1) = B_slice B_{slice-1} ⋯ B_1 = Ul Dl Tl
    @bm "forward B" begin
        # idx = div(slice-1, mc.parameters.safe_mult)
        idx = max(0, _find_range_with_value(mc, slice) - 1)

        lazy_build_forward!(mc, s, idx+1)
        copyto!(T, s.forward_u_stack[idx+1])
        target = idx > 0 ? last(mc.stack.ranges[idx]) + 1 : 1

        for l in target : slice
            multiply_slice_matrix_left!(mc, mc.model, l, T)
        end
        vmul!(tmp, T, Diagonal(s.forward_d_stack[idx+1]))
        udt_AVX_pivot!(U, D, tmp, mc.stack.pivot, mc.stack.tempv, apply_pivot)
        vmul!(T, tmp, s.forward_t_stack[idx+1])
    end
    return nothing
end

"""
    compute_backward_udt_block!(mc, s::UnequalTimeStack, slice[, ...])

Computes a UDT block `(U D T)^† = T^† D U^† = B_N B_{N-1} ⋯ B_slice`.

Default extra args are: (in order)
* `U = mc.stack.Ur`
* `D = mc.stack.Dr`
* `T = mc.stack.Tr`
* `tmp = mc.stack.tmp1`
"""
function compute_backward_udt_block!(
        mc, s, slice,
        U = mc.stack.Ur, D = mc.stack.Dr, T = mc.stack.Tr, tmp = mc.stack.tmp1
    )
    # B(N, slice) = B_N B_{N-1} ⋯ B_{slice+1} = (Ur Dr Tr)^† = Tr^† Dr^† Ur^†
    @bm "backward B" begin
        # idx = div.(slice + mc.parameters.safe_mult - 1, mc.parameters.safe_mult)
        idx = _find_range_with_value(mc, slice) + 1

        lazy_build_backward!(mc, s, idx)
        copyto!(U, s.backward_u_stack[idx])
        target = if idx <= length(mc.stack.ranges)
            first(mc.stack.ranges[idx]) - 1 else last(last(mc.stack.ranges))
        end

        for l in target : -1 : slice+1
            multiply_daggered_slice_matrix_left!(mc, mc.model, l, U)
        end
        vmul!(tmp, U, Diagonal(s.backward_d_stack[idx]))
        udt_AVX_pivot!(U, D, tmp, mc.stack.pivot, mc.stack.tempv)
        vmul!(T, tmp, s.backward_t_stack[idx])
    end
    return nothing
end

@bm function calculate_greens_full1!(mc, s, slice1, slice2)
    # stack = [0, Δτ, 2Δτ, ..., β] = [0, safe_mult, 2safe_mult, ... N]
    # @assert slice1 ≥ slice2

    # k ≥ l or slice1 ≥ slice2
    # B_{l+1}^-1 B_{l+2}^-1 ⋯ B_{k-1}^-1 B_k^-1
    compute_inverse_udt_block!(mc, s, slice2, slice1) # low high

    # B(slice2, 1) = Ul Dl Tl
    compute_forward_udt_block!(mc, s, slice2)

    # B(N, slice1) = (Ur Dr Tr)^† = Tr^† Dr^† Ur^†
    compute_backward_udt_block!(mc, s, slice1)

    # U D T remains valid here
    @bm "compute G" begin
        # [B_{l+1}^-1 B_{l+2}^-1 ⋯ B_k^-1 + B_l ⋯ B_1 B_N ⋯ B_{k+1}]^-1
        # [U D T + Ul (Dl Tl Tr^† Dr) Ur^†]^-1
        @bm "B1" begin
            vmul!(s.greens, mc.stack.Tl, adjoint(mc.stack.Tr))
            vmul!(mc.stack.tmp1, s.greens, Diagonal(mc.stack.Dr))
            vmul!(s.greens, Diagonal(mc.stack.Dl), mc.stack.tmp1)
        end
        # [U D T + Ul (G) Ur^†]^-1
        @bm "udt" begin
            udt_AVX_pivot!(mc.stack.Tr, mc.stack.Dr, s.greens, mc.stack.pivot, mc.stack.tempv, Val(false))
        end
        # [U D T + (Ul Tr) Dr (G Ur^†)]^-1
        @bm "B2" begin
            vmul!(mc.stack.Tl, mc.stack.Ul, mc.stack.Tr)
            # (G Ur^†) = (Ur / G)^-1
            # Ur := Ur / G
            rdivp!(mc.stack.Ur, s.greens, mc.stack.Ul, mc.stack.pivot)
        end
        # [U D T + Tl Dr Ur^-1]^-1
        # [U (D T Ur + U^† Tl Dr) Ur^-1]^-1
        # [U D_max (D_min T Ur 1/Dr_max + 1/D_max U^† Tl Dr_min) Dr_max Ur^-1]^-1
        @bm "B3" begin
            # 1/D_max U^† Tl Dr_min
            vmul!(mc.stack.Tr, adjoint(s.U), mc.stack.Tl)
            vmaxinv!(mc.stack.Dl, s.D) # mc.stack.Dl .= 1.0 ./ max.(1.0, s.D)
            vmul!(mc.stack.tmp1, Diagonal(mc.stack.Dl), mc.stack.Tr)
            vmin!(mc.stack.Dl, mc.stack.Dr) # mc.stack.Dl .= min.(1.0, mc.stack.Dr)
            vmul!(mc.stack.Tr, mc.stack.tmp1, Diagonal(mc.stack.Dl))
        end
        # [U D_max (D_min T Ur 1/Dr_max + Tr) Dr_max Ur^-1]^-1
        @bm "B4" begin
            # D_min T Ur 1/Dr_max
            vmul!(mc.stack.Tl, s.T, mc.stack.Ur)
            vmin!(mc.stack.Dl, s.D) # mc.stack.Dl .= min.(1.0, s.D)
            vmul!(mc.stack.tmp1, Diagonal(mc.stack.Dl), mc.stack.Tl)
            vmaxinv!(mc.stack.Dl, mc.stack.Dr) # mc.stack.Dl .= 1.0 ./ max.(1.0, mc.stack.Dr)
            vmul!(mc.stack.Tl, mc.stack.tmp1, Diagonal(mc.stack.Dl))
        end
        # [U D_max (Tl + Tr) Dr_max Ur^-1]^-1
        @bm "sum, UDT" begin
            rvadd!(mc.stack.Tl, mc.stack.Tr)
            udt_AVX_pivot!(mc.stack.Tr, mc.stack.Dl, mc.stack.Tl, mc.stack.pivot, mc.stack.tempv, Val(false))
        end
        # [U D_max (Tr Dl Tl) Dr_max Ur^-1]^-1
        # Ur 1/Dr_max Tl^-1 1/Dl Tr^† D_max U^†
        @bm "B5" begin
            # [[((1/Dr_max) / Tl) 1/Dl] Tr^†] D_max
            vmaxinv!(mc.stack.Dr, mc.stack.Dr) # mc.stack.Dr .= 1.0 ./ max.(1.0, mc.stack.Dr)
            copyto!(mc.stack.Ul, Diagonal(mc.stack.Dr))
            rdivp!(mc.stack.Ul, mc.stack.Tl, mc.stack.tmp1, mc.stack.pivot)
            vinv!(mc.stack.Dl) # mc.stack.Dl[i] .= 1.0 ./ mc.stack.Dl[i]
            vmul!(mc.stack.tmp1, mc.stack.Ul, Diagonal(mc.stack.Dl))
            vmul!(mc.stack.Ul, mc.stack.tmp1, adjoint(mc.stack.Tr))
            vmaxinv!(mc.stack.Dl, s.D) # mc.stack.Dl .= 1.0 ./ max.(1.0, s.D)
            vmul!(s.greens, mc.stack.Ul, Diagonal(mc.stack.Dl))
        end
        # Ur G U^†
        @bm "B6" begin
            vmul!(mc.stack.Tr, s.greens, adjoint(s.U))
            vmul!(s.greens, mc.stack.Ur, mc.stack.Tr)
        end
    end

    s.greens
end

@bm function calculate_greens_full2!(mc, s, slice1, slice2)
    # stack = [0, Δτ, 2Δτ, ..., β] = [0, safe_mult, 2safe_mult, ... N]
    # @assert slice1 ≤ slice2

    # k ≤ l or slice1 ≤ slice2
    # B_{k+1}^-1 B_{k+2}^-1 ⋯ B_{l-1}^-1 B_l^-1
    compute_inverse_udt_block!(mc, s, slice1, slice2) # low high

    # B(slice1, 1) = B_k ⋯ B_1 = Ul Dl Tl
    compute_forward_udt_block!(mc, s, slice1)

    # B(M, slice2) = B_M ⋯ B_slice2+1 = (Ur Dr Tr)^† = Tr^† Dr^† Ur^†
    compute_backward_udt_block!(mc, s, slice2)

    # U D T remains valid here
    @bm "compute G" begin
        # [B_{l} B_{l-1} ⋯ B_{k+1} + (B_k ⋯ B_1 B_M ⋯ B_{l+1})^-1]^-1
        # [T^-1 D^-1 U^† + (Ul Dl Tl Tr^† Dr Ur^†)^-1]^-1
        # [T^-1 D^-1 U^† + Ur (Dl Tl Tr^† Dr)^-1 Ul^†]^-1
        @bm "B1" begin
            vmul!(s.greens, mc.stack.Tl, adjoint(mc.stack.Tr))
            vmul!(mc.stack.tmp1, Diagonal(mc.stack.Dl), s.greens)
            vmul!(s.greens, mc.stack.tmp1, Diagonal(mc.stack.Dr))
        end
        # [T^-1 D^-1 U^† + Ur G^-1 Ul^†]^-1
        @bm "udt" begin
            udt_AVX_pivot!(mc.stack.Tr, mc.stack.Dr, s.greens, mc.stack.pivot, mc.stack.tempv, Val(false))
        end
        # [T^-1 D^-1 U^† + Ur G^-1 Dr^-1 Tr^† Ul^†]^-1
        # [T^-1 D_min^-1 (D_max^-1 U^† (Ul Tr) Dr_min + D_min T Ur G^-1 Dr_max^-1) Dr_min^-1 (Ul Tr)^†]^-1
        @bm "B2" begin
            # D_max^-1 U^† (Ul Tr) Dr_min
            vmul!(mc.stack.Tl, mc.stack.Ul, mc.stack.Tr) # keep me alive
            vmul!(mc.stack.Ul, adjoint(s.U), mc.stack.Tl)
            vmaxinv!(mc.stack.Dl, s.D)
            vmul!(s.U, Diagonal(mc.stack.Dl), mc.stack.Ul)
            vmin!(mc.stack.Dl, mc.stack.Dr)
            vmul!(mc.stack.Ul, s.U, Diagonal(mc.stack.Dl))
        end
        # [T^-1 D_min^-1 (Ul + D_min T Ur G^-1 Dr_max^-1) Dr_min^-1 Tl^†]^-1
        @bm "B3" begin
            # D_min T Ur G^-1 Dr_max^-1
            vmul!(s.U, s.T, mc.stack.Ur)
            rdivp!(s.U, s.greens, mc.stack.Ur, mc.stack.pivot)
            vmin!(mc.stack.Dl, s.D)
            vmul!(mc.stack.Ur, Diagonal(mc.stack.Dl), s.U)
            vmaxinv!(mc.stack.Dl, mc.stack.Dr)
            vmul!(mc.stack.Tr, mc.stack.Ur, Diagonal(mc.stack.Dl))
        end
        # [T^-1 D_min^-1 (Ul + Tr) Dr_min^-1 Tl^†]^-1
        @bm "sum, udt" begin
            rvadd!(mc.stack.Tr, mc.stack.Ul)
            udt_AVX_pivot!(mc.stack.Ul, mc.stack.Dl, mc.stack.Tr, mc.stack.pivot, mc.stack.tempv, Val(false))
        end
        # [T^-1 D_min^-1 Ul Dl Tr Dr_min^-1 Tl^-1]^-1
        # Tl ({[(Dr_min / Tr) Dl^-1] Ul^†} D_min) T
        @bm "B4" begin
            vmin!(mc.stack.Dr, mc.stack.Dr)
            copyto!(s.U, Diagonal(mc.stack.Dr))
            rdivp!(s.U, mc.stack.Tr, mc.stack.Ur, mc.stack.pivot)
            vinv!(mc.stack.Dl)
            vmul!(mc.stack.Ur, s.U, Diagonal(mc.stack.Dl))
            vmul!(s.U, mc.stack.Ur, adjoint(mc.stack.Ul))
            vmin!(s.D, s.D)
            vmul!(mc.stack.Ur, s.U, Diagonal(s.D))
        end
        # Tl Ur T
        @bm "B6" begin
            vmul!(mc.stack.Tr, mc.stack.Ur, s.T)
            vmul!(s.greens, mc.stack.Tl, mc.stack.Tr)
            rmul!(s.greens, -1.0)
        end
    end

    s.greens
end

################################################################################
### Iterators
################################################################################

# The method above provides the unequal-time Greens function in a stable,
# general but expensive way. The idea of iterators is to provide what is needed
# more quickly, without the UnequalTimeStack and without having to worry about
# managing stack variables.

# Currently we have:
# GreensIterator
#   - attempts to be general (but isn't yet)
#   - requires UnequalTimeStack
#   - should have high accuracy, but still faster than manual
# CombinedGreensIterator
#   - does just l = 0, k = 0:M-1, but returns G_kl and G_ll
#   - does not require UnequalTimeStack
#   - same accuracy as GreensIterator

abstract type AbstractGreensIterator end
abstract type AbstractUnequalTimeGreensIterator <: AbstractGreensIterator end

init!(::AbstractUnequalTimeGreensIterator) = nothing
"""
    verify(iterator::AbstractUnequalTimeGreensIterator[, maxerror=1e-6])

Returns true if the given `iterator` is accurate up to the given `maxerror`.

See also: (@ref)[`accuracy`]
"""
verify(it::AbstractUnequalTimeGreensIterator, maxerror=1e-6) = maximum(accuracy(it)) < maxerror

# Maybe split into multiple types?
struct GreensIterator{slice1, slice2, T <: DQMC} <: AbstractUnequalTimeGreensIterator
    mc::T
    recalculate::Int64
end

"""
    GreensIterator(mc::DQMC[, ks=Colon(), ls=0, recalculate=4mc.parameters.safe_mult])

Creates an Iterator which calculates all `[G[k <- l] for l in ls for k in ks]`
efficiently. Currently only allows `ks = :, ls::Integer`.

For stability it is necessary to fully recompute `G[k <- l]` every so often. one
can adjust the frequency of recalculations via `recalculate`. To estimate the
resulting accuracy, one may use `accuracy(GreensIterator(...))`.

This iterator requires the `UnequalTimeStack`. Iteration overwrites the
`DQMCStack` variables `curr_U`, `Ul`, `Dl`, `Tl`, `Ur`, `Dr`, `Tr`, `tmp1` and
`tmp2`, with `curr_U` acting as the output. For the iteration to remain valid,
the UnequalTimeStack must not be overwritten. As such:
- `greens!(mc)` can be called and remains valid
- `greens!(mc, slice)` can be called and remains valid
- `greens!(mc, k, l)` will break iteration but remains valid (call before iterating)
"""
function GreensIterator(
        mc::T, slice1 = Colon(), slice2 = 0, recalculate = 4mc.parameters.safe_mult
    ) where {T <: DQMC}
    GreensIterator{slice1, slice2, T}(mc, recalculate)
end
init!(it::GreensIterator) = initialize_stack(it.mc, it.mc.ut_stack)
Base.length(it::GreensIterator{:, i}) where {i} = it.mc.parameters.slices + 1 - i

# Slower, versatile version:
function Base.iterate(it::GreensIterator{:, i}) where {i}
    s = it.mc.ut_stack
    calculate_greens_full1!(it.mc, s, i, i)
    copyto!(s.T, s.greens)
    udt_AVX_pivot!(s.U, s.D, s.T, it.mc.stack.pivot, it.mc.stack.tempv)
    G = _greens!(it.mc, it.mc.stack.curr_U, s.greens, it.mc.stack.Ur)
    return (GreensMatrix(i, i, G), (i+1, i))
end
function Base.iterate(it::GreensIterator{:}, state)
    s = it.mc.ut_stack
    k, l = state
    if k > it.mc.parameters.slices
        return nothing
    elseif k % it.recalculate == 0
        # Recalculate
        calculate_greens_full1!(it.mc, s, k, l) # writes s.greens
        G = _greens!(it.mc, it.mc.stack.curr_U, s.greens, it.mc.stack.Tl)
        copyto!(s.T, s.greens)
        udt_AVX_pivot!(s.U, s.D, s.T, it.mc.stack.pivot, it.mc.stack.tempv)
        return (GreensMatrix(k, l, G), (k+1, l))
    elseif k % it.mc.parameters.safe_mult == 0
        # Stabilization
        multiply_slice_matrix_left!(it.mc, it.mc.model, k, s.U)
        vmul!(it.mc.stack.curr_U, s.U, Diagonal(s.D))
        vmul!(it.mc.stack.tmp1, it.mc.stack.curr_U, s.T)
        udt_AVX_pivot!(s.U, s.D, it.mc.stack.curr_U, it.mc.stack.pivot, it.mc.stack.tempv)
        vmul!(it.mc.stack.tmp2, it.mc.stack.curr_U, s.T)
        copyto!(s.T, it.mc.stack.tmp2)
        G = _greens!(it.mc, it.mc.stack.curr_U, it.mc.stack.tmp1, it.mc.stack.tmp2)
        return (GreensMatrix(k, l, G), (k+1, l))
    else
        # Quick advance
        multiply_slice_matrix_left!(it.mc, it.mc.model, k, s.U)
        vmul!(it.mc.stack.curr_U, s.U, Diagonal(s.D))
        vmul!(it.mc.stack.tmp1, it.mc.stack.curr_U, s.T)
        G = _greens!(it.mc, it.mc.stack.curr_U, it.mc.stack.tmp1, it.mc.stack.Ur)
        s.last_slices = (-1, -1) # for safety
        return (GreensMatrix(k, l, G), (k+1, l))
    end
end
"""
    accuracy(iterator::AbstractUnequalTimeGreensIterator)

Compares values from the given iterator to more verbose computations, returning
the maximum differences for each. This can be used to check numerical stability.
"""
function accuracy(iter::GreensIterator{:, l}) where {l}
    mc = iter.mc
    Gk0s = [deepcopy(greens(mc, k, l).val) for k in l:nslices(mc)]
    [maximum(abs.(Gk0s[i] .- G.val)) for (i, G) in enumerate(iter)]
end

"""
    CombinedGreensIterator(mc::DQMC)

Returns an iterator which efficiently calculates and returns triples
(G(0, l), G(l, 0), G(l, l)) for each iteration.

By default `l` runs from 1 (Δτ) to M (β), including both. This can be adjusted
by the keyword arguments `start = 1` and `stop = mc.parameters.slices`.

## Warning

Calling `greens(mc)`, `greens(mc, slice)`, `greens(mc, slice1, slice2)` or
similar functions will break the iteration sequence. If you need the result of
one of these, call it before starting the iteration. The result will remain
valid.

This happens because the iterator uses a bunch of temporary matrices as storage
between iterations. Unlike the functions mentioned above, the iterator does not
usually calculate greens function from scratch. Instead results from the last
iteration are used to calculate the next. Note that even with matrix
stabilization this only works for some time and therefore a full recalculation
happens every so often. This can be changed by adjusting the keyword argument
`recalculate = 2mc.paremeters.safe_mult`.

You can also set `recalculate = nothing` and pass a `max_delta = 1e-7`. In this
case the iterator will pick `recalculate` such that the difference between the
iterator and the respective results from `greens(...)` is at most `max_delta`.
Note that this is an experimental feature - the dependence of the error on the
state of the simulation has not been thoroughly investigated.

For reference - the iterator uses the stack matrices `tmp1` and `tmp2` as well
as the `UnequalTimeStack` matrix `greens` as temporary outputs. The UDT
decompositions are saved in the stack matrices `Ul`, `Dl`, `Tl`, `Ur`, `Dr`, `Tr`
and the `UnequalTimeStack` matrices `U, `D`, `T`. All of these matrices
(not outputs) need to remain valid between iterations. The stack matrix `curr_U`
and unequal time stack matrix `tmp` are used as temporary storage.
"""
struct CombinedGreensIterator <: AbstractUnequalTimeGreensIterator
    recalculate::Int
    start::Int
    stop::Int
end

function CombinedGreensIterator(
        mc::DQMC, model::Model;
        recalculate = 2mc.parameters.safe_mult, max_delta = 1e-7,
        start = 1, stop = mc.parameters.slices
    )
    if recalculate === nothing
        iter = CombinedGreensIterator(typemax(Int64), start, stop)
        recalc = estimate_recalculate(iter, max_delta)
        CombinedGreensIterator(recalc, start, stop)
    else
        CombinedGreensIterator(recalculate, start, stop)
    end
end
function CombinedGreensIterator(
        mc; recalculate = 2mc.parameters.safe_mult, max_delta=1e-7,
        start = 1, stop = mc.parameters.slices
    )
    if recalculate === nothing
        iter = CombinedGreensIterator(typemax(Int64), start, stop)
        recalc = estimate_recalculate(iter, max_delta)
        CombinedGreensIterator(recalc, start, stop)
    else
        CombinedGreensIterator(recalculate, start, stop)
    end
end
init!(it::CombinedGreensIterator) = error("Woops")
Base.length(it::CombinedGreensIterator) = it.stop - it.start + 1 # both included
function Base.:(==)(a::CombinedGreensIterator, b::CombinedGreensIterator)
    (a.recalculate == b.recalculate) && (a.start == b.start) && (a.stop == b.stop)
end

struct _CombinedGreensIterator{T <: DQMC}
    mc::T
    spec::CombinedGreensIterator
    # TODO rework measurements to work well with StructArrays and remove this
    copy_to_array::Bool
end

function init(mc::DQMC, it::CombinedGreensIterator)
    _CombinedGreensIterator(mc, it, isdefined(mc.stack, :complex_greens_temp))
end
# Base.iterate(it::CombinedGreensIterator, mc::DQMC) = iterate(_CombinedGreensIterator(it, mc))
# function Base.iterate(it::CombinedGreensIterator, mc::DQMC, idx)
#     iterate(_CombinedGreensIterator(it, mc), idx)
# end

Base.length(it::_CombinedGreensIterator) = length(it.spec)

# Fast specialized version
@bm function Base.iterate(it::_CombinedGreensIterator)
    s = it.mc.stack
    uts = it.mc.ut_stack

    # Need full built stack for curr_U to be save
    build_stack(it.mc, uts)
    # force invalidate stack because we modify uts.greens
    uts.last_slices = (-1, -1)

    if it.spec.start in (0, 1)
        # Get G00
        if current_slice(it.mc) == 1
            # The measure system triggers here to be a bit more performant
            copyto!(s.Tl, s.greens)
        else
            # calculate G00
            calculate_greens_full1!(it.mc, it.mc.ut_stack, 0, 0)
            copyto!(s.Tl, uts.greens)
        end
        copyto!(s.tmp1, s.Tl)

        # Prepare next iteration
        # G01 = (G00-I)B_1^-1; then G0l+1 = G0l B_{l+1}^-1
        vsub!(s.Tr, s.Tl, I)

        # prepare UDT stacks
        udt_AVX_pivot!(s.Ul, s.Dl, s.Tl, it.mc.stack.pivot, it.mc.stack.tempv)
        copyto!(uts.U, s.Ul)
        copyto!(uts.D, s.Dl)
        copyto!(uts.T, s.Tl)

        # G01 = (G00-I)B_1^-1; then G0l+1 = G0l B_{l+1}^-1
        udt_AVX_pivot!(s.Ur, s.Dr, s.Tr, it.mc.stack.pivot, it.mc.stack.tempv)

        if it.spec.start == 0
            # return for G00 iteration
            Gll = _greens!(it.mc, uts.greens, s.tmp1, s.tmp2)
            Gl0 = copyto!(s.tmp1, uts.greens)
            G0l = copyto!(s.tmp2, uts.greens)
            if it.copy_to_array
                return ((
                    GreensMatrix(0, 0, copyto!(uts.complex_greens_temp1, G0l)),
                    GreensMatrix(0, 0, copyto!(uts.complex_greens_temp2, Gl0)),
                    GreensMatrix(0, 0, copyto!(uts.complex_greens_temp3, Gll))
                ), 1)
            else
                return ((
                    GreensMatrix(0, 0, G0l),
                    GreensMatrix(0, 0, Gl0),
                    GreensMatrix(0, 0, Gll)
                ), 1)
            end
        elseif it.spec.start == 1
            # start with l = 1 iteration
            return iterate(it, 1)
        end
    else
        # See recalculate step below
        calculate_greens_full1!(it.mc, it.mc.ut_stack, it.spec.start, 0)
        copyto!(s.curr_U, uts.greens)
        calculate_greens_full2!(it.mc, it.mc.ut_stack, 0, it.spec.start)
        copyto!(uts.tmp, uts.greens)
        calculate_greens_full1!(it.mc, it.mc.ut_stack, it.spec.start, it.spec.start)

        copyto!(uts.T, uts.greens)
        Gll = _greens!(it.mc, uts.greens, uts.T, s.tmp2)
        udt_AVX_pivot!(uts.U, uts.D, uts.T, s.pivot, s.tempv)

        Gl0 = _greens!(it.mc, s.tmp1, s.curr_U, s.tmp2)
        copyto!(s.Tl, s.curr_U)
        udt_AVX_pivot!(s.Ul, s.Dl, s.Tl, s.pivot, s.tempv)

        G0l = _greens!(it.mc, s.tmp2, uts.tmp, s.curr_U)
        copyto!(s.Tr, uts.tmp)
        udt_AVX_pivot!(s.Ur, s.Dr, s.Tr, s.pivot, s.tempv)

        if it.copy_to_array
            l = it.spec.start
            return ((
                GreensMatrix(0, l, copyto!(uts.complex_greens_temp1, G0l)),
                GreensMatrix(l, 0, copyto!(uts.complex_greens_temp2, Gl0)),
                GreensMatrix(l, l, copyto!(uts.complex_greens_temp3, Gll))
            ), l+1)
        else
            return ((
                GreensMatrix(0, it.spec.start, G0l),
                GreensMatrix(it.spec.start, 0, Gl0),
                GreensMatrix(it.spec.start, it.spec.start, Gll)
            ), it.spec.start+1)
        end
    end
end
# probably need extra temp variables
@bm function Base.iterate(it::_CombinedGreensIterator, l)
    # l is 1-based
    s = it.mc.stack
    uts = it.mc.ut_stack
    # force invalidate stack because we modify uts.greens
    uts.last_slices = (-1, -1)
    # if start is 1 we need to stabilize like if it were 0
    shift = (it.spec.start != 1) * it.spec.start

    if l > it.spec.stop
        return nothing
    elseif (l - shift) % it.spec.recalculate == 0
        # Recalculation will overwrite
        # Ul, Dl, Tl, Ur, Dr, Tr, U, D, T, tmp1, tmp2, greens
        # curr_U is used only if the stack is rebuilt
        # this leaves uts.tmp (and curr_U)

        calculate_greens_full1!(it.mc, it.mc.ut_stack, l, 0)
        copyto!(s.curr_U, uts.greens)
        calculate_greens_full2!(it.mc, it.mc.ut_stack, 0, l)
        copyto!(uts.tmp, uts.greens)
        calculate_greens_full1!(it.mc, it.mc.ut_stack, l, l)

        # input and output can be the same here
        copyto!(uts.T, uts.greens)
        Gll = _greens!(it.mc, uts.greens, uts.T, s.tmp2)
        udt_AVX_pivot!(uts.U, uts.D, uts.T, s.pivot, s.tempv)

        Gl0 = _greens!(it.mc, s.tmp1, s.curr_U, s.tmp2)
        copyto!(s.Tl, s.curr_U)
        udt_AVX_pivot!(s.Ul, s.Dl, s.Tl, s.pivot, s.tempv)

        G0l = _greens!(it.mc, s.tmp2, uts.tmp, s.curr_U)
        copyto!(s.Tr, uts.tmp)
        udt_AVX_pivot!(s.Ur, s.Dr, s.Tr, s.pivot, s.tempv)

        if it.copy_to_array
            return ((
                GreensMatrix(0, l, copyto!(uts.complex_greens_temp1, G0l)),
                GreensMatrix(l, 0, copyto!(uts.complex_greens_temp2, Gl0)),
                GreensMatrix(l, l, copyto!(uts.complex_greens_temp3, Gll))
            ), l+1)
        else
            return ((
                GreensMatrix(0, l, G0l),
                GreensMatrix(l, 0, Gl0),
                GreensMatrix(l, l, Gll)
            ), l+1)
        end

    elseif ((l - shift) % it.spec.recalculate) % it.mc.parameters.safe_mult == 0
        # Stabilization
        # Reminder: These overwrite s.tmp1 and s.tmp2
        multiply_slice_matrix_left!(it.mc, it.mc.model, l, s.Ul)
        multiply_slice_matrix_inv_right!(it.mc, it.mc.model, l, s.Tr)
        multiply_slice_matrix_left!(it.mc, it.mc.model, l, uts.U)
        multiply_slice_matrix_inv_right!(it.mc, it.mc.model, l, uts.T)

        # Gl0
        vmul!(s.tmp1, s.Ul, Diagonal(s.Dl))
        vmul!(s.tmp2, s.tmp1, s.Tl)
        udt_AVX_pivot!(s.Ul, s.Dl, s.tmp1, s.pivot, s.tempv)
        vmul!(s.curr_U, s.tmp1, s.Tl)
        copyto!(s.Tl, s.curr_U)
        Gl0 = _greens!(it.mc, s.tmp1, s.tmp2, s.curr_U)
        # tmp1 must not change anymore

        # G0l
        vmul!(s.curr_U, Diagonal(s.Dr), s.Tr)
        vmul!(uts.greens, s.Ur, s.curr_U)
        copyto!(s.Tr, s.curr_U)
        udt_AVX_pivot!(s.tmp2, s.Dr, s.Tr, s.pivot, s.tempv)
        vmul!(s.curr_U, s.Ur, s.tmp2)
        copyto!(s.Ur, s.curr_U)
        G0l = _greens!(it.mc, s.tmp2, uts.greens, s.curr_U)
        # tmp1, tmp2 must not change anymore

        # Gll
        vmul!(uts.tmp, uts.U, Diagonal(uts.D))
        vmul!(uts.greens, uts.tmp, uts.T)
        udt_AVX_pivot!(s.curr_U, uts.D, uts.tmp, it.mc.stack.pivot, it.mc.stack.tempv)
        vmul!(uts.U, uts.tmp, uts.T)
        vmul!(uts.T, Diagonal(uts.D), uts.U)
        udt_AVX_pivot!(uts.tmp, uts.D, uts.T, it.mc.stack.pivot, it.mc.stack.tempv)
        vmul!(uts.U, s.curr_U, uts.tmp)
        Gll = _greens!(it.mc, uts.greens, uts.greens, s.curr_U)

        if it.copy_to_array
            return ((
                GreensMatrix(0, l, copyto!(uts.complex_greens_temp1, G0l)),
                GreensMatrix(l, 0, copyto!(uts.complex_greens_temp2, Gl0)),
                GreensMatrix(l, l, copyto!(uts.complex_greens_temp3, Gll))
            ), l+1)
        else
            return ((
                GreensMatrix(0, l, G0l),
                GreensMatrix(l, 0, Gl0),
                GreensMatrix(l, l, Gll)
            ), l+1)
        end

    else
        # Quick advance
        multiply_slice_matrix_left!(it.mc, it.mc.model, l, s.Ul)
        # multiply_slice_matrix_inv_left!(it.mc, it.mc.model, nslices(it.mc)-l+1, s.Ur)
        multiply_slice_matrix_inv_right!(it.mc, it.mc.model, l, s.Tr)
        multiply_slice_matrix_left!(it.mc, it.mc.model, l, uts.U)
        multiply_slice_matrix_inv_right!(it.mc, it.mc.model, l, uts.T)

        # Gl0
        vmul!(s.curr_U, s.Ul, Diagonal(s.Dl))
        vmul!(s.tmp2, s.curr_U, s.Tl)
        Gl0 = _greens!(it.mc, s.tmp1, s.tmp2, s.curr_U)

        # G0l
        vmul!(s.curr_U, s.Ur, Diagonal(s.Dr))
        vmul!(uts.greens, s.curr_U, s.Tr)
        G0l = _greens!(it.mc, s.tmp2, uts.greens, s.curr_U)

        # Gll
        vmul!(s.curr_U, uts.U, Diagonal(uts.D))
        vmul!(uts.greens, s.curr_U, uts.T)
        Gll = _greens!(it.mc, uts.greens, uts.greens, s.curr_U)

        if it.copy_to_array
            return ((
                GreensMatrix(0, l, copyto!(uts.complex_greens_temp1, G0l)),
                GreensMatrix(l, 0, copyto!(uts.complex_greens_temp2, Gl0)),
                GreensMatrix(l, l, copyto!(uts.complex_greens_temp3, Gll))
            ), l+1)
        else
            return ((
                GreensMatrix(0, l, G0l),
                GreensMatrix(l, 0, Gl0),
                GreensMatrix(l, l, Gll)
            ), l+1)
        end
    end
end

function accuracy(mc::DQMC, it::CombinedGreensIterator)
    iter = _CombinedGreensIterator(mc, it)
    Gk0s = [deepcopy(greens(mc, k, 0).val) for k in 1:nslices(mc)]
    G0ks = [deepcopy(greens(mc, 0, k).val) for k in 1:nslices(mc)]
    Gkks = [deepcopy(greens(mc, k, k).val) for k in 1:nslices(mc)]
    map(enumerate(iter)) do (i, Gs)
        (
            maximum(abs.(G0ks[i] .- Gs[1].val)),
            maximum(abs.(Gk0s[i] .- Gs[2].val)),
            maximum(abs.(Gkks[i] .- Gs[3].val))
        )
    end
end

function estimate_recalculate(iter::CombinedGreensIterator, max_delta=1e-7)
    deltas = accuracy(iter)
    idx = findfirst(ds -> any(ds .> max_delta), deltas)
    if idx === nothing
        return typemax(Int)
    else
        return idx
    end
end


text:
julia> sum([4,6,8])
18

julia> +((1:10)...)
55

julia +([1,2,3]...)
6

julia> prod([4,6,8])
192

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    SourceFile(code [, filename])

A UTF-8 source code string with associated file name and indexing structures.
"""
struct SourceFile
    # We use `code::String` for now but it could be some other UTF-8 based
    # string data structure with byte-based indexing.
    #
    # For example a rope data structure may be good for incremental editing
    # https://en.wikipedia.org/wiki/Rope_(data_structure)
    code::String
    filename::Union{Nothing,String}
    # String index of start of every line
    line_starts::Vector{Int}
end

function SourceFile(code::AbstractString; filename=nothing)
    line_starts = Int[1]
    for i in eachindex(code)
        # The line is considered to start after the `\n`
        # FIXME: \r and \n\r
        code[i] == '\n' && push!(line_starts, i+1)
    end
    if isempty(code) || last(code) != '\n'
        push!(line_starts, lastindex(code)+1)
    end
    SourceFile(code, filename, line_starts)
end

function SourceFile(; filename)
    SourceFile(read(filename, String); filename=filename)
end

# Get line number of the given byte within the code
function source_line(source::SourceFile, byte_index)
    searchsortedlast(source.line_starts, byte_index)
end

"""
Get line number and character within the line at the given byte index.
"""
function source_location(source::SourceFile, byte_index)
    line = searchsortedlast(source.line_starts, byte_index)
    i = source.line_starts[line]
    column = 1
    while i < byte_index
        i = nextind(source.code, i)
        column += 1
    end
    line, column
end

"""
Get byte range of the source line at byte_index, buffered by
`context_lines_before` and `context_lines_after` before and after.
"""
function source_line_range(source::SourceFile, byte_index;
                           context_lines_before=0, context_lines_after=0)
    line = searchsortedlast(source.line_starts, byte_index)
    fbyte = source.line_starts[max(line-context_lines_before, 1)]
    lbyte = source.line_starts[min(line+1+context_lines_after, end)] - 1
    fbyte,lbyte
end

function source_location(::Type{LineNumberNode}, source::SourceFile, byte_index)
    LineNumberNode(source_line(source, byte_index), source.filename)
end

function Base.show(io::IO, ::MIME"text/plain", source::SourceFile)
    fn = isnothing(source.filename) ? "" : " $(source.filename)"
    header = "## SourceFile$fn ##"
    print(io, header, "\n")
    heightlim = displaysize(io)[1] ÷ 2
    if !get(io, :limit, false) || length(source.line_starts) <= heightlim
        print(io, source.code)
    else
        r1 = source_line_range(source, 1, context_lines_after=heightlim-3)
        print(io, view(source, r1[1]:r1[2]))
        println(io, "⋮")
    end
end

function Base.getindex(source::SourceFile, rng::AbstractRange)
    i = first(rng)
    # Convert byte range into unicode String character range.
    # Assumes valid unicode! (SubString doesn't give us a reliable way to opt
    # out of the valid unicode check. The SubString{String} inner constructor
    # has some @boundscheck, but using @inbounds depends on inlining choices.)
    j = prevind(source.code, last(rng)+1)
    source.code[i:j]
end

# TODO: Change view() here to `sourcetext` ?
function Base.view(source::SourceFile, rng::AbstractRange)
    i = first(rng)
    j = prevind(source.code, last(rng)+1)
    SubString(source.code, i, j)
end

function Base.getindex(source::SourceFile, i::Int)
    source.code[i]
end

"""
    sourcetext(source::SourceFile)

Get the full source text of a `SourceFile` as a string.
"""
function sourcetext(source::SourceFile)
    return source.code
end


text:
# tests for the UAV type and related functions.

# using Base.Test
# using GridInterpolations

# push!(LOAD_PATH, ".")
# include("includes.jl")

function get_test_uav()
    actions = [-1., 1.]
    joint_actions = ones(2, 2)
    joint_actions[:, 1] = -1

    r = collect(0:10)
    q_table = zeros(10^7, size(joint_actions, 2))
    q_grid = RectangleGrid(r, r, r, r, r, r, r)
    policy = PolicyTable(q_table, q_grid, 10^7)

    c_table = zeros(10, 10, 10, 10, 10)
    c_grid = RectangleGrid(r, r, r, r, r)
    coord = CoordinationTable(c_table, c_grid, 10^5)

    strategy = coordinated_strategy
    act_select = "best"

    rng = MersenneTwister(1)
    sensor_range_std = 15.
    sensor_theta_std = deg2rad(1.)
    sensor_psi_std = deg2rad(1.)
    sensor_velocity_std = 1.
    use_sensor_noise = false

    uav = UAV(1, actions, joint_actions, policy, coord, false, strategy, act_select, rng, sensor_range_std, sensor_theta_std, sensor_psi_std, sensor_velocity_std, use_sensor_noise)
    return uav
end

# UAV tests
function test_uav_constructor()
   get_test_uav()
end

function test_get_polar_state()
    uav = get_test_uav()
    uav_state = UAVState(ones(Float64, 8)...)
    state = State([uav_state, uav_state])
    state.responding_states = [1,0]
    polar_state = get_polar_state(uav, state)
    if uav.is_policy_joint == true
        @test polar_state == [0,0,0,sqrt(2),sqrt(2),1,0]
    else
        @test polar_state == [0,0,0,sqrt(2),sqrt(2),1]
    end
end

function test_get_belief_state()
    uav = get_test_uav()
    uav_state = UAVState(ones(Float64, 8)...)
    state = State([uav_state, uav_state])
    state.responding_states = [1,1]
    belief_state = get_belief_state(uav, state)
    @test length(belief_state) == uav.policy_table.num_states
end

function test_get_action()
    uav = get_test_uav()
    uav_state = UAVState(ones(Float64, 8)...)
    state = State([uav_state, uav_state])
    state.responding_states = [1,1]
    action, signals, _ = get_action(uav, state)
    @test action == -1.0
    @test signals == [:no_signal, :no_signal]
end
tic()
# UAV tests
test_uav_constructor()

# function tests
test_get_polar_state()
test_get_belief_state()
test_get_action()
toc()


text:
@testset "Interpolate 1D" begin
	df = DataFrame(Temp=[10,11,12,14],Humi=[40.,NaN,50,60],
	      datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
	      DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
	dfi = interpdf(df,[DateTime(2010,1,1,0,30,0),DateTime(2010,1,1,12,0,0)]);
	@test dfi[:Temp][1] ≈ 10.5
	@test isnan(dfi[:Humi][1])
	@test isnan(dfi[:Humi][2])
	@test isnan(dfi[:Temp][2])
	@test interp1([1,2,3,4],[10,20,30,40],[1.5]) ≈ [15.]

	@test interp1([DateTime(2000,1,1),DateTime(2000,1,2)],[0.0,1.0],
					[DateTime(2000,1,1,12,0,0)]) ≈ [0.5]
end

@testset "Interpolate 2D" begin
	# meshgrid
	x = [1,2,3,4];
	y = [10,20,30,40];
	xi,yi = meshgrid(x,y)
	@test xi[2,3] == 3.0
	@test yi[2,3] == 20.0

	# mesh2vec
	xi = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4];
	yi = [10 10 10 10; 20 20 20 20; 30 30 30 30; 40 40 40 40];
	x,y = mesh2vec(xi,yi);
	@test x[end] == 4
	@test y[1] == 10

	# interp2
	x = [1.,2.,3.,4.,5.,];
	y = [2.,3.,4.,5.,6.,7.,8.];
	z = ones(Float64,(length(y),length(x)));
	z[2,3] = 2.;
	xi = [4.1,3.0,3.0,3.5,5.0,5.01]
	yi = [2.1,3.0,3.5,3.5,2.0,3.30]
	zi = interp2(x,y,z,xi,yi);
	zi_true = [1.0,2.0,1.5,1.25,1.0,NaN]
	for (i,v) in enumerate(zi_true)
		if isnan(v)
			@test isnan(zi[i])
		else
			@test zi[i] ≈ v;
		end
	end

	# scalar input
	@test isnan(interp2(x,y,z,xi,yi)[end])

	# meshgrid input
	xi,yi = meshgrid([1.,2.],[2.,3.])
	zi = interp2(x,y,z,xi,yi);
	@test zi ≈ [1.0 1.0; 1.0 1.0];
end


text:
using Combinatorics

mutable struct SimpleMatrix{I, M, T, J}
    rows::Dict{MonSigPair{I, M}, Polynomial{J, T}}
    tbl::EasyTable{M, J}
end

mutable struct F5matrix{I, M, T, J,
                        SΓ <: SigPolynomialΓ{I, M, T},
                        O <: Base.Order.Ordering}
    sigtail_mat::SimpleMatrix{I, M, T, J}
    sigs_rows::SortedDict{MonSigPair{I, M}, Polynomial{J, T}, O}
    tbl::EasyTable{M, J}
    ctx::SΓ
    max_pos::I
    max_posit_key::I
    tag::Symbol
    trace_sig_tail_tags::Vector{Symbol}
end

function f5matrix(ctx::SigPolynomialΓ{I, M, T},
                  mons::Vector{M},
                  row_sigs::Vector{MonSigPair{I, M}},
                  max_pos = zero(I),
                  max_posit_key = zero(I),
                  tag = :f;
                  trace_sig_tail_tags = Symbol[],
                  enable_lower_pos_rewrite = true,
                  interreduction_step = false) where {I, M, T}

    get_orig_elem = p -> interreduction_step || (!(enable_lower_pos_rewrite) && pos(ctx, p) < max_pos)
    tbl = easytable(mons)
    sigtail_mons = M[]

    sigs_rows = Array{Tuple{MonSigPair{I, M}, Polynomial{ind_type(tbl), T}}}(undef, length(row_sigs))
    sigtail_sigs_rows_unind = Tuple{MonSigPair{I, M}, Polynomial{M, T}}[]
    for (i, sig) in enumerate(row_sigs)
        sig_dat = ctx(sig..., orig_elem = get_orig_elem(sig))
        sigs_rows[i] = (sig, indexpolynomial(tbl, sig_dat[:poly]))
        if !(interreduction_step) && pos(ctx, sig) == max_pos && gettag(ctx, sig) in trace_sig_tail_tags
            sig_tail = project(mul(ctx, sig...), sig_dat[:sigtail])
            append!(sigtail_mons, sig_tail.mo)
            push!(sigtail_sigs_rows_unind, (sig, sig_tail))
        end
    end

    sigtail_tbl = easytable(sort(unique(sigtail_mons),
                                 rev=true,
                                 order=order(ctx.po.mo)))
    sigtail_sigs_rows = Dict([(sig, indexpolynomial(sigtail_tbl, p))
                              for (sig, p) in sigtail_sigs_rows_unind])

    if interreduction_step
        sigs_lms = Dict(zip(row_sigs, [leadingmonomial(sig_row[2]) for sig_row in sigs_rows]))
        row_order = interredorder(ctx, sigs_lms)
    else
        row_order = mpairordering(ctx)
    end

    F5matrix(SimpleMatrix(sigtail_sigs_rows, sigtail_tbl),
             SortedDict(sigs_rows, row_order),
             tbl,
             ctx,
             max_pos,
             max_posit_key,
             tag,
             trace_sig_tail_tags)
end

rows(mat::F5matrix) = values(mat.sigs_rows)
sigs(mat::F5matrix) = keys(mat.sigs_rows)

function check_row_echelon(mat::F5matrix)
    for (p1, p2) in combinations(rows(mat), 2)
        (isempty(p1) || isempty(p2)) && continue
        leadingmonomial(p1) == leadingmonomial(p2) && return false
    end
    return true
end

# for printing matrices
function mat_show(mat)
    mat_vis = zeros(Int, length(mat.sigs_rows), length(mat.tbl.val))
    for (i, (sig, row)) in enumerate(mat.sigs_rows)
        if typeof(mat) <: SimpleMatrix
            for (j, c) in row[2]
                mat_vis[i, j] = Int(c)
            end
        else
            for (j, c) in row
                mat_vis[i, j] = Int(c)
            end
        end
    end
    mat_vis
end

Base.show(io::IO, mat::F5matrix) = Base.show(io, mat_show(mat))

function reduction!(mat::F5matrix{I, M, T, J};
                    interreduction_step = false) where {I, M, T, J, Tbuf}

    ctx = mat.ctx
    n_cols = length(mat.tbl)
    pivots = collect(Base.Iterators.repeated(nullmonsigpair(ctx), n_cols))
    buffer = zeros(buftype(ctx.po.co), n_cols)
    sig_tail_buffer = zeros(buftype(ctx.po.co), length(mat.sigtail_mat.tbl))
    arit_operations_groebner = 0
    arit_operations_module_overhead = 0

    for (sig, row) in mat.sigs_rows
        should_add_sig_tails = pos(ctx, sig) == mat.max_pos && gettag(ctx, sig) in mat.trace_sig_tail_tags
        l = leadingmonomial(row)
        dont_red_cond = interreduction_step ? sig[1] != one(ctx.po.mo) : isnull(pivots[l])
        if dont_red_cond
            pivots[l] = sig
            continue
        end
        buffer!(row, buffer)
        should_add_sig_tails && buffer!(mat.sigtail_mat.rows[sig], sig_tail_buffer)
        for (k, c) in enumerate(buffer)
            (iszero(c) || isnull(pivots[k])) && continue
            arit_operations_groebner += length(mat.sigs_rows[pivots[k]])
            mult = critical_loop!(buffer, mat.sigs_rows[pivots[k]], ctx.po.co)

            # add sig tails
            if should_add_sig_tails && pos(ctx, pivots[k]) == mat.max_pos
                arit_operations_module_overhead += length(mat.sigtail_mat.rows[pivots[k]])
                sub_row!(sig_tail_buffer, mat.sigtail_mat.rows[pivots[k]], mult, ctx.po.co)
            end

        end
        first_nz, new_row = unbuffer!(buffer, ctx.po.co, J)
        if !(iszero(first_nz))
            pivots[first_nz] = sig
        end

        if should_add_sig_tails
            _, new_sig_tail = unbuffer!(sig_tail_buffer, ctx.po.co, J)
            mat.sigtail_mat.rows[sig] = new_sig_tail
        end

        mat.sigs_rows[sig] = new_row
    end
    @debug "matrix num of arit ops in reduction:" arit_operations_groebner
    arit_operations_groebner, arit_operations_module_overhead
end

function buffer!(row::Polynomial{J, T},
                 buffer::Vector{Tbuf}) where {J, M, T, Tbuf}

    [buffer[j] = Tbuf(c) for (j, c) in row]
end

function unbuffer!(buffer::Vector{Tbuf},
                   ctx::NmodLikeΓ{T, Tbuf},
                   ::Type{J}) where {T, Tbuf, J}

    coeffs = T[]
    mons = J[]
    first_nz = 0
    for (j, c) in enumerate(buffer)
        mod_coeff = deflate(ctx, normal(ctx, c))
        iszero(mod_coeff) && continue
        if iszero(first_nz)
            first_nz = J(j)
        end
        push!(coeffs, mod_coeff)
        push!(mons, J(j))
        buffer[j] = zero(Tbuf)
    end
    row = Polynomial(mons, coeffs)
    first_nz, row
end

function sub_row!(buffer::Vector{Tbuf},
                  pivot_row::Polynomial{J, T},
                  mult::T,
                  ctx::NmodLikeΓ{T, Tbuf}) where {J, T, Tbuf}

    for (j, c) in pivot_row
        buffer[j] = submul(ctx, buffer[j], mult, c)
    end
end

function critical_loop!(buffer::Vector{Tbuf},
                        pivot::Polynomial{J, T},
                        ctx::NmodLikeΓ{T, Tbuf}) where {J, T, Tbuf}

    mult1 = deflate(ctx, normal(ctx, buffer[monomial(pivot, 1)]))
    mult2 = inv(ctx, coefficient(pivot, 1))
    mult = mul(ctx, mult1, mult2)
    buffer[monomial(pivot, 1)] = zero(Tbuf)
    try
        for (j, c) in pivot[2:end]
            buffer[j] = submul(ctx, buffer[j], mult, c)
        end
        return mult
    catch BoundsError
        return mult
    end
end

# managing the basis and syzygies after reduction

function new_syz!(ctx::SΓ,
                  sig::MonSigPair{I, M},
                  sig_tail::Polynomial{M, T},
                  pairs::PairSet{I, M, SΓ},
                  H::Syz{I, M}) where {I, M, T, SΓ <: SigPolynomialΓ{I, M, T}}

    new_sig = mul(ctx, sig...)
    ctx(new_sig, zero(eltype(ctx.po)), sig_tail)
    push!(H[new_sig[1]], new_sig[2])
    if gettag(ctx, new_sig) in [:g_gen, :p_gen]
        gettag(ctx, new_sig) == :g_gen ? look_for_tag = :g : look_for_tag = :p
        attached_to = getattkey(ctx, new_sig)
        for (i, Hi) in H
            if ctx.ord_indices[i][:att_key] == attached_to && ctx.ord_indices[i][:tag] == look_for_tag && ctx.ord_indices[i][:position] > pos(ctx, new_sig)
                push!(Hi, new_sig[2])
            end
        end
    end
    new_rewriter!(ctx, pairs, new_sig)
end

function new_basis!(ctx::SΓ,
                    sig::MonSigPair{I, M},
                    poly::Polynomial{M, T},
                    sig_tail::Polynomial{M, T},
                    pairs::PairSet{I, M, SΓ},
                    G::Basis{I, M},
                    H::Syz{I, M};
                    enable_lower_pos_rewrite = true) where {I, M, T, SΓ <: SigPolynomialΓ{I, M, T}}

    m, (pos_key, t) = sig
    new_sig = mul(ctx, sig...)
    # add element to basis if any of the following two conditions hold:
    # reductions of initial generators are added
    add_cond_1 = isone(ctx.po.mo[m]) && isone(ctx.po.mo[t]) && !(t in keys(G[pos_key]))
    # leading term dropped during reduction
    add_cond_2 = lt(ctx.po.mo, leadingmonomial(poly), leadingmonomial(ctx(sig..., orig_elem = true)[:poly]))
    if add_cond_1 || add_cond_2
        ctx(new_sig, poly, sig_tail)
        lm = leadingmonomial(poly)
        new_rewriter!(ctx, pairs, new_sig)
        push!(G[pos_key], (new_sig[2], lm))
        pairs!(ctx, pairs, new_sig, lm, G, H, enable_lower_pos_rewrite = enable_lower_pos_rewrite)
    end
end

function new_elems_f5!(ctx::SΓ,
                       mat::F5matrix{I, M, T},
                       pairs::PairSet{I, M, SΓ},
                       G::Basis{I, M},
                       H::Syz{I, M},
                       info_hashmap::Dict{I, Info},
                       non_zero_cond::Vector{I};
                       enable_lower_pos_rewrite = true) where {I, M, T, SΓ <: SigPolynomialΓ{I, M, T}}

    max_degree = zero(I)
    for (sig, row) in mat.sigs_rows
        @inbounds begin
            if pos(ctx, sig) == mat.max_pos
                gettag(ctx, sig) in mat.trace_sig_tail_tags ? sig_tail = tail(unindexpolynomial(mat.sigtail_mat.tbl, mat.sigtail_mat.rows[sig])) : sig_tail = zero(Polynomial{M, T})
                if isempty(row)
                    new_syz!(ctx, sig, sig_tail, pairs, H)
                else
                    p = unindexpolynomial(mat.tbl, row)
                    d = degree(ctx.po, p)
                    if d > max_degree
                        max_degree = d
                    end
                    if p.mo == [one(ctx.po.mo)]
                        println("constant in basis")
                        if mat.tag == :h
                            filter!(j -> j != mat.max_posit_key, non_zero_cond)
                            while !(isempty(pairs))
                                alpha = first(pairs)[1]
                                if alpha[2][1] == mat.max_posit_key
                                    pop!(pairs)
                                else
                                    break
                                end
                            end
                        end
                    end
                    new_basis!(ctx, sig, p, sig_tail, pairs, G, H, enable_lower_pos_rewrite = enable_lower_pos_rewrite)
                end
            end
        end
    end
    if max_degree > info_hashmap[mat.max_posit_key]["max_deg_reached"]
        info_hashmap[mat.max_posit_key]["max_deg_reached"] = max_degree
    end
    pairs
end

function new_elems_decomp!(ctx::SΓ,
                           mat::F5matrix{I, M, T},
                           pairs::PairSet{I, M, SΓ},
                           G::Basis{I, M},
                           H::Syz{I, M},
                           info_hashmap::Dict{I, Info},
                           non_zero_cond::Vector{I};
                           enable_lower_pos_rewrite = true) where {I, M, T, SΓ <: SigPolynomialΓ{I, M, T}}

    !(mat.tag in [:f, :h]) && return new_elems_f5!(ctx, mat, pairs, G, H, info_hashmap, non_zero_cond, enable_lower_pos_rewrite = enable_lower_pos_rewrite)

    zero_red = filter(sig_row -> isempty(sig_row[2]) && pos(ctx, sig_row[1]) == mat.max_pos && gettag(ctx, sig_row[1]) == mat.tag,
                      mat.sigs_rows)
    isempty(zero_red) && return new_elems_f5!(ctx, mat, pairs, G, H, info_hashmap, non_zero_cond, enable_lower_pos_rewrite = enable_lower_pos_rewrite)

    # insert g's s.t. g*f in I
    pols_to_insert = [unindexpolynomial(mat.sigtail_mat.tbl, mat.sigtail_mat.rows[sig]) for (sig, _) in zero_red]
    g_gen = random_lin_comb(ctx.po, pols_to_insert)

    new_gen!(ctx, info_hashmap, mat.max_pos,
             mat.max_posit_key, mat.tag == :f ? :g_gen : :p_gen, g_gen)
    for (j, (sig, _)) in enumerate(zero_red)
        ctx(mul(ctx, sig...), zero(eltype(ctx.po)), tail(pols_to_insert[j]))
        if j < length(zero_red)
            new_gen!(ctx, info_hashmap, pos(ctx, sig), mat.max_posit_key, mat.tag == :f ? :g : :p, pols_to_insert[j])
        end
    end

    non_triv_syz = [g[2] for g in G[mat.max_posit_key]]
    filter!(pos_elems -> ctx.ord_indices[pos_elems[1]][:position] < mat.max_pos, G)
    pairs = pairset(ctx)
    for posit_key in keys(ctx.ord_indices)
        sg = (posit_key, one(ctx.po.mo))
        if pos(ctx, sg) >= mat.max_pos
            pair!(ctx, pairs, sg)
            G[posit_key] = Tuple{M, M}[]
            if !(posit_key in keys(H))
                H[posit_key] = copy(non_triv_syz)
            end
        end
    end
    pairs
end


text:
module Tuneables

include("types.jl")

end # module

max_stars_repo_name: jonniedie/Tuneables

max_issues_repo_name: jonniedie/Tuneables

max_forks_repo_name: jonniedie/Tuneables


text:
using ConvOptDL
using Test
using StatsBase

@testset "data_loader.jl" begin
    dloader = FewShotDataLoader("./test_data.jls")
    batch = sample(dloader, 2, support_n_ways=3, support_k_shots=5)
    @test size(batch) == 2
    @test size(batch.support_samples) == (32, 32, 3, 15, 2)
    @test size(batch.support_labels) == (15, 2)
    try
        _ = sample(dloader, 5, support_n_ways=3, support_k_shots=5, query_k_shots=2)
    catch err
        @test err.msg == false
    end
end