text:
using DimensionalData, Test, BenchmarkTools, Dates, Statistics
using DimensionalData.LookupArrays, DimensionalData.Dimensions

@testset "dims2indices" begin
    a = [1 2 3; 4 5 6]
    da = DimArray(a, (X(143.0:2:145.0), Y(-38.0:-36.0)))
    dimz = dims(da)

    @test Dimensions._dims2indices(dimz[1], Y) == Colon()
    @test dims2indices(dimz, (Y(),)) == (Colon(), Colon())
    @test (@ballocated dims2indices($dimz, (Y(),))) == 0
    @test dims2indices(dimz, (Y(1),)) == (Colon(), 1)
    @test (@ballocated dims2indices($dimz, (Y(1),))) == 0
    @test dims2indices(dimz, (Ti(4), X(2))) == (2, Colon())
    @test dims2indices(dimz, (Y(2), X(3:7))) == (3:7, 2)
    @test (@ballocated dims2indices($dimz, (Y(2), X(3:7)))) == 0
    @test dims2indices(dimz, (X(2), Y([1, 3, 4]))) == (2, [1, 3, 4])
    @test dims2indices(da, (X(2), Y([1, 3, 4]))) == (2, [1, 3, 4])
    v = [1, 3, 4]
    @test (@ballocated dims2indices($da, (X(2), Y($v)))) == 0
    @test_throws ArgumentError dims2indices(nothing, (Y(2), X(3:7)))
    @test dims2indices(dimz, CartesianIndex(1, 1)) == (CartesianIndex(1, 1),)
    @test dims2indices(dimz[1], 1) == 1
    @test dims2indices(dimz[1], X(2)) == 2
end

@testset "lookup" begin
    l = Sampled(2.0:2.0:10, ForwardOrdered(), Regular(2.0), Points(), nothing)
    @test l[:] == l
    @test l[1:5] == l
    @test l[1:5] isa typeof(l)
    @test l[Near(2.1)] == 2.0
    # TODO properly handle index mashing arrays: here Regular should become Irregular
    # @test d[[1, 3, 4]] == X(Sampled([2.0, 6.0, 8.0], ForwardOrdered(), Regular(2.0), Points(), nothing))
    # @test d[[true, false, false, false, true]] == X(Sampled([2.0, 10.0], ForwardOrdered(), Regular(2.0), Points(), nothing))
    @test l[2] === 4.0
    @test l[CartesianIndex((4,))] == 8.0
    @test l[2:2:4] == Sampled(4.0:4.0:8.0, ForwardOrdered(), Regular(2.0), Points(), nothing)
    l = NoLookup(1:100)
    l = NoLookup(1:100)
    @test l[100] == 100
    @test l[1:5] isa typeof(l)
    @test l[CartesianIndex((1,))] == 1
    @test view(l, 50) === view(1:100, 50)
    @test Base.dotview(l, 50) === 50
end

@testset "dimension" begin
    d = X(Sampled(2.0:2.0:10, ForwardOrdered(), Regular(2.0), Points(), nothing))
    @test d[:] == d
    @test d[1:5] == d
    @test d[1:5] isa typeof(d)
    # TODO properly handle index mashing arrays: here Regular should become Irregular
    # @test d[[1, 3, 4]] == X(Sampled([2.0, 6.0, 8.0], ForwardOrdered(), Regular(2.0), Points(), nothing))
    # @test d[[true, false, false, false, true]] == X(Sampled([2.0, 10.0], ForwardOrdered(), Regular(2.0), Points(), nothing))
    @test d[2] === 4.0
    @test d[CartesianIndex((4,))] == 8.0
    @test d[2:2:4] == X(Sampled(4.0:4.0:8.0, ForwardOrdered(), Regular(2.0), Points(), nothing))
    d = Y(NoLookup(1:100))
    @test d[100] == 100
    @test d[1:5] isa typeof(d)
    @test d[CartesianIndex((1,))] == 1
    @test view(d, 50) === view(1:100, 50)
    @test Base.dotview(d, 50) === 50
end

@testset "array" begin
    a = [1 2; 3 4]
    xmeta = Metadata(:meta => "X")
    ymeta = Metadata(:meta => "Y")
    ameta = Metadata(:meta => "da")
    dimz = (X(Sampled(143.0:2:145.0; order=ForwardOrdered(), metadata=xmeta)),
            Y(Sampled(-38.0:2:-36.0; order=ForwardOrdered(), metadata=ymeta)))
    refdimz = (Ti(1:1),)
    da = @test_nowarn DimArray(a, dimz; refdims=refdimz, name=:test, metadata=ameta)

    @testset "getindex for single integers returns values" begin
        @test da[X(1), Y(2)] == 2
        @test da[X(2), Y(2)] == 4
        @test da[1, 2] == 2
        @test da[2] == 3
        @inferred getindex(da, X(2), Y(2))
    end

    @testset "LinearIndex getindex returns an Array, except Vector" begin
        @test da[1:2] isa Array
        @test x = da[1, :][1:2] isa DimArray
    end

    @testset "mixed CartesianIndex indexing works" begin

        da3 = cat(da, 10da; dims=Z)
        @test da3[1, CartesianIndex(1, 2)] == 10
        @test view(da3, 1:2, CartesianIndex(1, 2)) == [10, 30]
    end

    @testset "getindex returns DimensionArray slices with the right dimensions" begin
        a = da[X(1:2), Y(1)]
        @test a == [1, 3]
        @test typeof(a) <: DimArray{Int,1}
        @test dims(a) == (X(Sampled(143.0:2.0:145.0, ForwardOrdered(), Regular(2.0), Points(), xmeta)),)
        @test refdims(a) ==
            (Ti(1:1), Y(Sampled(-38.0:2.0:-38.0, ForwardOrdered(), Regular(2.0), Points(), ymeta)),)
        @test name(a) == :test
        @test metadata(a) === ameta
        @test metadata(a, X) === xmeta
        @test bounds(a) === ((143.0, 145.0),)
        @test bounds(a, X) === (143.0, 145.0)
        @test locus(da, X) == Center()

        a = da[X(1), Y(1:2)]
        @test a == [1, 2]
        @test typeof(a) <: DimArray{Int,1}
        @test typeof(parent(a)) <: Array{Int,1}
        @test dims(a) == (Y(Sampled(-38.0:2.0:-36.0, ForwardOrdered(), Regular(2.0), Points(), ymeta)),)
        @test refdims(a) ==
            (Ti(1:1), X(Sampled(143.0:2.0:143.0, ForwardOrdered(), Regular(2.0), Points(), xmeta)),)
        @test name(a) == :test
        @test metadata(a) == ameta
        @test bounds(a) == ((-38.0, -36.0),)
        @test bounds(a, Y()) == (-38.0, -36.0)

        a = da[X(:), Y(:)]
        @test a == [1 2; 3 4]
        @test typeof(a) <: DimArray{Int,2}
        @test typeof(parent(a)) <: Array{Int,2}
        @test typeof(dims(a)) <: Tuple{<:X,<:Y}
        @test dims(a) == (X(Sampled(143.0:2.0:145.0,
                            ForwardOrdered(), Regular(2.0), Points(), xmeta)),
                          Y(Sampled(-38.0:2.0:-36.0,
                            ForwardOrdered(), Regular(2.0), Points(), ymeta)))
        @test refdims(a) == (Ti(1:1),)
        @test name(a) == :test
        @test bounds(a) == ((143.0, 145.0), (-38.0, -36.0))
        @test bounds(a, X) == (143.0, 145.0)

        # Indexing with array works
        a = da[X([2, 1]), Y([2, 1])]
        @test a == [4 3; 2 1]
    end

    @testset "selectors work" begin
        @test da[At(143), -38.0..36.0] == [1, 2]
        @test da[144.0..146.0, Near(-37.1)] == [3]
        @test da[X=At(143), Y=-38.0..36.0] == [1, 2]
        @test da[X=144.0..146.0, Y=Near(-37.1)] == [3]
    end

    @testset "view DimensionArray containing views" begin
        v = view(da, Y(1), X(1))
        @test v[] == 1
        @test typeof(v) <: DimArray{Int,0}
        @test typeof(parent(v)) <:SubArray{Int,0}
        @test typeof(dims(v)) == Tuple{}
        @test dims(v) == ()
        @test refdims(v) ==
            (Ti(1:1), X(Sampled(143.0:2.0:143, ForwardOrdered(), Regular(2.0), Points(), xmeta)),
             Y(Sampled(-38.0:2.0:-38.0, ForwardOrdered(), Regular(2.0), Points(), ymeta)))
        @test name(v) == :test
        @test metadata(v) == ameta
        @test bounds(v) == ()

        v = view(da, Y(1), X(1:2))
        @test v == [1, 3]
        @test typeof(v) <: DimArray{Int,1}
        @test typeof(parent(v)) <: SubArray{Int,1}
        @test typeof(dims(v)) <: Tuple{<:X}
        @test dims(v) == (X(Sampled(143.0:2.0:145.0, ForwardOrdered(), Regular(2.0), Points(), xmeta)),)
        @test refdims(v) ==
            (Ti(1:1), Y(Sampled(-38.0:-38.0, ForwardOrdered(), Regular(2.0), Points(), ymeta)),)
        @test name(v) == :test
        @test metadata(v) == ameta
        @test bounds(v) == ((143.0, 145.0),)

        # Test that dims are actually views using a vector
        da_vec = rebuild(da; dims=map(d -> rebuild(d, rebuild(lookup(d); data=Array(d))), dims(da)))
        v = view(da_vec, Y(1:2), X(1:1))
        @test v == [1 2]
        @test typeof(v) <: DimArray{Int,2}
        @test typeof(parent(v)) <: SubArray{Int,2}
        @test typeof(dims(v)) <: Tuple{<:X,<:Y}
        testdims = (X(Sampled(view([143.0, 143.0], 1:1), ForwardOrdered(), Regular(2.0), Points(), xmeta)),
             Y(Sampled(view([-38.0, -36.0], 1:2), ForwardOrdered(), Regular(2.0), Points(), ymeta)))
        @test typeof(dims(v)) == typeof(testdims)
        @test dims(v) == testdims
        @test bounds(v) == ((143.0, 143.0), (-38.0, -36.0))

        v = view(da, Y(Base.OneTo(2)), X(1))
        @test v == [1, 2]
        @test typeof(parent(v)) <: SubArray{Int,1}
        @test typeof(dims(v)) <: Tuple{<:Y}
        @test dims(v) == (Y(Sampled(-38.0:2.0:-36.0, ForwardOrdered(), Regular(2.0), Points(), ymeta)),)
        @test refdims(v) ==
            (Ti(1:1), X(Sampled(143.0:2.0:143.0, ForwardOrdered(), Regular(2.0), Points(), xmeta)),)
        @test bounds(v) == ((-38.0, -36.0),)

        @test view(da, 1) == fill(1)
    end

    @testset "setindex!" begin
        da_set = copy(da)
        da_set[X(2), Y(2)] = 77
        @test da_set == [1 2; 3 77]
        da_set[X(1:2), Y(1)] .= 99
        @test da_set == [99 2; 99 77]
        da_set[1] = 55
        @test da_set == [55 2; 99 77]
        da_set[1, 1, 1, 1] = 66
        @test da_set == [66 2; 99 77]
        @test_throws BoundsError da_set[2, 2, 2, 2] = 66
        da_set[:, :] = [1 2; 3 4]
        @test da_set == [1 2; 3 4]
    end

    @testset "logical indexing" begin
        A = DimArray(zeros(40, 50), (X, Y));
        I = rand(40, 50) .< 0.5
        @test all(A[I] .== 0.0)
        A[I] .= 3
        @test all(A[I] .== 3.0)
        @test all(view(A, I .== 3.0))
    end

    @testset "Empty getindedex/view/setindex throws a BoundsError" begin
        @test_throws BoundsError da[]
        @test_throws BoundsError view(da)
        @test_throws BoundsError da[] = 1
    end

    @testset "Cartesian indices work as usual" begin
        @test da[CartesianIndex(2, 2)] == 4
        @test view(da, CartesianIndex(2, 2)) == view(parent(da), 2, 2)
        da_set = deepcopy(da)
        da_set[CartesianIndex(2, 2)] = 5
        @test da_set[2, 2] == 5
    end

    a2 = [1 2 3 4
          3 4 5 6
          4 5 6 7]
    b2 = [4 4 4 4
          4 4 4 4
          4 4 4 4]

    @testset "indexing into NoLookup dims is just regular indexing" begin
        ida = DimArray(a2, (X(), Y()))
        ida[Y(3:4), X(2:3)] = [5 6; 6 7]
    end

    dimz2 = (Dim{:row}(10:10:30), Dim{:column}(-20:10:10))
    da2 = DimArray(a2, dimz2; refdims=refdimz, name=:test2)

    @testset "lookup step is updated when indexed with a range" begin
        @test step.(lookup(da2)) == (10.0, 10.0)
        @test step.(lookup(da2[1:3, 1:4])) == (10.0, 10.0)
        @test step.(lookup(da2[1:2:3, 1:3:4])) == (20.0, 30.0)
        @test step.(lookup(da2[column=1:2:4, row=1:3:3])) == (30.0, 20.0)
    end

    @testset "Symbol dimension names also work for indexing" begin
        @test da2[Dim{:row}(2)] == [3, 4, 5, 6]
        @test da2[Dim{:column}(4)] == [4, 6, 7]
        @test da2[column=4] == [4, 6, 7]
        @test da2[Dim{:column}(1), Dim{:row}(3)] == 4
        @test da2[column=1, Dim{:row}(3)] == 4
        @test da2[Dim{:column}(1), row=3] == 4
        @test da2[column=1, row=3] == 4
        @test view(da2, column=1, row=3) == fill(4)
        @test view(da2, column=1, Dim{:row}(1)) == fill(1)
        da2_set = deepcopy(da2)
        da2_set[Dim{:column}(1), Dim{:row}(1)] = 99
        @test da2_set[1, 1] == 99
        da2_set[column=2, row=2] = 88
        @test da2_set[2, 2] == 88
        da2_set[Dim{:row}(3), column=3] = 77
        @test da2_set[3, 3] == 77

        # We can also construct without using `Dim{X}`
        @test dims(DimArray(a2, (:a, :b))) == dims(DimArray(a2, (Dim{:a}, Dim{:b})))
    end

    @testset "Type inference holds indexing with Symbol dimension names" begin
        da2_set = deepcopy(da2)
        # Inrerence
        @inferred getindex(da2, column=1, row=3)
        @inferred view(da2, column=1, row=3)
        @inferred setindex!(da2_set, 77, Dim{:row}(1), column=2)

        # With a large type
        if VERSION >= v"1.5"
            da4 = DimArray(zeros(1, 2, 3, 4, 5, 6, 7, 8), (:a, :b, :c, :d, :d, :f, :g, :h))
            @inferred getindex(da2, a=1, b=2, c=3, d=4, e=5)
            # Type inference breaks with 6 arguments.
            # @inferred getindex(da2, a=1, b=2, c=3, d=4, e=5, f=6)
            # @code_warntype getindex(da2, a=1, b=2, c=3, d=4, e=5, f=6)
        end
    end
end

@testset "stack" begin
    A = [1.0 2.0 3.0;
         4.0 5.0 6.0]
    dimz = (X([:a, :b]), Y(10.0:10.0:30.0))
    da1 = DimArray(A, dimz; name=:one)
    da2 = DimArray(Float32.(2A), dimz; name=:two)
    da3 = DimArray(Int.(3A), dimz; name=:three)

    s = DimStack((da1, da2, da3))

    @testset "getindex" begin
        @test s[1, 1] === (one=1.0, two=2.0f0, three=3)
        @test s[X(2), Y(3)] === (one=6.0, two=12.0f0, three=18)
        @test s[X=At(:b), Y=At(10.0)] === (one=4.0, two=8.0f0, three=12)
        slicedds = s[At(:a), :]
        @test slicedds[:one] == [1.0, 2.0, 3.0]
        @test slicedds.data == (one=[1.0, 2.0, 3.0], two=[2.0f0, 4.0f0, 6.0f0], three=[3, 6, 9])
        @testset "linear indices" begin
            linear2d = s[1:2]
            @test linear2d isa NamedTuple
            @test linear2d == (one=[1.0, 4.0], two=[2.0f0, 8.0f0], three=[3, 12])
            linear1d = s[Y(1)][1:2]
            @test linear1d isa DimStack
            @test linear1d.data == (one=[1.0, 4.0], two=[2.0f0, 8.0f0], three=[3, 12])
        end
    end

    @testset "getindex Tuple" begin
        st1 = s[(:three, :one)]
        @test keys(st1) === (:three, :one)
        @test values(st1) == (da3, da1)
    end

    @testset "view" begin
        sv = view(s, 1, 1)
        @test sv.data == (one=fill(1.0), two=fill(2.0f0), three=fill(3))
        @test dims(sv) == ()
        sv = view(s, X(1:2), Y(3:3))
        @test sv.data == (one=[3.0 6.0]', two=[6.0f0 12.0f0]', three=[9 18]')
        slicedds = view(s, X=At(:a), Y=:)
        @test slicedds[:one] == [1.0, 2.0, 3.0]
        @test slicedds.data == (one=[1.0, 2.0, 3.0], two=[2.0f0, 4.0f0, 6.0f0], three=[3, 6, 9])
        @testset "linear indices" begin
            linear2d = view(s, 1)
            @test linear2d isa NamedTuple
            @test linear2d == (one=fill(1.0), two=fill(2.0f0), three=fill(3))
            linear1d = view(s[X(1)], 1)
            @test linear1d == (one=fill(1.0), two=fill(2.0f0), three=fill(3))
            # Its not clear if this should work or not
            @test_broken linear1d isa DimStack
        end
    end

    @testset "setindex!" begin
        s_set = deepcopy(s)
        s_set[1, 1] = (one=9, two=10, three=11)
        @test s_set[1, 1] === (one=9.0, two=10.0f0, three=11)
        s_set[X=At(:b), Y=At(10.0)] = (one=7, two=11, three=13)
        @test s_set[2, 1] === (one=7.0, two=11.0f0, three=13)
    end

    @testset "Empty getindedex/view/setindex throws a BoundsError" begin
        @test_throws BoundsError s[]
        @test_throws BoundsError view(s)
        @test_throws BoundsError s[] = 1
    end

    @testset "Cartesian indices work as usual" begin
        @test s[CartesianIndex(2, 2)] == (one=5.0, two=10.0, three=15.0)
        @test view(s, CartesianIndex(2, 2)) == map(d -> view(d, 2, 2), s.data)
        s_set = deepcopy(s)
        s_set[CartesianIndex(2, 2)] = (one=5, two=6, three=7)
        @test s_set[2, 2] === (one=5.0, two=6.0f0, three=7)
        s_set[CartesianIndex(2, 2)] = (9, 10, 11)
        @test s_set[2, 2] === (one=9.0, two=10.0f0, three=11)
        @test_throws ArgumentError s_set[CartesianIndex(2, 2)] = (seven=5, two=6, three=7)
    end
end

@testset "indexing irregular with bounds slice" begin
    total_time = 5
    t1 = Ti([mod1(i, 12) for i in 1:total_time]; sampling=Intervals()) # centurial averages with seasonal cycle
    t2 = Ti([Date(-26000 + ((i-1)÷12)*100, mod1(i, 12), 1) for i in 1:total_time]; sampling=Intervals()) # centurial averages with seasonal cycle
    t3 = Ti([DateTime(-26000 + ((i-1)÷12)*100, mod1(i, 12), 1) for i in 1:total_time]; sampling=Intervals()) # centurial averages with seasonal cycle
    t4 = Ti([DateTime(-26000 + ((i-1)÷12)*100, mod1(i, 12), 1) for i in 1:total_time]; sampling=Points()) # centurial averages with seasonal cycle

    p = reshape([1, 2, 3, 4, 5], 1, 1, 5)
    A1 = DimArray(p, (X, Y, t1))
    A2 = DimArray(p, (X, Y, t2))
    A3 = DimArray(p, (X, Y, t3))

    @test view(A1, Ti(5)) == [5;;]
    @test view(A2, Ti(5)) == [5;;]
    @test view(A3, Ti(5)) == [5;;]
end


text:
module Sciplotlib

using Reexport
@reexport using Colors
@reexport using PyPlot: PyPlot, plt, matplotlib as mpl
#    `plt` is set to mpl.pyplot in PyPlot's __init__. `PyPlot as plt` doesn't work.

using PyCall
using Unitful
using Unitful: Units
using ColorVectorSpace
using Printf
using PartialFunctions: $

include("plot.jl")
export plot

include("set.jl")
export set, legend, hylabel

include("colors.jl")
export mix, lighten, darken, deemph, toRGBAtuple,
       black, white, lightgrey, mplcolors,
       C0, C1, C2, C3, C4, C5, C6, C7, C9, C10

include("typecompat.jl")
export as_mpl_type

include("ticks.jl")
include("util.jl")

include("style.jl")
export sciplotlib_style, set_mpl_style!

include("init.jl")
export rcParams, rcParams_original

end

ext: jl


text:
@testset "$(rpad("HDF5 Utils",80))" begin

    using HDF5

    ntime = 10
    nsave = 1
    dim   = 1

    sol = SolutionODE{Float64}(dim, ntime, nsave)

    for i in 1:ntime
        sol[1,i] = i
    end

    tmp = mktempdir()
    file = tmp * "/test.h5"
    writeSolutionToHDF5(sol, file)
    h5 = h5open(file, "r")
    @test h5["x"][:,:] == sol.x
    close(h5)
    rm(tmp, recursive=true)

    tmp = mktempdir()
    h5 = createHDF5(sol, tmp * "/test.h5", 2ntime+1)
    writeSolutionToHDF5(sol, h5)
    writeSolutionToHDF5(sol, h5, ntime)
    close(h5)
    rm(tmp, recursive=true)

end


text:
"""
Problema de asignación. (Vanderbei LP, pg. 227). Dado un conjunto S de m personas, un conjunto D de m tareas, y para cada iS,jD un costo cij de asignar la tarea j a la persona i, el problema de asignación consiste en asignar a cada persona una y solo una tarea, de manera tal que se minimize el costo total. Formule este problema como un problema de programación lineal con enteros. Escriba una función para resolverlo, que tome como entrada una matriz de costos cij arbitraria.

"""

export solve_with_lineal_programming;
function solve_with_lineal_programming(m::Matrix{T}) where {T<:Number}

	#Checking matrix
	checkmatrix(m);
	ml = size(m,1);

	#Create model
	local model = Model(with_optimizer(GLPK.Optimizer));

	#Variables
	@variable(model, x[1:ml,1:ml], Bin);

	#Objective
	@objective(model, Min, sum(m[i,j] * x[i,j] for i in 1:ml, j in 1:ml));

	#Constrains
	@constraint(model, [i=1:ml], sum(x[i,j] for j in 1:ml) <= 1);
	@constraint(model, [j=1:ml], sum(x[i,j] for i in 1:ml) >= 1);

	#Optimize
	optimize!(model);

	return pick_results(value.(x),m);
end

function pick_results(bm,m)
	r = Vector{eltype(m)}();
	zer0 = zero(eltype(m));
	for i in 1:length(bm)
		bm[i] != zer0 && push!(r,m[i]);
	end
	return push!(Vector{typeof(r)}(),r);
end


text:
using Distances
dist = Euclidean()

tic()
tic()
# file lines iterator object
input = eachline(STDIN)

state = start(input)
# read data size N
(n, state) = next(input, state)
n = parse(Int16, n)

# 2xN node's coordinates matrix
# Data[:, k] = [x_k, y_k]
data = Array{Float64}(2, n)
# NxN node's distances matric
d = Array{Float64}(n, n)

# read file data to Data matrix
while !done(input, state)
    (v, state) = next(input, state)
    (i_v, x_v, y_v) = split(v)
    ii = parse(Int32, i_v)
    xx = parse(Float64, x_v)
    yy = parse(Float64, y_v)
    data[:,ii] = [xx, yy]
end
# calculate nodes distances as NxN matrix
pairwise!(d, dist, data)

# show(d)

d_orig = copy(d)
dim1 = size(d,1)
dim12 = size(d)

for i in 1:dim1
    d[i,i] = 10.0e6
end

d1 = copy(d)
tour = zeros(dim12)
cost = 0
min_dist = Array{Float32}(dim1)
short_path = Array{Int32}(dim1)
best_nbr_cost = 0
best_nbr = Array{Float64}(dim1)

k = 1
for i in 1:dim1-1
    min_dist[i] = min(d1[k,:]...)
    short_path[i] = findfirst(d1[k,:].==min(d1[k,:]...))
    cost = cost + min_dist[i]
    k = short_path[i]
    d1[k,1] = 10.0e6
    for visited_node in 1:i
        d1[k, short_path[visited_node]] = 10.0e6
    end
end

tour[1, short_path[1]] = 1
for i in 2:dim1-1
    tour[short_path[i-1], short_path[i]] = 1
end

short_path[dim1] = 1
tour[k, short_path[dim1]] = 1
cost = cost + d[k,1]

crnt_tour = copy(short_path)
best_tour = copy(short_path)
best_obj = cost
crnt_tour_cost = cost

print("initial cost: ")
println(best_obj)
println(string("initialization time: ", toq(), "s"))

nbr_cost = Array{Float64}(dim12)

tabu_tenure = zeros(dim12)
max_tabu_tenure = round(sqrt(dim1))

penalty = zeros(1, Integer((dim1-1)*(dim1-2)/2))
frequency = zeros(dim12)
frequency[1,:] = 100_000
frequency[:,1] = 100_000

for i in 1:dim1
    frequency[i,i] = 100_000
end

iter_snc_last_imprv = 0

best_nbr = copy(crnt_tour)
best_i = 0
best_j = 0

for iter in 1:Integer(floor((10_000_000*sqrt(n))/n^2))
    # @printf("Iteration %d\n", iter)
    if crnt_tour_cost < 0
        @printf("ERROR! crnt_tour_cost < 0\n")
        println(crnt_tour)
        return
    end
    nbr_cost = fill(Inf, dim12)
    for i in 1:dim1-2
        for j in i+1:dim1-1
            if i == 1
                if j-1 == 1
                    nbr_cost[crnt_tour[i], crnt_tour[j]] = crnt_tour_cost - d[1,crnt_tour[i]] + d[1,crnt_tour[j]] - d[crnt_tour[j],crnt_tour[j+1]] + d[crnt_tour[i],crnt_tour[j+1]]
                    best_i = i
                    best_j = j
                    best_nbr_cost = nbr_cost[crnt_tour[i],crnt_tour[j]]
                    tabu_node1 = crnt_tour[i]
                    tabu_node2 = crnt_tour[j]
                else
                    nbr_cost[crnt_tour[i], crnt_tour[j]] = crnt_tour_cost - d[1,crnt_tour[i]] + d[1,crnt_tour[j]] - d[crnt_tour[j],crnt_tour[j+1]] + d[crnt_tour[i],crnt_tour[j+1]] - d[crnt_tour[i],crnt_tour[i+1]] + d[crnt_tour[j],crnt_tour[i+1]] - d[crnt_tour[j-1],crnt_tour[j]] + d[crnt_tour[j-1],crnt_tour[i]]
                end
            else
                if j-1 == 1
                    nbr_cost[crnt_tour[i],crnt_tour[j]] = crnt_tour_cost - d[crnt_tour[i-1],crnt_tour[i]] + d[crnt_tour[i-1],crnt_tour[j]] - d[crnt_tour[j],crnt_tour[j+1]] + d[crnt_tour[i],crnt_tour[j+1]]
                else
                    nbr_cost[crnt_tour[i],crnt_tour[j]] = crnt_tour_cost - d[crnt_tour[i-1],crnt_tour[i]] + d[crnt_tour[i-1],crnt_tour[j]] - d[crnt_tour[j],crnt_tour[j+1]] + d[crnt_tour[i],crnt_tour[j+1]] -
                    d[crnt_tour[i],crnt_tour[i+1]] + d[crnt_tour[j],crnt_tour[i+1]] - d[crnt_tour[j-1],crnt_tour[j]] + d[crnt_tour[j-1],crnt_tour[i]]
                end
            end

            if nbr_cost[crnt_tour[i],crnt_tour[j]] < best_nbr_cost
                best_nbr_cost = nbr_cost[crnt_tour[i],crnt_tour[j]]
                best_i = i
                best_j = j
                tabu_node1 = crnt_tour[i]
                tabu_node2 = crnt_tour[j]
                if tabu_node1 > tabu_node2
                    tabu_node1, tabu_node2 = tabu_node2, tabu_node1
                end
            end
        end
    end

    best_nbr[best_i] = crnt_tour[best_j]
    best_nbr[best_j] = crnt_tour[best_i]

    while tabu_tenure[tabu_node1,tabu_node2] > 0
        if best_nbr_cost < best_obj
            break
        else
            nbr_cost[tabu_node1,tabu_node2] = nbr_cost[tabu_node1,tabu_node2] * 1000
            best_nbr_cost_col = minimum(nbr_cost,2)
            best_nbr_cost = minimum(best_nbr_cost_col,1)[1]
            (R,C) = collect(zip(ind2sub(size(nbr_cost),find(nbr_cost.==best_nbr_cost))...))[1]
            tabu_node1 = R
            tabu_node2 = C
        end
    end

    if best_nbr_cost > crnt_tour_cost
        min_d_col = minimum(d,2)
        penal_nbr_cost = nbr_cost + minimum(min_d_col,1)[1]*frequency
        penal_best_nbr_cost_col = minimum(penal_nbr_cost,2)
        penal_best_nbr_cost = minimum(penal_best_nbr_cost_col,1)[1]
        (Rp,Cp) = collect(zip(ind2sub(size(penal_nbr_cost),find(penal_nbr_cost.==penal_best_nbr_cost))...))[1]
        tabu_node1 = Rp
        tabu_node2 = Cp
        best_nbr_cost = nbr_cost[tabu_node1,tabu_node2]
    end

    for row in 1:dim1-1
        for col in row+1:dim1
            if tabu_tenure[row,col] > 0
                tabu_tenure[row,col] = tabu_tenure[row,col] - 1
                tabu_tenure[col,row] = tabu_tenure[row,col]
            end
        end
    end

    tabu_tenure[tabu_node1,tabu_node2] = max_tabu_tenure
    tabu_tenure[tabu_node2,tabu_node1] = tabu_tenure[tabu_node1,tabu_node2]

    frequency[tabu_node1,tabu_node2] = frequency[tabu_node1,tabu_node2] + 1

    crnt_tour = copy(best_nbr)
    crnt_tour_cost = best_nbr_cost

    if crnt_tour_cost < best_obj
        best_obj = crnt_tour_cost
        # @printf("Iteration %d\n", iter)
        # @printf("best obj = %d\n", best_obj)
        best_tour = copy(crnt_tour)
        iter_snc_last_imprv = 0
    else
        iter_snc_last_imprv += 1
        if iter_snc_last_imprv >= 30*log2(n)
            # @printf("Iterations since last improvement exceeded...\n")
            min_freq_col = minimum(frequency,2)
            min_freq = minimum(min_freq_col,1)[1]
            (R,C) = collect(zip(ind2sub(size(frequency),find(frequency.==minimum(minimum(frequency,2),1)[1]))...))[1]
            freq_indx1 = R
            freq_indx2 = C
            indx_in_crnt_tour1 = findfirst(crnt_tour.==R)
            indx_in_crnt_tour2 = findfirst(crnt_tour.==C)
            # Diversify using a move that has the lowest frequency
            crnt_tour[indx_in_crnt_tour1], crnt_tour[indx_in_crnt_tour2] = crnt_tour[indx_in_crnt_tour2], crnt_tour[indx_in_crnt_tour1]
            tabu_tenure = zeros(dim12)
            frequency = zeros(dim12)
            frequency[1,:] = 100_000
            frequency[:,1] = 100_000
            for i in 1:dim1
                frequency[i,i] = 100_000
            end
            tabu_tenure[R,C] = max_tabu_tenure
            tabu_tenure[C,R] = max_tabu_tenure
            frequency[R,C] = frequency[R,C] + 1
            frequency[C,R] = frequency[R,C]
            # Re-calculare crnt tour cost
            crnt_tour_cost = d[1,crnt_tour[1]]
            for i in 1:dim1-1
                crnt_tour_cost = crnt_tour_cost + d[crnt_tour[i],crnt_tour[i+1]]
            end
            iter_snc_last_imprv = 0
            if crnt_tour_cost < best_obj
                best_obj = crnt_tour_cost
                best_tour = copy(crnt_tour)
            end
        end
    end
    # @printf("current tour cost = %d\t", crnt_tour_cost)
    # @printf("best obj = %d\n", best_obj)
end
# println(STDOUT, best_obj)
# for n in best_tour
#     print(STDERR, n)
#     print(STDERR, " ")
# end
# println(STDERR)

println(STDOUT, string("best cost: ", best_obj))
println(string("total time: ", toq(), "s"))

max_stars_repo_name: AfroMetal/algorithms-tsp-tabusearch

max_issues_repo_name: AfroMetal/algorithms-tsp-tabusearch

max_forks_repo_name: AfroMetal/algorithms-tsp-tabusearch


text:
"""
    CurvedMesh{NDIMS, RealT<:Real} <: AbstractMesh{NDIMS}

A structured curved mesh.

Different numbers of cells per dimension are possible and arbitrary functions
can be used as domain faces.

!!! warning "Experimental code"
    This mesh type is experimental and can change any time.
"""
mutable struct CurvedMesh{NDIMS, RealT<:Real} <: AbstractMesh{NDIMS}
  cells_per_dimension::NTuple{NDIMS, Int}
  mapping::Any # Not relevant for performance
  mapping_as_string::String
  periodicity::NTuple{NDIMS, Bool}
  current_filename::String
  unsaved_changes::Bool
end

"""
    CurvedMesh(cells_per_dimension, mapping, RealT; unsaved_changes=true, mapping_as_string=mapping2string(mapping, length(cells_per_dimension)))

Create a CurvedMesh of the given size and shape that uses `RealT` as coordinate type.

# Arguments
- `cells_per_dimension::NTupleE{NDIMS, Int}`: the number of cells in each dimension.
- `mapping`: a function of `NDIMS` variables to describe the mapping, which transforms
             the reference mesh to the physical domain.
- `RealT::Type`: the type that should be used for coordinates.
- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}`
                 deciding for each dimension if the boundaries in this dimension are periodic.
- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file.
- `mapping_as_string::String`: the code that defines the `mapping`.
                               If `CodeTracking` can't find the function definition, it can be passed directly here.
"""
function CurvedMesh(cells_per_dimension, mapping; RealT=Float64, periodicity=true, unsaved_changes=true,
                    mapping_as_string=mapping2string(mapping, length(cells_per_dimension)))
  NDIMS = length(cells_per_dimension)

  # Convert periodicity to a Tuple of a Bool for every dimension
  if all(periodicity)
    # Also catches case where periodicity = true
    periodicity = ntuple(_->true, NDIMS)
  elseif !any(periodicity)
    # Also catches case where periodicity = false
    periodicity = ntuple(_->false, NDIMS)
  else
    # Default case if periodicity is an iterable
    periodicity = Tuple(periodicity)
  end

  return CurvedMesh{NDIMS, RealT}(Tuple(cells_per_dimension), mapping, mapping_as_string, periodicity, "", unsaved_changes)
end

"""
    CurvedMesh(cells_per_dimension, faces, RealT; unsaved_changes=true, faces_as_string=faces2string(faces))

Create a CurvedMesh of the given size and shape that uses `RealT` as coordinate type.

# Arguments
- `cells_per_dimension::NTupleE{NDIMS, Int}`: the number of cells in each dimension.
- `faces::NTuple{2*NDIMS}`: a tuple of `2 * NDIMS` functions that describe the faces of the domain.
                            Each function must take `NDIMS-1` arguments.
                            `faces[1]` describes the face onto which the face in negative x-direction
                            of the unit hypercube is mapped. The face in positive x-direction of
                            the unit hypercube will be mapped onto the face described by `faces[2]`.
                            `faces[3:4]` describe the faces in positive and negative y-direction respectively
                            (in 2D and 3D).
                            `faces[5:6]` describe the faces in positive and negative z-direction respectively (in 3D).
- `RealT::Type`: the type that should be used for coordinates.
- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for
                 each dimension if the boundaries in this dimension are periodic.
"""
function CurvedMesh(cells_per_dimension, faces::Tuple; RealT=Float64, periodicity=true)
  NDIMS = length(cells_per_dimension)

  validate_faces(faces)

  # Use the transfinite mapping with the correct number of arguments
  mapping = transfinite_mapping(faces)

  # Collect definitions of face functions in one string (separated by semicolons)
  face2substring(face) = code_string(face, ntuple(_ -> Float64, NDIMS-1))
  join_semicolon(strings) = join(strings, "; ")

  faces_definition = faces .|> face2substring .|> string |> join_semicolon

  # Include faces definition in `mapping_as_string` to allow for evaluation
  # without knowing the face functions
  mapping_as_string = "$faces_definition; mapping = transfinite_mapping(faces)"

  return CurvedMesh(cells_per_dimension, mapping; RealT=RealT, periodicity=periodicity, mapping_as_string=mapping_as_string)
end

"""
    CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

Create a CurvedMesh that represents a uncurved structured mesh with a rectangular domain.

# Arguments
- `cells_per_dimension::NTuple{NDIMS, Int}`: the number of cells in each dimension.
- `coordinates_min::NTuple{NDIMS, RealT}`: coordinate of the corner in the negative direction of each dimension.
- `coordinates_max::NTuple{NDIMS, RealT}`: coordinate of the corner in the positive direction of each dimension.
- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for
                 each dimension if the boundaries in this dimension are periodic.
"""
function CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max; periodicity=true)
  NDIMS = length(cells_per_dimension)
  RealT = promote_type(eltype(coordinates_min), eltype(coordinates_max))

  mapping = coordinates2mapping(coordinates_min, coordinates_max)
  mapping_as_string = "coordinates_min = $coordinates_min; " *
                      "coordinates_max = $coordinates_max; " *
                      "mapping = coordinates2mapping(coordinates_min, coordinates_max)"
  return CurvedMesh(cells_per_dimension, mapping; RealT=RealT, periodicity=periodicity, mapping_as_string=mapping_as_string)
end

# Extract a string of the code that defines the mapping function
mapping2string(mapping, ndims) = string(code_string(mapping, ntuple(_ -> Float64, ndims)))

# Interpolate linearly between left and right value where s should be between -1 and 1
linear_interpolate(s, left_value, right_value) = 0.5 * ((1 - s) * left_value + (1 + s) * right_value)

# Convert min and max coordinates of a rectangle to the corresponding transformation mapping
function coordinates2mapping(coordinates_min::NTuple{1}, coordinates_max::NTuple{1})
  mapping(xi) = linear_interpolate(xi, coordinates_min[1], coordinates_max[1])
end

function coordinates2mapping(coordinates_min::NTuple{2}, coordinates_max::NTuple{2})
  mapping(xi, eta) = SVector(linear_interpolate(xi,  coordinates_min[1], coordinates_max[1]),
                             linear_interpolate(eta, coordinates_min[2], coordinates_max[2]))
end

function coordinates2mapping(coordinates_min::NTuple{3}, coordinates_max::NTuple{3})
  mapping(xi, eta, zeta) = SVector(linear_interpolate(xi,   coordinates_min[1], coordinates_max[1]),
                                   linear_interpolate(eta,  coordinates_min[2], coordinates_max[2]),
                                   linear_interpolate(zeta, coordinates_min[3], coordinates_max[3]))
end

# In 1D
# Linear mapping from the reference element to the domain described by the faces
function linear_mapping(x, faces)
  return linear_interpolate(x, faces[1](), faces[2]())
end

# In 2D
# Bilinear mapping from the reference element to the domain described by the faces
function bilinear_mapping(x, y, faces)
  x1 = faces[1](-1) # Bottom left
  x2 = faces[2](-1) # Bottom right
  x3 = faces[1](1) # Top left
  x4 = faces[2](1) # Top right

  return 0.25 * (x1 * (1 - x) * (1 - y) +
                 x2 * (1 + x) * (1 - y) +
                 x3 * (1 - x) * (1 + y) +
                 x4 * (1 + x) * (1 + y))
end

# In 3D
# Trilinear mapping from the reference element to the domain described by the faces
function trilinear_mapping(x, y, z, faces)
  x1 = faces[1](-1, -1) # mapped from (-1,-1,-1)
  x2 = faces[2](-1, -1) # mapped from ( 1,-1,-1)
  x3 = faces[1]( 1, -1) # mapped from (-1, 1,-1)
  x4 = faces[2]( 1, -1) # mapped from ( 1, 1,-1)
  x5 = faces[1](-1,  1) # mapped from (-1,-1, 1)
  x6 = faces[2](-1,  1) # mapped from ( 1,-1, 1)
  x7 = faces[1]( 1,  1) # mapped from (-1, 1, 1)
  x8 = faces[2]( 1,  1) # mapped from ( 1, 1, 1)

  return 0.125 * (x1 * (1 - x) * (1 - y) * (1 - z) +
                  x2 * (1 + x) * (1 - y) * (1 - z) +
                  x3 * (1 - x) * (1 + y) * (1 - z) +
                  x4 * (1 + x) * (1 + y) * (1 - z) +
                  x5 * (1 - x) * (1 - y) * (1 + z) +
                  x6 * (1 + x) * (1 - y) * (1 + z) +
                  x7 * (1 - x) * (1 + y) * (1 + z) +
                  x8 * (1 + x) * (1 + y) * (1 + z) )
end

# Use linear mapping in 1D
transfinite_mapping(faces::NTuple{2, Any}) = x -> linear_mapping(x, faces)

# In 2D
# Transfinite mapping from the reference element to the domain described by the faces
function transfinite_mapping(faces::NTuple{4, Any})
  mapping(x, y) = (linear_interpolate(x, faces[1](y), faces[2](y)) +
                   linear_interpolate(y, faces[3](x), faces[4](x)) -
                   bilinear_mapping(x, y, faces))
end

# In 3D
# Correction term for the Transfinite mapping
function correction_term_3d(x, y, z, faces)
  # Correction for x-terms
  c_x = linear_interpolate(x, linear_interpolate(y, faces[3](-1, z), faces[4](-1, z)) +
                              linear_interpolate(z, faces[5](-1, y), faces[6](-1, y)),
                              linear_interpolate(y, faces[3]( 1, z), faces[4]( 1, z)) +
                              linear_interpolate(z, faces[5]( 1, y), faces[6]( 1, y)) )

  # Correction for y-terms
  c_y = linear_interpolate(y, linear_interpolate(x, faces[1](-1,  z), faces[2](-1,  z)) +
                              linear_interpolate(z, faces[5]( x, -1), faces[6]( x, -1)),
                              linear_interpolate(x, faces[1]( 1,  z), faces[2]( 1,  z)) +
                              linear_interpolate(z, faces[5]( x,  1), faces[6]( x,  1)) )

  # Correction for x-terms
  c_z = linear_interpolate(z, linear_interpolate(x, faces[1](y, -1), faces[2](y, -1)) +
                              linear_interpolate(y, faces[3](x, -1), faces[4](x, -1)),
                              linear_interpolate(x, faces[1](y,  1), faces[2](y,  1)) +
                              linear_interpolate(y, faces[3](x,  1), faces[4](x,  1)) )

  return 0.5 * (c_x + c_y + c_z)
end

# In 3D
# Transfinite mapping from the reference element to the domain described by the faces
function transfinite_mapping(faces::NTuple{6, Any})
  mapping(x, y, z) =  (linear_interpolate(x, faces[1](y, z), faces[2](y, z)) +
                       linear_interpolate(y, faces[3](x, z), faces[4](x, z)) +
                       linear_interpolate(z, faces[5](x, y), faces[6](x, y)) -
                       correction_term_3d(x, y, z, faces) +
                       trilinear_mapping(x, y, z, faces))
end

function validate_faces(faces::NTuple{2, Any}) end

function validate_faces(faces::NTuple{4, Any})
  @assert faces[1](-1) ≈ faces[3](-1) "faces[1](-1) needs to match faces[3](-1) (bottom left corner)"
  @assert faces[2](-1) ≈ faces[3]( 1) "faces[2](-1) needs to match faces[3](1) (bottom right corner)"
  @assert faces[1]( 1) ≈ faces[4](-1) "faces[1](1) needs to match faces[4](-1) (top left corner)"
  @assert faces[2]( 1) ≈ faces[4]( 1) "faces[2](1) needs to match faces[4](1) (top right corner)"
end

function validate_faces(faces::NTuple{6, Any})
  @assert (faces[1](-1, -1) ≈
           faces[3](-1, -1) ≈
           faces[5](-1, -1)) "faces[1](-1, -1), faces[3](-1, -1) and faces[5](-1, -1) need to match at (-1, -1, -1) corner"

  @assert (faces[2](-1, -1) ≈
           faces[3]( 1, -1) ≈
           faces[5]( 1, -1)) "faces[2](-1, -1), faces[3](1, -1) and faces[5](1, -1) need to match at (1, -1, -1) corner"

  @assert (faces[1]( 1, -1) ≈
           faces[4](-1, -1) ≈
           faces[5](-1,  1)) "faces[1](1, -1), faces[4](-1, -1) and faces[5](-1, 1) need to match at (-1, 1, -1) corner"

  @assert (faces[2]( 1, -1) ≈
           faces[4]( 1, -1) ≈
           faces[5]( 1,  1)) "faces[2](1, -1), faces[4](1, -1) and faces[5](1, 1) need to match at (1, 1, -1) corner"

  @assert (faces[1](-1,  1) ≈
           faces[3](-1,  1) ≈
           faces[6](-1, -1)) "faces[1](-1, 1), faces[3](-1, 1) and faces[6](-1, -1) need to match at (-1, -1, 1) corner"

  @assert (faces[2](-1,  1) ≈
           faces[3]( 1,  1) ≈
           faces[6]( 1, -1)) "faces[2](-1, 1), faces[3](1, 1) and faces[6](1, -1) need to match at (1, -1, 1) corner"

  @assert (faces[1]( 1,  1) ≈
           faces[4](-1,  1) ≈
           faces[6](-1,  1)) "faces[1](1, 1), faces[4](-1, 1) and faces[6](-1, 1) need to match at (-1, 1, 1) corner"

  @assert (faces[2]( 1,  1) ≈
           faces[4]( 1,  1) ≈
           faces[6]( 1,  1)) "faces[2](1, 1), faces[4](1, 1) and faces[6](1, 1) need to match at (1, 1, 1) corner"
end

# Check if mesh is periodic
isperiodic(mesh::CurvedMesh) = all(mesh.periodicity)
isperiodic(mesh::CurvedMesh, dimension) = mesh.periodicity[dimension]

@inline Base.ndims(::CurvedMesh{NDIMS}) where {NDIMS} = NDIMS
@inline Base.real(::CurvedMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
Base.size(mesh::CurvedMesh) = mesh.cells_per_dimension
Base.size(mesh::CurvedMesh, i) = mesh.cells_per_dimension[i]
Base.axes(mesh::CurvedMesh) = map(Base.OneTo, mesh.cells_per_dimension)
Base.axes(mesh::CurvedMesh, i) = Base.OneTo(mesh.cells_per_dimension[i])

function Base.show(io::IO, ::CurvedMesh{NDIMS, RealT}) where {NDIMS, RealT}
  print(io, "CurvedMesh{", NDIMS, ", ", RealT, "}")
end

function Base.show(io::IO, ::MIME"text/plain", mesh::CurvedMesh{NDIMS, RealT}) where {NDIMS, RealT}
  if get(io, :compact, false)
    show(io, mesh)
  else
    summary_header(io, "CurvedMesh{" * string(NDIMS) * ", " * string(RealT) * "}")
    summary_line(io, "size", size(mesh))

    summary_line(io, "mapping", "")
    # Print code lines of mapping_as_string
    mapping_lines = split(mesh.mapping_as_string, ";")
    for i in eachindex(mapping_lines)
      summary_line(increment_indent(io), "line $i", strip(mapping_lines[i]))
    end
    summary_footer(io)
  end
end


text:
# ---
# title: 1238. Circular Permutation in Binary Representation
# id: problem1238
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/circular-permutation-in-binary-representation/description/>
# hidden: true
# ---
#
# Given 2 integers `n` and `start`. Your task is return **any** permutation `p`
# of `(0,1,2.....,2^n -1) `such that :
#
#   * `p[0] = start`
#   * `p[i]` and `p[i+1]` differ by only one bit in their binary representation.
#   * `p[0]` and `p[2^n -1]` must also differ by only one bit in their binary representation.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 2, start = 3
#     Output: [3,2,0,1]
#     Explanation: The binary representation of the permutation is (11,10,00,01).
#     All the adjacent element differ by one bit. Another valid permutation is [3,1,0,2]
#
#
# **Example 2:**
#
#
#
#     Input: n = 3, start = 2
#     Output: [2,6,7,5,4,0,1,3]
#     Explanation: The binary representation of the permutation is (010,110,111,101,100,000,001,011).
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 16`
#   * `0 <= start < 2 ^ n`
#
#
## @lc code=start
using LeetCode

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


text:
### Thermal Generation Formulations

abstract type AbstractThermalFormulation end

abstract type AbstractThermalDispatchForm <: AbstractThermalFormulation end

abstract type AbstractThermalCommitmentForm <: AbstractThermalFormulation end

struct KenuvenThermalCommitment <: AbstractThermalCommitmentForm end

struct StandardThermalCommitment <: AbstractThermalCommitmentForm end

struct ThermalDispatch <: AbstractThermalDispatchForm end

struct ThermalRampLimitDispatch <: AbstractThermalDispatchForm end

include("thermal_generation/output_constraints.jl")
include("thermal_generation/ramping_constraints.jl")
include("thermal_generation/thermal_variables.jl")
include("thermal_generation/thermalgencommitment_cost.jl")
include("thermal_generation/thermalgenvariable_cost.jl")
include("thermal_generation/unitcommitment_constraints.jl")


text:
import Base.LinAlg.BLAS: asum, dot, blascopy!, nrm2, scal, scal!, gemm, gemm!, gemv, gemv!,
    syrk, symm, symm!, symv, symv!, trmm, trsm, trmv, trsv, trsv!, ger!

################################## Level 1 ##################################

# Unit-stride `dot`.
push!(ops, DiffOp(:(Base.dot),
    :(Tuple{DLA.SA{<:DLA.BF}, DLA.SA{<:DLA.BF}}),
    [true, true]
))
∇(::typeof(dot), ::Arg1, p, z::BF, z̄::BF, x::SA{<:BF}, y::SA{<:BF}) = z̄ .* y
∇(::typeof(dot), ::Arg2, p, z::BF, z̄::BF, x::SA{<:BF}, y::SA{<:BF}) = z̄ .* x
∇(x̄, ::typeof(dot), ::Arg1, p, z::BF, z̄::BF, x::SA{<:BF}, y::SA{<:BF}) = (x̄ .= x̄ .+ z̄ .* y)
∇(ȳ, ::typeof(dot), ::Arg2, p, z::BF, z̄::BF, x::SA{<:BF}, y::SA{<:BF}) = (ȳ .= ȳ .+ z̄ .* x)

# Arbitrary-stride `dot`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.dot),
    :(Tuple{Int, DLA.SA{<:DLA.BF}, Int, DLA.SA{<:DLA.BF}, Int}),
    [false, true, false, true, false],
))
∇(::typeof(dot), ::Arg2, p, z::BF, z̄::BF, n::Int, x::SA{<:BF}, ix::Int, y::SA{<:BF}, iy::Int) =
    scal!(n, z̄, blascopy!(n, y, iy, zeros(x), ix), ix)
∇(::typeof(dot), ::Arg4, p, z::BF, z̄::BF, n::Int, x::SA{<:BF}, ix::Int, y::SA{<:BF}, iy::Int) =
    scal!(n, z̄, blascopy!(n, x, ix, zeros(y), iy), iy)
∇(x̄, ::typeof(dot), ::Arg2, p, z::BF, z̄::BF, n::Int, x::SA{<:BF}, ix::Int, y::SA{<:BF}, iy::Int) =
    (x̄ .= x̄ .+ scal!(n, z̄, blascopy!(n, y, iy, zeros(x), ix), ix))
∇(ȳ, ::typeof(dot), ::Arg4, p, z::BF, z̄::BF, n::Int, x::SA{<:BF}, ix::Int, y::SA{<:BF}, iy::Int) =
    (ȳ .= ȳ .+ scal!(n, z̄, blascopy!(n, x, ix, zeros(y), iy), iy))

# Unit-stride `nrm2`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.nrm2),
    :(Tuple{DLA.SA{<:DLA.BF}}),
    [true]
))
∇(::typeof(nrm2), ::Arg1, p, y::BF, ȳ::BF, x::SA{<:BF}) = x .* (ȳ / y)
∇(x̄::AA, ::typeof(nrm2), ::Arg1, p, y::BF, ȳ::BF, x::SA{<:BF}) = (x̄ .= x̄ .+ x .* (ȳ / y))

# Arbitrary-stride `nrm2`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.nrm2),
    :(Tuple{Integer, DLA.SA{<:DLA.BF}, Integer}),
    [false, true, false]
))
∇(::typeof(nrm2), ::Arg2, p, y::BF, ȳ::BF, n::Integer, x::SA{<:BF}, inc::Integer) =
    scal!(n, ȳ / y, blascopy!(n, x, inc, zeros(x), inc), inc)
∇(x̄::SA{<:BF}, ::typeof(nrm2), ::Arg2, p, y::BF, ȳ::BF, n::Integer, x::SA{<:BF}, inc::Integer) =
    (x̄ .= x̄ .+ scal!(n, ȳ / y, blascopy!(n, x, inc, zeros(x), inc), inc))

# Unit-stride `asum`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.asum),
    :(Tuple{DLA.SA{<:DLA.BF}}),
    [true]
))
∇(::typeof(asum), ::Arg1, p, y::BF, ȳ::BF, x::SA{<:BF}) = ȳ .* sign.(x)
∇(x̄::AA, ::typeof(asum), ::Arg1, p, y::BF, ȳ::BF, x::SA{<:BF}) = (x̄ .= x̄ .+ ȳ .* sign.(x))

# Arbitrary-stride `asum`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.asum),
    :(Tuple{Integer, DLA.SA{<:DLA.BF}, Integer}),
    [false, true, false]
))
∇(::typeof(asum), ::Arg2, p, y::BF, ȳ::BF, n::Integer, x, inc::Integer) =
    scal!(n, ȳ, blascopy!(n, sign.(x), inc, zeros(x), inc), inc)
∇(x̄::SA{<:BF}, ::typeof(asum), ::Arg2, p, y::BF, ȳ::BF, n::Integer, x::SA{<:BF}, inc::Integer) =
    (x̄ .= x̄ .+ scal!(n, ȳ, blascopy!(n, sign.(x), inc, zeros(x), inc), inc))

# Some weird stuff going on that I haven't figured out yet. This is a very old attempt.
# let f = :(scal{T <: AbstractArray, V <: AbstractFloat})
#     ā = :(blascopy!(n, z̄, inc, zeros(X), inc) .* X)
#     X̄ = :(scal!(n, a, z̄, inc))
# end

################################## Level 2 ##################################

# `gemv` sensitivities implementation.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.gemv),
    :(Tuple{Char, T, DLA.SM{T}, DLA.SV{T}} where T<:DLA.BF),
    [false, true, true, true],
))
∇(::typeof(gemv), ::Arg2, _, y::SV{T}, ȳ::SV, tA::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    dot(ȳ, y) / α
∇(::typeof(gemv), ::Arg3, _, y::SV{T}, ȳ::SV, tA::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    uppercase(tA) == 'N' ? α * ȳ * x' : α * x * ȳ'
∇(Ā::SM{T}, ::typeof(gemv), ::Arg3, _, y::SV{T}, ȳ::SV{T}, tA::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    uppercase(tA) == 'N' ? ger!(α, ȳ, x, Ā) : ger!(α, x, ȳ, Ā)
∇(::typeof(gemv), ::Arg4, _, y::SV{T}, ȳ::SV{T}, tA::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    gemv(uppercase(tA) == 'N' ? 'T' : 'N', α, A, ȳ)
∇(x̄::SV{T}, ::typeof(gemv), ::Arg4, _, y::SV{T}, ȳ::SV{T}, tA::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    gemv!(uppercase(tA) == 'N' ? 'T' : 'N', α, A, ȳ, one(T), x̄)

# `gemv` sensitivities implementation with `α = 1`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.gemv),
    :(Tuple{Char, DLA.SM{T}, DLA.SV{T}} where T<:DLA.BF),
    [false, true, true],
))
∇(::typeof(gemv), ::Arg2, p, y::SV{T}, ȳ::SV{T}, tA::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(gemv, Val{3}, p, y, ȳ, tA, one(T), A, x)
∇(Ā::SM{T}, ::typeof(gemv), ::Arg2, p, y::SV{T}, ȳ::SV{T}, tA::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(Ā, gemv, Val{3}, p, y, ȳ, tA, one(T), A, x)
∇(::typeof(gemv), ::Arg3, p, y::SV{T}, ȳ::SV{T}, tA::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(gemv, Val{4}, p, y, ȳ, tA, one(T), A, x)
∇(x̄::SV{T}, ::typeof(gemv), ::Arg3, p, y::SV{T}, ȳ::SV{T}, tA::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(x̄, gemv, Val{4}, p, y, ȳ, tA, one(T), A, x)

# `symv` sensitivity implementations.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.symv),
    :(Tuple{Char, T, DLA.SM{T}, DLA.SV{T}} where T<:DLA.BF),
    [false, true, true, true],
))
∇(::typeof(symv), ::Arg2, p, y::SV{T}, ȳ::SV{T}, ul::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    dot(ȳ, y) / α
function ∇(::typeof(symv), ::Arg3, p, y::SV{T}, ȳ::SV{T}, ul::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF
    Y, Ȳ, X = reshape(y, length(y), 1), reshape(ȳ, length(ȳ), 1), reshape(x, length(x), 1)
    return ∇(symm, Val{4}, p, Y, Ȳ, 'L', ul, α, A, X)
end
function ∇(Ā::SM{T}, ::typeof(symv), ::Arg3, p, y::SV{T}, ȳ::SV{T}, ul::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF
    Y, Ȳ, X = reshape(y, length(y), 1), reshape(ȳ, length(ȳ), 1), reshape(x, length(x), 1)
    return ∇(Ā, symm, Val{4}, p, Y, Ȳ, 'L', ul, α, A, X)
end
∇(::typeof(symv), ::Arg4, p, y::SV{T}, ȳ::SV{T}, ul::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    symv(ul, α, A, ȳ)
∇(x̄::SV{T}, ::typeof(symv), ::Arg4, p, y::SV{T}, ȳ::SV{T}, ul::Char, α::T, A::SM{T}, x::SV{T}) where T<:BF =
    symv!(ul, α, A, ȳ, one(T), x̄)

# `symv` sensitivity implementations for `α=1`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.symv),
    :(Tuple{Char, DLA.SM{T}, DLA.SV{T}} where T<:DLA.BF),
    [false, true, true],
))
∇(::typeof(symv), ::Arg2, p, y::SV{T}, ȳ::SV{T}, ul::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(symv, Val{3}, p, y, ȳ, ul, one(T), A, x)
∇(Ā::SM{T}, ::typeof(symv), ::Arg2, p, y::SV{T}, ȳ::SV{T}, ul::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(Ā, symv, Val{3}, p, y::SV{T}, ȳ::SV{T}, ul, one(T), A, x)
∇(::typeof(symv), ::Arg3, p, y::SV{T}, ȳ::SV{T}, ul::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(symv, Val{4}, p, y, ȳ, ul, one(T), A, x)
∇(B̄::SV{T}, ::typeof(symv), ::Arg3, p, y::SV{T}, ȳ::SV{T}, ul::Char, A::SM{T}, x::SV{T}) where T<:BF =
    ∇(B̄, symv, Val{4}, p, y, ȳ, ul, one(T), A, x)

# `trmv` sensitivity implementations.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.trmv),
    :(Tuple{Char, Char, Char, DLA.SM{T}, DLA.SV{T}} where T<:DLA.BF),
    [false, false, false, true, true],
))
function ∇(::typeof(trmv), ::Arg4, p, y::SV{T}, ȳ::SV{T},
    ul::Char, ta::Char, dA::Char,
    A::SM{T},
    b::SV{T},
) where T<:BF
    Ā = (uppercase(ul) == 'L' ? tril! : triu!)(uppercase(ta) == 'N' ? ȳ * b' : b * ȳ')
    dA == 'U' && fill!(view(Ā, diagind(Ā)), zero(T))
    return Ā
end
∇(::typeof(trmv), ::Arg5, p, y::SV{T}, ȳ::SV{T},
    ul::Char, ta::Char, dA::Char,
    A::SM{T},
    b::SV{T},
) where T<:BF = trmv(ul, uppercase(ta) == 'N' ? 'T' : 'N', dA, A, ȳ)

# `trsv` sensitivity implementations.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.trsv),
    :(Tuple{Char, Char, Char, DLA.SM{T}, DLA.SV{T}} where T<:DLA.BF),
    [false, false, false, true, true],
))
function ∇(::typeof(trsv), ::Arg4, p, y::SV{T}, ȳ::SV{T},
    ul::Char, ta::Char, dA::Char,
    A::SM{T},
    x::SV{T},
) where T<:BF
    Y, Ȳ, X = reshape(y, length(y), 1), reshape(ȳ, length(ȳ), 1), reshape(x, length(x), 1)
    Ā = ∇(trsm, Val{6}, p, Y, Ȳ, 'L', ul, ta, dA, one(T), A, X)
    dA == 'U' && fill!(view(Ā, diagind(Ā)), zero(T))
    return Ā
end
∇(::typeof(trsv), ::Arg5, p, y::SV{T}, ȳ::SV{T},
    ul::Char, ta::Char, dA::Char,
    A::SM{T},
    x::SV{T},
) where T<:BF = trsv(ul, uppercase(ta) == 'N' ? 'T' : 'N', dA, A, ȳ)

################################## Level 3 ##################################

# `gemm` sensitivities implementation.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.gemm),
    :(Tuple{Char, Char, T, DLA.SM{T}, DLA.SM{T}} where T<:DLA.BF),
    [false, false, true, true, true],
))
∇(::typeof(gemm), ::Arg3, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF = sum(Ȳ .* Y) / α
∇(::typeof(gemm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF =
    uppercase(tA) == 'N' ?
        uppercase(tB) == 'N' ?
            gemm('N', 'T', α, Ȳ, B) :
            gemm('N', 'N', α, Ȳ, B) :
        uppercase(tB) == 'N' ?
            gemm('N', 'T', α, B, Ȳ) :
            gemm('T', 'T', α, B, Ȳ)
∇(Ā::SM{T}, ::typeof(gemm), ::Arg4, _, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF =
    uppercase(tA) == 'N' ?
        uppercase(tB) == 'N' ?
            gemm!('N', 'T', α, Ȳ, B, one(T), Ā) :
            gemm!('N', 'N', α, Ȳ, B, one(T), Ā) :
        uppercase(tB) == 'N' ?
            gemm!('N', 'T', α, B, Ȳ, one(T), Ā) :
            gemm!('T', 'T', α, B, Ȳ, one(T), Ā)
∇(::typeof(gemm), ::Arg5, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF =
    uppercase(tA) == 'N' ?
        uppercase(tB) == 'N' ?
            gemm('T', 'N', α, A, Ȳ) :
            gemm('T', 'N', α, Ȳ, A) :
        uppercase(tB) == 'N' ?
            gemm('N', 'N', α, A, Ȳ) :
            gemm('T', 'T', α, Ȳ, A)
∇(B̄::SM{T}, ::typeof(gemm), ::Arg5, _, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF =
    uppercase(tA) == 'N' ?
        uppercase(tB) == 'N' ?
            gemm!('T', 'N', α, A, Ȳ, one(T), B̄) :
            gemm!('T', 'N', α, Ȳ, A, one(T), B̄) :
        uppercase(tB) == 'N' ?
            gemm!('N', 'N', α, A, Ȳ, one(T), B̄) :
            gemm!('T', 'T', α, Ȳ, A, one(T), B̄)

# `gemm` sensitivities implementation for `α = 1`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.gemm),
    :(Tuple{Char, Char, DLA.SM{T}, DLA.SM{T}} where T<:DLA.BF),
    [false, false, true, true],
))
∇(::typeof(gemm), ::Arg3, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    A::SM{T},
    B::SM{T}
) where T<:BF = ∇(gemm, Val{4}, p, Y, Ȳ, tA, tB, one(T), A, B)
∇(Ā::SM{T}, ::typeof(gemm), ::Arg3, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(Ā, gemm, Val{4}, p, Y, Ȳ, tA, tB, one(T), A, B)
∇(::typeof(gemm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(gemm, Val{5}, p, Y, Ȳ, tA, tB, one(T), A, B)
∇(B̄::SM{T}, ::typeof(gemm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    tA::Char,
    tB::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(B̄, gemm, Val{5}, p, Y, Ȳ, tA, tB, one(T), A, B)

# # `syrk` sensitivity implementations.
# @explicit_intercepts(
#     syrk,
#     Tuple{Char, Char, ∇Scalar, StridedVecOrMat{<:∇Scalar}},
#     [false, false, true, true],
# )
# function ∇(::typeof(syrk), ::Type{Arg{3}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     α::∇Scalar,
#     A::StridedVecOrMat{<:∇Scalar},
# )
#     g! = uppercase(uplo) == 'L' ? tril! : triu!
#     return sum(g!(Ȳ .* Y)) / α
# end
# function ∇(::typeof(syrk), ::Type{Arg{4}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     α::∇Scalar,
#     A::StridedVecOrMat{<:∇Scalar},
# )
#     triȲ = uppercase(uplo) == 'L' ? tril(Ȳ) : triu(Ȳ)
#     out = gemm('N', trans, α, triȲ .+ triȲ', A)
#     return uppercase(trans) == 'N' ? out : out'
# end
# function ∇(Ā::StridedVecOrMat{T}, ::typeof(syrk), ::Type{Arg{4}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     α::∇Scalar,
#     A::StridedVecOrMat{T},
# ) where T<:∇Scalar
#     triȲ = uppercase(uplo) == 'L' ? tril(Ȳ) : triu(Ȳ)
#     out = gemm('N', trans, α, triȲ .+ triȲ', A)
#     return broadcast!((ā, δā)->ā+δā, Ā, Ā, uppercase(trans) == 'N' ? out : out')
# end

# # `syrk` sensitivity implementations for `α=1`.
# @explicit_intercepts(
#     syrk,
#     Tuple{Char, Char, StridedVecOrMat{<:∇Scalar}},
#     [false, false, true],
# )
# ∇(::typeof(syrk), ::Type{Arg{3}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     A::StridedVecOrMat{<:∇Scalar},
# ) = ∇(syrk, Arg{4}, p, Y, Ȳ, uplo, trans, one(eltype(A)), A)
# ∇(Ā::StridedVecOrMat{T}, ::typeof(syrk), ::Type{Arg{4}}, p, Y, Ȳ,
#     uplo::Char,
#     trans::Char,
#     A::StridedVecOrMat{T},
# ) where T<:∇Scalar = ∇(Ā, syrk, Arg{4}, p, Y, Ȳ, uplo, char, one(eltype(A)), A)

# `symm` sensitivity implementations.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.symm),
    :(Tuple{Char, Char, T, DLA.SM{T}, DLA.SM{T}} where T<:DLA.BF),
    [false, false, true, true, true],
))
∇(::typeof(symm), ::Arg3, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF = sum(Ȳ .* Y) / α
function ∇(::typeof(symm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF
    tmp = uppercase(side) == 'L' ? Ȳ * B' : B'Ȳ
    g! = uppercase(ul) == 'L' ? tril! : triu!
    return α * g!(tmp + tmp' - Diagonal(tmp))
end
function ∇(Ā::SM{T}, ::typeof(symm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF
    tmp = uppercase(side) == 'L' ? Ȳ * B' : B'Ȳ
    g! = uppercase(ul) == 'L' ? tril! : triu!
    return broadcast!((ā, δā)->ā + δā, Ā, Ā, α * g!(tmp + tmp' - Diagonal(tmp)))
end
∇(::typeof(symm), ::Arg5, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF = symm(side, ul, α, A, Ȳ)
∇(B̄::SM{T}, ::typeof(symm), ::Arg5, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF = symm!(side, ul, α, A, Ȳ, one(T), B̄)

# `symm` sensitivity implementations for `α=1`.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.symm),
    :(Tuple{Char, Char, DLA.SM{T}, DLA.SM{T}} where T<:DLA.BF),
    [false, false, true, true],
))
∇(::typeof(symm), ::Arg3, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(symm, Val{4}, p, Y, Ȳ, side, ul, one(T), A, B)
∇(Ā::SM{T}, ::typeof(symm), ::Arg3, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(Ā, symm, Val{4}, p, Y, Ȳ, side, ul, one(T), A, B)
∇(::typeof(symm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(symm, Val{5}, p, Y, Ȳ, side, ul, one(T), A, B)
∇(B̄::SM{T}, ::typeof(symm), ::Arg4, p, Y::SM{T}, Ȳ::SM{T},
    side::Char,
    ul::Char,
    A::SM{T},
    B::SM{T},
) where T<:BF = ∇(B̄, symm, Val{5}, p, Y, Ȳ, side, ul, one(T), A, B)

# `trmm` sensitivity implementations.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.trmm),
    :(Tuple{Char, Char, Char, Char, T, DLA.SM{T}, DLA.SM{T}} where T<:DLA.BF),
    [false, false, false, false, true, true, true],
))
∇(::typeof(trmm), ::Arg5, p, Y::SM{T}, Ȳ::SM{T},
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF = sum(Ȳ .* Y) / α
function ∇(::typeof(trmm), ::Arg6, p, Y::SM{T}, Ȳ::SM{T},
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF
    Ā_full = uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            gemm('N', 'T', α, Ȳ, B) :
            gemm('N', 'T', α, B, Ȳ) :
        uppercase(ta) == 'N' ?
            gemm('T', 'N', α, B, Ȳ) :
            gemm('T', 'N', α, Ȳ, B)
    dA == 'U' && fill!(view(Ā_full, diagind(Ā_full)), zero(T))
    return (uppercase(ul) == 'L' ? tril! : triu!)(Ā_full)
end
∇(::typeof(trmm), ::Type{Val{7}}, p, Y::SM{T}, Ȳ::SM{T},
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::SM{T},
    B::SM{T},
) where T<:BF =
    uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            trmm('L', ul, 'T', dA, α, A, Ȳ) :
            trmm('L', ul, 'N', dA, α, A, Ȳ) :
        uppercase(ta) == 'N' ?
            trmm('R', ul, 'T', dA, α, A, Ȳ) :
            trmm('R', ul, 'N', dA, α, A, Ȳ)

# `trsm` sensitivity implementations.
push!(ops, DiffOp(:(Base.LinAlg.BLAS.trsm),
    :(Tuple{Char, Char, Char, Char, T, DLA.SM{T}, DLA.SM{T}} where T<:DLA.BF),
    [false, false, false, false, true, true, true],
))
∇(::typeof(trsm), ::Arg5, p, Y::SM{T}, Ȳ::SM{T},
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::SM{T},
    X::SM{T},
) where T<:BF = sum(Ȳ .* Y) / α
function ∇(::typeof(trsm), ::Arg6, p, Y::SM{T}, Ȳ::SM{T},
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::SM{T},
    X::SM{T},
) where T<:BF
    Ā_full = uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            trsm('L', ul, 'T', dA, -one(T), A, Ȳ * Y') :
            trsm('R', ul, 'T', dA, -one(T), A, Y * Ȳ') :
        uppercase(ta) == 'N' ?
            trsm('R', ul, 'T', dA, -one(T), A, Y'Ȳ) :
            trsm('L', ul, 'T', dA, -one(T), A, Ȳ'Y)
    dA == 'U' && fill!(view(Ā_full, diagind(Ā_full)), zero(T))
    return (uppercase(ul) == 'L' ? tril! : triu!)(Ā_full)
end
∇(::typeof(trsm), ::Type{Val{7}}, p, Y::SM{T}, Ȳ::SM{T},
    side::Char, ul::Char, ta::Char, dA::Char,
    α::T,
    A::SM{T},
    X::SM{T},
) where T<:BF =
    uppercase(side) == 'L' ?
        uppercase(ta) == 'N' ?
            trsm('L', ul, 'T', dA, α, A, Ȳ) :
            trsm('L', ul, 'N', dA, α, A, Ȳ) :
        uppercase(ta) == 'N' ?
            trsm('R', ul, 'T', dA, α, A, Ȳ) :
            trsm('R', ul, 'N', dA, α, A, Ȳ)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Zenith"
version = v"0.12.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/bvaisvil/zenith/archive/refs/tags/$(version).tar.gz",
                  "2b33892be95149550c84179b341e304c4222e3489bc121ea8c8346e075433aa6"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/zenith*/
# Get rid of misleading settings
rm -rf .cargo
cargo build --release
install -Dvm 755 "target/${rust_target}/release/zenith${exeext}" "${bindir}/zenith${exeext}"
"""

# This works only on Unix systems, but not FreeBSD:
# https://github.com/bvaisvil/zenith/issues/8
platforms = supported_platforms(; exclude=p -> Sys.iswindows(p) || Sys.isfreebsd(p))

# The products that we will ensure are always built
products = [
    ExecutableProduct("zenith", :zenith),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", compilers=[:c, :rust])

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
function print_matrix(M, precision)
  # This function prints a 2D matrix
  a, b = size(M)
  println("Hamiltonian:")
  for i=1:b
  	for j=1:a
  		print(round(M[j, i], precision), "\t")
  	end
  	println()
  end
end


text:
module KernelPerf
import Perftests: @perf, meta

abstract List{T}

type Nil{T} <: List{T}
end

type Cons{T} <: List{T}
    head::T
    tail::List{T}
end

function listn1n2(n1::Int,n2::Int)
    l1 = Nil{Int}()
    for i=n2:-1:n1
        l1 = Cons{Int}(i,l1)
    end
    l1
end

@perf listn1n2(1,10^6) meta("cons", "List concatenation")

# issue #1211
include("ziggurat.jl")
a = Array(Float64, 1000000)
@perf randn_zig!(a) meta("randn_zig", "Ziggurat gaussian number generator", issue=1211)

# issue #950
include("gk.jl")
@perf gk(350,[0.1]) meta("gk", "Grigoriadis Khachiyan matrix games", issue=950)

# issue #942
x = sparse(ones(280,280));
@perf (x*x) meta("sparsematmul", "normal", "Sparse matrix multiplication", issue=942)
s2 = sparse(rand(1:2000,10^5), kron([1:10^4;],ones(Int,10)), ones(Int,10^5), 2000, 10^4);
@perf s2*s2' meta("sparsematmul", "fillin", "Sparse matrix multiplication with fill-in", issue=942)

# issue #938
x = 1:600000;
@perf sparse(x,x,x) meta("sparserange", "Construction of a sparse matrix from ranges", issue=938)

# issue 4707
include("getdivgrad.jl")
A = getDivGrad(64,64,64)
v = rand(64^3)
@perf A*v meta("sparsemat_densevec", "Sparse matrix - dense vector multiplication", issue=4707)

# issue #939
y = [500000:-1:1;]
@perf sortperm(y) meta("sortperm", "Sorting of a worst-case vector", issue=939)

# issue #445
include("stockcorr.jl")
@perf stockcorr() meta("stockcorr", "Correlation analysis of random matrices", issue=445)

include("bench_eu.jl")
@perf bench_eu_vec(10000) meta("finance", "vectorized", "Vectorized Monte Carlo financial simulation")
@perf bench_eu_devec(10000) meta("finance", "devectorized", "Devectorized Monte Carlo financial simulation")

# issue #1163
include("actor_centrality.jl")
@perf actor_centrality() meta("actorgraph", "Centrality of actors in IMDB database", issue=1163)

# issue #1168
include("laplace.jl")
@perf laplace_vec() meta("laplace", "vectorized", "Vectorized Laplacian", issue=1168)
@perf laplace_devec() meta("laplace", "devectorized", "Devectorized Laplacian", issue=1168)

# issue #1169
include("go_benchmark.jl")
@perf benchmark(10) meta("go","Simulation of random games of Go", issue=1169)

# issue #3142
include("simplex.jl")
@perf doTwoPassRatioTest() meta("simplex", "Dual simplex algorithm for Linear Programming", issue=3142)

# issue #3811
include("raytracer.jl")
@perf Raytracer(5, 256, 4) meta("raytracer", "Simple native Julia raytracer", issue=3811)

function cmp_with_func(x::Vector, f::Function)
    count::Int = 0
    for i = 1:length(x)-1
      if f(x[i], x[i+1]) count += 1 end
    end
    count
end
x = randn(200_000)
@perf cmp_with_func(x, isless) meta("funarg", "Function argument benchmark")

arith_vectorized(b,c,d) = b.*c + d .+ 1.0
len = 1_000_000
b = randn(len)
c = randn(len)
d = randn(len)
@perf arith_vectorized(b,c,d) meta("vectorize", "Vectorized arithmetic")

open("random.csv","w") do io
    writecsv(io, rand(100000,4))
end
function parsecsv()
    for line in EachLine(open("random.csv"))
        split(line, ',')
    end
end
@perf parsecsv() meta("splitline", "CSV parsing")
rm("random.csv")

include("json.jl")
@perf parse_json(_json_data) meta("json", "JSON parsing")

include("indexing.jl")
x = [1:100_000;]
y = filter(iseven, 1:length(x))
logical_y = map(iseven, 1:length(x))
@perf add1!(x,y) meta("indexing", "add1", "Increment vector x at locations y")
@perf devec_add1!(x,y) meta("indexing", "devectorized_add1", "Devectorized increment vector x at locations y")
@perf add1!(x,logical_y) meta("indexing", "add1_logical", "Increment x_i if y_i is true")
@perf devec_add1_logical!(x,logical_y) meta("indexing", "devectorized_add1_logical", "Devectorized increment x_i if y_i is true")

end #module


text:
using Test, Base.CoreLogging
import Base.CoreLogging: Info

import InteractiveUtils
InteractiveUtils.versioninfo(verbose=true)

using GPUCompiler

using LLVM, LLVM.Interop

include("util.jl")

@testset "GPUCompiler" begin

GPUCompiler.reset_runtime()

GPUCompiler.enable_timings()

Base.libllvm_version < v"12" && include("native.jl")    # TODO(vchuravy): wrap ORCv2
include("ptx.jl")
Base.libllvm_version < v"12" && include("spirv.jl")     # TODO(maleadt): rebuild SPIRV tools
Base.libllvm_version < v"12" && include("gcn.jl")       # TODO(jpsamaroo): test failures
include("bpf.jl")

include("examples.jl")

haskey(ENV, "CI") && GPUCompiler.timings()

end


text:
include("abstract.jl")

# For recording error information
mutable struct MagnitudeStats
    max::Float64
    min::Float64
    sum::Float64
    count::Int64
end

MagnitudeStats() = MagnitudeStats(-Inf, +Inf, 0.0, 0)

function Base.push!(stat::MagnitudeStats, value)
    v = log10(abs(value))
    stat.max = max(stat.max, v)
    stat.min = min(stat.min, v)
    stat.sum += v
    stat.count += 1
end

Base.min(s::MagnitudeStats) = s.count > 0 ? 10.0^(s.min) : 0.0
Base.max(s::MagnitudeStats) = s.count > 0 ? 10.0^(s.max) : 0.0
Statistics.mean(s::MagnitudeStats) = s.count > 0 ? 10.0^(s.sum / s.count) : 0.0
Base.length(s::MagnitudeStats) = s.count

function Base.show(io::IO, s::MagnitudeStats)
    println(io, "MagnitudeStats: ($(s.count) Values)")
    println(io, "\tmin = $(min(s))")
    println(io, "\tmean = $(mean(s))")
    println(io, "\tmax = $(max(s))")
end

"""
Analysis data of determinant quantum Monte Carlo (DQMC) simulation
"""
@with_kw mutable struct DQMCAnalysis
    acc_rate::Float64 = 0.
    prop_local::Int = 0
    acc_local::Int = 0
    acc_rate_global::Float64 = 0.
    prop_global::Int = 0
    acc_global::Int = 0

    imaginary_probability::MagnitudeStats = MagnitudeStats()
    negative_probability::MagnitudeStats = MagnitudeStats()
    propagation_error::MagnitudeStats = MagnitudeStats()
end

"""
Parameters of determinant quantum Monte Carlo (DQMC)
"""
struct DQMCParameters
    global_moves::Bool
    global_rate::Int

    thermalization::Int
    sweeps::Int

    silent::Bool
    check_sign_problem::Bool
    check_propagation_error::Bool

    safe_mult::Int
    delta_tau::Float64
    beta::Float64
    slices::Int

    measure_rate::Int
end

function DQMCParameters(;
        global_moves::Bool  = false,
        global_rate::Int    = 5,
        thermalization::Int = 100,
        sweeps::Int         = 100,
        silent::Bool        = false,
        check_sign_problem::Bool = true,
        check_propagation_error::Bool = true,
        safe_mult::Int      = 10,
        measure_rate::Int   = 10,
        warn_round::Bool    = true,
        kwargs...
    )
    nt = (;kwargs...)
    keys(nt) == (:beta,) && (nt = (;beta=nt.beta, delta_tau=0.1))
    @assert length(nt) >= 2 "Invalid keyword arguments to DQMCParameters: $nt"
    if (Set ∘ keys)(nt) == Set([:delta_tau, :beta, :slices])
        delta_tau, beta = nt.delta_tau, nt.beta
        slices = round(Int, beta/delta_tau)
        if slices != nt.slices
            error(
                "Given slices ($(nt.slices)) does not match calculated slices" *
                " beta/delta_tau ≈ $(slices)"
            )
        end
    elseif (Set ∘ keys)(nt) == Set([:beta, :slices])
        beta, slices = nt.beta, nt.slices
        delta_tau = beta / slices
    elseif (Set ∘ keys)(nt) == Set([:delta_tau, :slices])
        delta_tau, slices = nt.delta_tau, nt.slices
        beta = delta_tau * slices
    elseif (Set ∘ keys)(nt) == Set([:delta_tau, :beta])
        delta_tau, beta = nt.delta_tau, nt.beta
        slices = round(beta/delta_tau)
        if warn_round && !(slices ≈ beta/delta_tau)
            @warn "beta/delta_tau = $(beta/delta_tau) not an integer. Rounded to $slices"
        end
    else
        error("Invalid keyword arguments to DQMCParameters $nt")
    end
    DQMCParameters(
        global_moves,
        global_rate,
        thermalization,
        sweeps,
        silent,
        check_sign_problem,
        check_propagation_error,
        safe_mult,
        delta_tau,
        beta,
        slices,
        measure_rate
    )
end

"""
Determinant quantum Monte Carlo (DQMC) simulation
"""
mutable struct DQMC{
        M <: Model, CB <: Checkerboard, ConfType <: Any, RT <: AbstractRecorder,
        Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack
    } <: MonteCarloFlavor

    model::M
    conf::ConfType
    last_sweep::Int

    s::Stack
    ut_stack::UTStack
    p::DQMCParameters
    a::DQMCAnalysis

    configs::RT
    thermalization_measurements::Dict{Symbol, AbstractMeasurement}
    measurements::Dict{Symbol, AbstractMeasurement}

    function DQMC{M, CB, ConfType, RT, Stack, UTStack}(; kwargs...) where {
            M <: Model, CB <: Checkerboard, ConfType <: Any,
            RT <: AbstractRecorder,
            Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack
        }
        complete!(new{M, CB, ConfType, RT, Stack, UTStack}(), kwargs)
    end
    function DQMC(CB::Type{<:Checkerboard}; kwargs...)
        DQMC{
            Model, CB, Any, AbstractRecorder,
            AbstractDQMCStack, AbstractDQMCStack
        }(; kwargs...)
    end
    function DQMC{CB}(
            model::M,
            conf::ConfType,
            last_sweep::Int,
            s::Stack,
            ut_stack::UTStack,
            p::DQMCParameters,
            a::DQMCAnalysis,
            configs::RT,
            thermalization_measurements::Dict{Symbol, AbstractMeasurement},
            measurements::Dict{Symbol, AbstractMeasurement}
        ) where {
            M <: Model, CB <: Checkerboard, ConfType <: Any,
            RT <: AbstractRecorder,
            Stack <: AbstractDQMCStack, UTStack <: AbstractDQMCStack
        }
        new{M, CB, ConfType, RT, Stack, UTStack}(
            model, conf, last_sweep, s, ut_stack, p, a, configs,
            thermalization_measurements, measurements
        )
    end
    function DQMC{M, CB, C, RT, S, UTS}(args...) where {M, CB, C, RT, S, UTS}
        DQMC{CB}(args...)
    end
end

function complete!(a::DQMC, kwargs)
    for (field, val) in kwargs
        setfield!(a, field, val)
    end
    make_concrete!(a)
end

function make_concrete!(a::DQMC{M, CB, C, RT, S, UTS}) where {M, CB, C, RT, S, UTS}
    Ts = (
        isdefined(a, :model) ? typeof(a.model) : M,
        CB,
        isdefined(a, :conf) ? typeof(a.conf) : C,
        isdefined(a, :configs) ? typeof(a.configs) : RT,
        isdefined(a, :s) ? typeof(a.s) : S,
        isdefined(a, :ut_stack) ? typeof(a.ut_stack) : UTS
    )
    all(Ts .== (M, CB, C, RT, S, UTS)) && return a
    data = [(f, getfield(a, f)) for f in fieldnames(DQMC) if isdefined(a, f)]
    DQMC{Ts...}(; data...)
end

include("stack.jl")
include("unequal_time_stack.jl")
include("slice_matrices.jl")

"""
    DQMC(m::M; kwargs...) where M<:Model

Create a determinant quantum Monte Carlo simulation for model `m` with
keyword parameters `kwargs`.

### Keyword Arguments:
- `seed`: The random seed used by the simulation.
- `checkerboard=false`: If true, the simulation uses a generic checkerboard
decomposition.
- `thermalization_measurements::Dict{Symbol, AbstractMeasurement}`: A collection
of measurements run during the thermalization stage. By default, none are used.
- `measurements::Dict{Symbol, AbstractMeasurement}`: A collection of measurements
run during the measurement stage. Calls `default_measurements` if not specified.
- `recorder = ConfigRecorder`: Type of recorder used for saving configurations
generated during the simulation. Used (by default) when `replay!`ing simulations.
 (`Discarder` or `ConfigRecorder`)
- `recording_rate = measure_rate`: Rate at which configurations are recorded.
- `thermalization = 100`: Number of thermalization sweeps
- `sweeps`: Number of measurement sweeps
- `all_checks = true`: Check for Propagation instabilities and sign problems.
- `safe_mult = 10`: Number of "safe" matrix multiplications. Every `safe_mult`
multiplications, a UDT decomposition is used to stabilize the product.
- `delta_tau = 0.1`: Time discretization of the path integral
- `beta::Float64`: Inverse temperature used in the simulation
- `slices::Int = beta / delta_tau`: Number of imaginary time slice in the
simulation
- `measure_rate = 10`: Number of sweeps discarded between every measurement.
- `global_moves = false`:: Currently not used
- `global_rate = 5`: Currently not used
- `last_sweep = 0`: Sets the index of the last finished sweep. The simulation
will start with sweep `last_sweep + 1`.
"""
function DQMC(m::M;
        seed::Int=-1,
        checkerboard::Bool=false,
        thermalization_measurements = Dict{Symbol, AbstractMeasurement}(),
        measurements = Dict{Symbol, AbstractMeasurement}(),
        last_sweep = 0,
        recorder = ConfigRecorder,
        measure_rate = 10,
        recording_rate = measure_rate,
        kwargs...
    ) where M<:Model
    # default params
    # paramskwargs = filter(kw->kw[1] in fieldnames(DQMCParameters), kwargs)
    p = DQMCParameters(measure_rate = measure_rate; kwargs...)

    HET = hoppingeltype(DQMC, m)
    GET = greenseltype(DQMC, m)
    HMT = hopping_matrix_type(DQMC, m)
    GMT = greens_matrix_type(DQMC, m)
    IMT = interaction_matrix_type(DQMC, m)
    stack = DQMCStack{GET, HET, GMT, HMT, IMT}()
    ut_stack = UnequalTimeStack{GET, GMT}()

    conf = rand(DQMC, m, p.slices)
    analysis = DQMCAnalysis()
    CB = checkerboard ? CheckerboardTrue : CheckerboardFalse

    mc = DQMC{M, CB, typeof(conf), recorder, DQMCStack, AbstractDQMCStack}(
        model = m, conf = conf, last_sweep = last_sweep, s = stack,
        ut_stack = ut_stack, p = p, a = analysis,
        thermalization_measurements = thermalization_measurements,
        measurements = measurements
    )

    mc.configs = recorder(mc, m, recording_rate)
    seed == -1 || Random.seed!(seed)
    mc.conf = conf
    init!(mc)
    return make_concrete!(mc)
end

"""
    DQMC(m::M, params::Dict)
    DQMC(m::M, params::NamedTuple)

Create a determinant quantum Monte Carlo simulation for model `m` with
(keyword) parameters as specified in the dictionary/named tuple `params`.
"""
DQMC(m::Model, params::Dict{Symbol, T}) where T = DQMC(m; params...)
DQMC(m::Model, params::NamedTuple) = DQMC(m; params...)

# convenience
@inline beta(mc::DQMC) = mc.p.beta
@inline nslices(mc::DQMC) = mc.p.slices
@inline model(mc::DQMC) = mc.model
@inline conf(mc::DQMC) = mc.conf
@inline current_slice(mc::DQMC) = mc.s.current_slice
@inline last_sweep(mc::DQMC) = mc.last_sweep
@inline configurations(mc::DQMC) = mc.configs
@inline lattice(mc::DQMC) = lattice(mc.model)
@inline parameters(mc::DQMC) = merge(parameters(mc.p), parameters(mc.model))
@inline parameters(p::DQMCParameters) = (
    beta = p.beta, delta_tau = p.delta_tau, thermalization = p.thermalization, sweeps = p.sweeps
)

# cosmetics
import Base.summary
import Base.show
Base.summary(mc::DQMC) = "DQMC simulation of $(summary(mc.model))"
function Base.show(io::IO, mc::DQMC)
    print(io, "Determinant quantum Monte Carlo simulation\n")
    print(io, "Model: ", mc.model, "\n")
    print(io, "Beta: ", mc.p.beta, " (T ≈ $(round(1/mc.p.beta, sigdigits=3)))\n")
    N_th_meas = length(mc.thermalization_measurements)
    N_me_meas = length(mc.measurements)
    print(io, "Measurements: ", N_th_meas + N_me_meas, " ($N_th_meas + $N_me_meas)")
end
Base.show(io::IO, m::MIME"text/plain", mc::DQMC) = print(io, mc)

function ConfigRecorder(mc::DQMC, model::Model, rate = 10)
    ConfigRecorder{typeof(compress(mc, model, conf(mc)))}(rate)
end

function init!(mc::DQMC)
    init_hopping_matrices(mc, mc.model)
    initialize_stack(mc, mc.s)
    nothing
end
@deprecate resume_init!(mc::DQMC) init!(mc) false

"""
    run!(mc::DQMC[; kwargs...])

Runs the given Monte Carlo simulation `mc`. Returns true if the run finished and
false if it cancelled early to generate a resumable save-file.

### Keyword Arguments:
- `verbose = true`: If true, print progress messaged to stdout.
- `thermalization`: Number of thermalization sweeps. Uses the value passed to
`DQMC` by default.
- `sweeps`: Number of measurement sweeps. Uses the value passed to `DQMC` by
default.
- `safe_every::TimePeriod`: Set the interval for regularly scheduled saves.
- `safe_before::Date`: If this date is passed, `run!` will generate a resumable
save file and exit
- `grace_period = Minute(5)`: Buffer between the current time and `safe_before`.
The time required to generate a save file should be included here.
- `resumable_filename`: Name of the resumable save file. The default is based on
`safe_before`.
- `overwrite = false`: If set to true a file with the same name as
`resumable_filename` will be overwritten. (This will create a temporary backup)

See also: [`resume!`](@ref)
"""
@bm function run!(
        mc::DQMC;
        verbose::Bool = true,
        sweeps::Int = mc.p.sweeps,
        thermalization = mc.p.thermalization,
        safe_before::TimeType = now() + Year(100),
        safe_every::TimePeriod = Hour(10000),
        grace_period::TimePeriod = Minute(5),
        resumable_filename::String = "resumable_$(Dates.format(safe_before, "d_u_yyyy-HH_MM")).jld",
        overwrite = false
    )

    # Update number of sweeps
    if (mc.p.thermalization != thermalization) || (mc.p.sweeps != sweeps)
        verbose && println("Rebuilding DQMCParameters with new number of sweeps.")
        p = DQMCParameters(
            mc.p.global_moves,
            mc.p.global_rate,
            thermalization,
            sweeps,
            mc.p.silent,
            mc.p.check_sign_problem,
            mc.p.check_propagation_error,
            mc.p.safe_mult,
            mc.p.delta_tau,
            mc.p.beta,
            mc.p.slices,
            mc.p.measure_rate
        )
        mc.p = p
    end
    total_sweeps = sweeps + thermalization

    # Generate measurement groups
    groups = generate_groups(mc, mc.model, collect(values(mc.measurements)))

    start_time = now()
    last_checkpoint = now()
    max_sweep_duration = 0.0
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    # fresh stack
    verbose && println("Preparing Green's function stack")
    init!(mc)
    build_stack(mc, mc.s)
    propagate(mc)

    _time = time()
    verbose && println("\n\nThermalization stage - ", thermalization)
    # do_th_measurements && prepare!(mc.thermalization_measurements, mc, mc.model)
    # mc.last_sweep:total_sweeps won't change when last_sweep is changed
    for i in mc.last_sweep+1:total_sweeps
        verbose && (i == thermalization + 1) && println("\n\nMeasurement stage - ", sweeps)
        for u in 1:2 * nslices(mc)
            update(mc, i)

            if current_slice(mc) == 1 && i ≤ thermalization &&
                mc.s.direction == 1 && iszero(i % mc.p.measure_rate)
                measure!(mc.thermalization_measurements, mc, mc.model, i)
            end
            if current_slice(mc) == 1 && mc.s.direction == 1 && i > thermalization
                push!(mc.configs, mc, mc.model, i)
                if iszero(i % mc.p.measure_rate)
                    for (requirement, group) in groups
                        apply!(requirement, group, mc, mc.model, i)
                    end
                end
            end
        end
        mc.last_sweep = i

        if mod(i, 10) == 0
            mc.a.acc_rate = mc.a.acc_rate / (10 * 2 * nslices(mc))
            mc.a.acc_rate_global = mc.a.acc_rate_global / (10 / mc.p.global_rate)
            sweep_dur = (time() - _time)/10
            max_sweep_duration = max(max_sweep_duration, sweep_dur)
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
                @printf("\t\tacc rate (local) : %.1f%%\n", mc.a.acc_rate*100)
                if mc.p.global_moves
                  @printf("\t\tacc rate (global): %.1f%%\n", mc.a.acc_rate_global*100)
                  @printf("\t\tacc rate (global, overall): %.1f%%\n",
                    mc.a.acc_global/mc.a.prop_global*100)
                end
            end

            mc.a.acc_rate = 0.0
            mc.a.acc_rate_global = 0.0
            flush(stdout)
            _time = time()
        end

        if safe_before - now() < Millisecond(grace_period) +
                Millisecond(round(Int, 2e3max_sweep_duration))

            println("Early save initiated for sweep #$i.\n")
            verbose && println("Current time: ", Dates.format(now(), "d.u yyyy HH:MM"))
            verbose && println("Target time:  ", Dates.format(safe_before, "d.u yyyy HH:MM"))
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
            verbose && println("\nEarly save finished")

            return false
        elseif (now() - last_checkpoint) > safe_every
            verbose && println("Performing scheduled save.")
            last_checkpoint = now()
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
        end
    end

    mc.a.acc_rate = mc.a.acc_local / mc.a.prop_local
    mc.a.acc_rate_global = mc.a.acc_global / mc.a.prop_global

    if verbose
        if length(mc.a.imaginary_probability) > 0
            s = mc.a.imaginary_probability
            println("\nImaginary Probability Errors: ($(s.count))")
            @printf("\tmax  = %0.3e\n", max(s))
            @printf("\tmean = %0.3e\n", mean(s))
            @printf("\tmin  = %0.3e\n\n", min(s))
        end
        if length(mc.a.negative_probability) > 0
            s = mc.a.negative_probability
            println("\nNegative Probability Errors: ($(s.count))")
            @printf("\tmax  = %0.3e\n", max(s))
            @printf("\tmean = %0.3e\n", mean(s))
            @printf("\tmin  = %0.3e\n\n", min(s))
        end
        if length(mc.a.propagation_error) > 0
            s = mc.a.propagation_error
            println("\nPropagation Errors: ($(s.count))")
            @printf("\tmax  = %0.3e\n", max(s))
            @printf("\tmean = %0.3e\n", mean(s))
            @printf("\tmin  = %0.3e\n\n", min(s))
        end
    end

    end_time = now()
    if verbose
        println("\nEnded: ", Dates.format(end_time, "d.u yyyy HH:MM"))
        @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)
        println()
    end

    return true
end

"""
    update(mc::DQMC, i::Int)

Propagates the Green's function and performs local and global updates at
current imaginary time slice.
"""
function update(mc::DQMC, i::Int)
    propagate(mc)

    # global move
    # if mc.p.global_moves && (current_slice(mc) == mc.p.slices &&
    #        mc.s.direction == -1 && iszero(mod(i, mc.p.global_rate)))
    #     mc.a.prop_global += 1
    #     b = global_move(mc, mc.model, mc.conf) # not yet in DQMC_optional, i.e. unsupported
    #     mc.a.acc_global += b
    # end

    # local moves
    sweep_spatial(mc)

    nothing
end

"""
    sweep_spatial(mc::DQMC)

Performs a sweep of local moves along spatial dimension at current
imaginary time slice.
"""
@bm function sweep_spatial(mc::DQMC)
    m = model(mc)
    N = size(conf(mc), 1)

    @inbounds for i in 1:N
        detratio, ΔE_boson, passthrough = propose_local(mc, m, i, current_slice(mc), conf(mc))
        mc.a.prop_local += 1

        if mc.p.check_sign_problem
            if abs(imag(detratio)) > 1e-6
                push!(mc.a.imaginary_probability, abs(imag(detratio)))
                mc.p.silent || @printf(
                    "Did you expect a sign problem? imag. detratio:  %.9e\n",
                    abs(imag(detratio))
                )
            end
            if real(detratio) < 0.0
                push!(mc.a.negative_probability, real(detratio))
                mc.p.silent || @printf(
                    "Did you expect a sign problem? negative detratio %.9e\n",
                    real(detratio)
                )
            end
        end
        p = real(exp(- ΔE_boson) * detratio)

        # Metropolis
        if p > 1 || rand() < p
            accept_local!(mc, m, i, current_slice(mc), conf(mc), detratio, ΔE_boson, passthrough)
            # Δ, detratio,ΔE_boson)
            mc.a.acc_rate += 1.0
            mc.a.acc_local += 1
        end
    end
    mc.a.acc_rate /= N
    nothing
end

"""
    replay(mc::DQMC[; configurations::Iterable = mc.configs; kwargs...])

Replays previously generated configurations and measures observables along the
way.

### Keyword Arguments (both):
- `verbose = true`: If true, print progress messaged to stdout.
- `safe_every::TimePeriod`: Set the interval for regularly scheduled saves.
- `safe_before::Date`: If this date is passed, `replay!` will generate a
resumable save file and exit
- `grace_period = Minute(5)`: Buffer between the current time and `safe_before`.
The time required to generate a save file should be included here.
- `filename`: Name of the save file. The default is based on `safe_before`.
- `start=1`: The first sweep in the simulation. This will be changed when using
`resume!(save_file)`.
- `ignore`: A collection of measurement keys to ignore. Defaults to the key of
the configuration measurement.
- `measure_rate = 1`: Rate at which measurements are taken. Note that this is
based on the recorded configurations, not actual sweeps.
"""
function replay!(
        mc::DQMC, configurations = mc.configs;
        ignore = tuple(),
        verbose::Bool = true,
        safe_before::TimeType = now() + Year(100),
        safe_every::TimePeriod = Hour(10000),
        grace_period::TimePeriod = Minute(5),
        resumable_filename::String = "resumable_$(Dates.format(safe_before, "d_u_yyyy-HH_MM")).jld",
        overwrite = false,
        measure_rate = 1
    )
    start_time = now()
    last_checkpoint = now()
    max_sweep_duration = 0.0
    verbose && println("Started: ", Dates.format(start_time, "d.u yyyy HH:MM"))

    # Check for measurements
    isempty(mc.thermalization_measurements) && @debug(
        "There is no thermalization process in a replayed simulation."
    )
    isempty(mc.measurements) && @warn(
        "There are no measurements set up for this simulation!"
    )
    # Generate measurement groups
    groups = generate_groups(
        mc, mc.model,
        [mc.measurements[k] for k in keys(mc.measurements) if !(k in ignore)]
    )

    if measure_rate != mc.p.measure_rate
        mc.p = DQMCParameters(
            mc.p.global_moves, mc.p.global_rate,
            mc.p.thermalization, mc.p.sweeps,
            mc.p.silent, mc.p.check_sign_problem,mc.p.check_propagation_error,
            mc.p.safe_mult, mc.p.delta_tau, mc.p.beta, mc.p.slices,
            measure_rate
        )
    end

    verbose && println("Preparing Green's function stack")
    init!(mc)
    build_stack(mc, mc.s)
    propagate(mc)
    mc.s.current_slice = 1
    mc.conf = rand(DQMC, mc.model, nslices(mc))

    _time = time()
    verbose && println("\n\nReplaying measurement stage - ", length(configurations))
    prepare!(mc.measurements, mc, mc.model)
    for i in mc.last_sweep+1:mc.p.measure_rate:length(configurations)
        copyto!(mc.conf, decompress(mc, mc.model, configurations[i]))
        calculate_greens(mc, 0) # outputs to mc.s.greens
        for (requirement, group) in groups
            apply!(requirement, group, mc, mc.model, i)
        end
        mc.last_sweep = i

        if mod(i, 10) == 0
            sweep_dur = (time() - _time)/10
            max_sweep_duration = max(max_sweep_duration, sweep_dur)
            if verbose
                println("\t", i)
                @printf("\t\tsweep dur: %.3fs\n", sweep_dur)
            end
            flush(stdout)
            _time = time()
        end

        if safe_before - now() < Millisecond(grace_period) +
                Millisecond(round(Int, 2e3max_sweep_duration))

            println("Early save initiated for sweep #$i.\n")
            verbose && println("Current time: ", Dates.format(now(), "d.u yyyy HH:MM"))
            verbose && println("Target time:  ", Dates.format(safe_before, "d.u yyyy HH:MM"))
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
            verbose && println("\nEarly save finished")

            return false
        elseif (now() - last_checkpoint) > safe_every
            verbose && println("Performing scheduled save.")
            last_checkpoint = now()
            save(resumable_filename, mc, overwrite = overwrite, rename = false)
        end
    end
    finish!(mc.measurements, mc, mc.model)

    end_time = now()
    verbose && println("Ended: ", Dates.format(end_time, "d.u yyyy HH:MM"))
    verbose && @printf("Duration: %.2f minutes", (end_time - start_time).value/1000. /60.)

    return true
end

"""
    greens(mc::DQMC)

Obtain the current equal-time Green's function, i.e. the fermionic expectation
value of `Gᵢⱼ = ⟨cᵢcⱼ^†⟩`. The indices relate to sites and flavors, but the
exact meanign depends on the model. For the attractive Hubbard model
`G[i, j] = ⟨c_{i, ↑} c_{j, ↑}^†⟩ = ⟨c_{i, ↓} c_{j, ↓}^†⟩` due to symmetry.

Internally, `mc.s.greens` is an effective Green's function. This method
transforms it to the actual Green's function by multiplying hopping matrix
exponentials from left and right.
"""
@bm greens(mc::DQMC) = copy(_greens!(mc))
"""
    greens!(mc::DQMC[; output=mc.s.greens_temp, input=mc.s.greens, temp=mc.s.Ur])
"""
@bm function greens!(mc::DQMC; output=mc.s.greens_temp, input=mc.s.greens, temp=mc.s.Ur)
    _greens!(mc, output, input, temp)
end
"""
    _greens!(mc::DQMC[, output=mc.s.greens_temp, input=mc.s.greens, temp=mc.s.Ur])
"""
function _greens!(
        mc::DQMC_CBFalse, target::AbstractMatrix = mc.s.greens_temp,
        source::AbstractMatrix = mc.s.greens, temp::AbstractMatrix = mc.s.Ur
    )
    eThalfminus = mc.s.hopping_matrix_exp
    eThalfplus = mc.s.hopping_matrix_exp_inv
    vmul!(temp, source, eThalfminus)
    vmul!(target, eThalfplus, temp)
    return target
end
function _greens!(
        mc::DQMC_CBTrue, target::AbstractMatrix = mc.s.greens_temp,
        source::AbstractMatrix = mc.s.greens, temp::AbstractMatrix = mc.s.Ur
    )
    chkr_hop_half_minus = mc.s.chkr_hop_half
    chkr_hop_half_plus = mc.s.chkr_hop_half_inv
    copyto!(target, source)

    @inbounds @views begin
        for i in reverse(1:mc.s.n_groups)
            vmul!(temp, target, chkr_hop_half_minus[i])
            copyto!(target, temp)
        end
        for i in reverse(1:mc.s.n_groups)
            vmul!(temp, chkr_hop_half_plus[i], target)
            copyto!(target, temp)
        end
    end
    return target
end
"""
    greens(mc::DQMC, l::Integer)

Calculates the equal-time greens function at a given slice index `l`, i.e.
`G_{ij}(l, l) = G_{ij}(l⋅Δτ, l⋅Δτ) = ⟨cᵢ(l⋅Δτ)cⱼ(l⋅Δτ)^†⟩`.

Note: This internally overwrites the stack variables `Ul`, `Dl`, `Tl`, `Ur`,
`Dr`, `Tr`, `curr_U`, `tmp1` and `tmp2`. All of those can be used as temporary
or output variables here, however keep in mind that other results may be
invalidated. (I.e. `G = greens!(mc)` would get overwritten.)
"""
@bm greens(mc::DQMC, slice::Integer) = copy(_greens!(mc, slice))
"""
    greens!(mc::DQMC, l::Integer[; output=mc.s.greens_temp, temp1=mc.s.tmp1, temp2=mc.s.tmp2])
"""
@bm function greens!(
        mc::DQMC, slice::Integer;
        output=mc.s.greens_temp, temp1 = mc.s.tmp1, temp2 = mc.s.tmp2
    )
    _greens!(mc, slice, output, temp1, temp2)
end
"""
    _greens!(mc::DQMC, l::Integer[, output=mc.s.greens_temp, temp1=mc.s.tmp1, temp2=mc.s.tmp2])
"""
function _greens!(
        mc::DQMC, slice::Integer, output::AbstractMatrix = mc.s.greens_temp,
        temp1::AbstractMatrix = mc.s.tmp1, temp2::AbstractMatrix = mc.s.tmp2
    )
    calculate_greens(mc, slice, temp1)
    _greens!(mc, output, temp1, temp2)
end

################################################################################
### FileIO
################################################################################

#     save_mc(filename, mc, entryname)
#
# Saves (minimal) information necessary to reconstruct a given `mc::DQMC` to a
# JLD-file `filename` under group `entryname`.
#
# When saving a simulation the default `entryname` is `MC`
function save_mc(file::JLDFile, mc::DQMC, entryname::String="MC")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(mc))
    save_parameters(file, mc.p, entryname * "/Parameters")
    save_analysis(file, mc.a, entryname * "/Analysis")
    write(file, entryname * "/conf", mc.conf)
    _save(file, mc.configs, entryname * "/configs")
    write(file, entryname * "/last_sweep", mc.last_sweep)
    save_measurements(file, mc, entryname * "/Measurements")
    save_model(file, mc.model, entryname * "/Model")
    nothing
end

CB_type(T::UnionAll) = T.body.parameters[2]
CB_type(T::DataType) = T.parameters[2]

#     load_mc(data, ::Type{<: DQMC})
#
# Loads a DQMC from a given `data` dictionary produced by `JLD.load(filename)`.
function _load(data, ::Type{T}) where T <: DQMC
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $T version $(data["VERSION"])"))
    end

    CB = CB_type(data["type"])
    @assert CB <: Checkerboard
    mc = DQMC(CB)
    mc.p = _load(data["Parameters"], data["Parameters"]["type"])
    mc.a = _load(data["Analysis"], data["Analysis"]["type"])
    mc.conf = data["conf"]
    mc.configs = _load(data["configs"], data["configs"]["type"])
    mc.last_sweep = data["last_sweep"]
    mc.model = _load(data["Model"], data["Model"]["type"])

    measurements = _load(data["Measurements"], Measurements)
    mc.thermalization_measurements = measurements[:TH]
    mc.measurements = measurements[:ME]
    HET = hoppingeltype(DQMC, mc.model)
    GET = greenseltype(DQMC, mc.model)
    HMT = hopping_matrix_type(DQMC, mc.model)
    GMT = greens_matrix_type(DQMC, mc.model)
    IMT = interaction_matrix_type(DQMC, mc.model)
    mc.s = DQMCStack{GET, HET, GMT, HMT, IMT}()
    mc.ut_stack = UnequalTimeStack{GET, GMT}()

    make_concrete!(mc)
end

#   save_parameters(file::JLDFile, p::DQMCParameters, entryname="Parameters")
#
# Saves (minimal) information necessary to reconstruct a given
# `p::DQMCParameters` to a JLD-file `filename` under group `entryname`.
#
# When saving a simulation the default `entryname` is `MC/Parameters`
function save_parameters(file::JLDFile, p::DQMCParameters, entryname::String="Parameters")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(p))

    write(file, entryname * "/global_moves", Int(p.global_moves))
    write(file, entryname * "/global_rate", p.global_rate)
    write(file, entryname * "/thermalization", p.thermalization)
    write(file, entryname * "/sweeps", p.sweeps)
    write(file, entryname * "/silent", Int(p.silent))
    write(file, entryname * "/check_sign_problem", Int(p.check_sign_problem))
    write(file, entryname * "/check_propagation_error", Int(p.check_propagation_error))
    write(file, entryname * "/safe_mult", p.safe_mult)
    write(file, entryname * "/delta_tau", p.delta_tau)
    write(file, entryname * "/beta", p.beta)
    write(file, entryname * "/slices", p.slices)
    write(file, entryname * "/measure_rate", p.measure_rate)

    nothing
end

#     load_parameters(data, ::Type{<: DQMCParameters})
#
# Loads a DQMCParameters object from a given `data` dictionary produced by
# `JLD.load(filename)`.
function _load(data, ::Type{T}) where T <: DQMCParameters
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $T version $(data["VERSION"])"))
    end

    data["type"](
        Bool(data["global_moves"]),
        data["global_rate"],
        data["thermalization"],
        data["sweeps"],
        Bool(data["silent"]),
        Bool(data["check_sign_problem"]),
        Bool(data["check_propagation_error"]),
        data["safe_mult"],
        data["delta_tau"],
        data["beta"],
        data["slices"],
        data["measure_rate"],
    )
end

function save_analysis(file::JLDFile, a::DQMCAnalysis, entryname::String="Analysis")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/type", typeof(a))

    save_stats(file, a.imaginary_probability, entryname * "/imag_prob")
    save_stats(file, a.negative_probability, entryname * "/neg_prob")
    save_stats(file, a.propagation_error, entryname * "/propagation")
end
function save_stats(file::JLDFile, ms::MagnitudeStats, entryname::String="MStats")
    write(file, entryname * "/max", ms.max)
    write(file, entryname * "/min", ms.min)
    write(file, entryname * "/sum", ms.sum)
    write(file, entryname * "/count", ms.count)
end

function _load(data, ::Type{T}) where T <: DQMCAnalysis
    if !(data["VERSION"] == 1)
        throw(ErrorException("Failed to load $T version $(data["VERSION"])"))
    end

    data["type"](
        imaginary_probability = load_stats(data["imag_prob"]),
        negative_probability = load_stats(data["neg_prob"]),
        propagation_error = load_stats(data["propagation"])
    )
end
function load_stats(data)
    MagnitudeStats(data["max"], data["min"], data["sum"], data["count"])
end

include("DQMC_mandatory.jl")
include("DQMC_optional.jl")
include("measurements/generic.jl")


text:
using Test, KrigingModel
using Statistics
using GaussianProcesses

include("scaler.jl")
include("optimizer.jl")
include("kriging.jl")