text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Unit test file for hdrimage.jl

c1 = RGB( 1f0,  2f0,  3f0)
c2 = RGB( 4f0,  5f0,  6f0)
c3 = RGB( 7f0,  8f0,  9f0)
c4 = RGB(10f0, 11f0, 12f0)
c5 = RGB(13f0, 14f0, 15f0)
c6 = RGB(16f0, 17f0, 18f0)
c7 = RGB(19f0, 20f0, 21f0)
c8 = RGB(22f0, 23f0, 24f0)

@testset "constructors" begin
    # test constructor from matrix
    @testset "from matrix" begin
        rgb_pixel_matrix = [c1 c3 c5
                            c2 c4 c6]
        image = HdrImage(rgb_pixel_matrix)
        @test all(image.pixel_matrix .=== rgb_pixel_matrix)
    end

    # test custom constructors
    @testset "size constructor given type" begin
        img_width, img_height = 3, 2
        rgb_zeros = zeros(RGB{Float32}, img_width, img_height)

        image = HdrImage(img_width, img_height)
        @test image.pixel_matrix == rgb_zeros
        @test all(image.pixel_matrix .=== rgb_zeros)

        # convenience alias
        image = HdrImage(img_width, img_height)
        @test image.pixel_matrix == rgb_zeros
        @test all(image.pixel_matrix .=== rgb_zeros)
    end

    @testset "size constructor default type" begin
        img_width, img_height = 3, 2
        rgb_zeros = zeros(RGB{Float32}, img_width, img_height)
        image = HdrImage(img_width, img_height)
        @test image.pixel_matrix == rgb_zeros
        @test all(image.pixel_matrix .=== rgb_zeros)

        # test correct defaulting
        @test image.pixel_matrix == HdrImage(img_width, img_height).pixel_matrix
        @test all(image.pixel_matrix .=== HdrImage(img_width, img_height).pixel_matrix)
    end

    @testset "array and size" begin
        img_width, img_height = 3, 2
        arr = collect(RGB(1f0, 2f0, 3f0) .+ 3f0i for i ∈ 0:img_width*img_height-1)
        image = HdrImage(arr, img_width, img_height)
        @test image.pixel_matrix == reshape(arr, img_width, img_height)
        @test all(image.pixel_matrix .=== reshape(arr, img_width, img_height))
    end

    @testset "array and shape" begin
        img_width, img_height = 3, 2
        shape = (img_width, img_height)
        arr = collect(RGB(1f0, 2f0, 3f0) .+ 3f0i for i ∈ 0:img_width*img_height-1)
        image = HdrImage(arr, shape)
        @test image.pixel_matrix == reshape(arr, shape)
        @test all(image.pixel_matrix .=== reshape(arr, shape))
    end
end

@testset "iterations" begin
    # testset variables
    rgb_pixel_matrix = [c1 c3 c5
                        c2 c4 c6]
    image = HdrImage(rgb_pixel_matrix)

    # test indexing properties
    @testset "indexing properties" begin
        @test length(image) === 6
        @test firstindex(image) === 1
        @test lastindex(image) === 6
    end

    # test get value
    @testset "get value" begin
        # linear indexing
        @test c1 == image[begin]
        @test c1 == image[1]
        @test c2 == image[2]
        @test c3 == image[3]
        @test c4 == image[4]
        @test c5 == image[5]
        @test c6 == image[end]
        @test c6 == image[6]

        # test exceptions
        @test_throws BoundsError image[7]

        # cartesian indexing
        @test c1 == image[1,1]
        @test c2 == image[2,1]
        @test c3 == image[1,2]
        @test c4 == image[2,2]
        @test c5 == image[1,3]
        @test c6 == image[2,3]

        @test c3 == image[CartesianIndex(1,2)]

        # test exceptions
        @test_throws BoundsError image[0,1]
        @test_throws BoundsError image[1,0]
        @test_throws BoundsError image[CartesianIndex(6,7)]
    end

    # test set value
    @testset "set value" begin
        _image = HdrImage(1, 2)
        _image[1] = c1
        _image[2] = c2
        @test _image.pixel_matrix == [c1 c2]
    end

    # test iterability
    @testset "iterability" begin
        @test all(i == j for (i, j) in zip((c1, c2, c3, c4, c5, c6), image))
    end
end

@testset "Broadcasting" begin
    a = 2
    rgb_pixel_matrix_1 = [c1 c3
                          c2 c4]
    img_1 = HdrImage(rgb_pixel_matrix_1)
    rgb_pixel_matrix_2 = [c5 c7
                          c6 c8]
    img_2 = HdrImage(rgb_pixel_matrix_2)

    # testing equivalence to custom defined methods
    @test all(img_1 .+ img_2 .== img_1.pixel_matrix + img_2.pixel_matrix)
    @test all(img_1 .- img_2 .== img_1.pixel_matrix - img_2.pixel_matrix)
    @test all(img_1 .* img_2 .== img_1.pixel_matrix .* img_2.pixel_matrix)
    @test all(a .* img_1 .== a * img_1.pixel_matrix)

    # broadcasting operators can be applied between any broadcastable type instances
    img_3 = HdrImage([c1 c1])
    @test all(((Ref(RGB(.1, .2, .3)) .+ img_3) .≈ Ref(RGB(1.1, 2.2, 3.3))))
    @test all(Ref(c1) .== img_3)

    # # it works for any operator valid for the types of the elements
    @test all((img_1 .* img_2) .≈ HdrImage([c1*c5 c3*c7; c2*c6 c4*c8]))
end

@testset "Tone mapping" begin
    @testset "normalize" begin
        img = HdrImage([RGB(  5f0,   10f0,   15f0)
                        RGB(500f0, 1000f0, 1500f0)], 2, 1)

        img = normalize(img, 1000f0, luminosity=100f0)
        @test all(img .≈ [RGB(0.5f2, 1.0f2, 1.5f2), RGB(0.5f4, 1.0f4, 1.5f4)])
    end

    @testset "clamp" begin
        img = HdrImage([RGB(0.5f1, 1.0f1, 1.5f1)
                        RGB(0.5f3, 1.0f3, 1.5f3)], 2, 1)
        img = clamp(img)
        @test all(0 <= col <= 1 for pix ∈ img for col ∈ pix)
    end

    @testset "luminosity" begin
        image = HdrImage([RGB(  5f0,   10f0,   15f0) RGB(500f0, 1000f0, 1500f0)])
        @test luminosity(image) ≈ 100f0
    end

end

@testset "IO" begin
    # test color pretty printing
    @testset "show" begin
        io = IOBuffer()
        rgb_pixel_matrix = [c1 c4
                            c2 c5
                            c3 c6]
        image = HdrImage(rgb_pixel_matrix)
        # compact
        show(io, image)
        @test String(take!(io)) == " (1.0 2.0 3.0)  (10.0 11.0 12.0)\n (4.0 5.0 6.0)  (13.0 14.0 15.0)\n (7.0 8.0 9.0)  (16.0 17.0 18.0)"
        # extended
        show(io, "text/plain", image)
        @test String(take!(io)) == "3x2 HdrImage:\n (1.0 2.0 3.0)  (10.0 11.0 12.0)\n (4.0 5.0 6.0)  (13.0 14.0 15.0)\n (7.0 8.0 9.0)  (16.0 17.0 18.0)"
    end
end

@testset "Other" begin
    # testset variables
    img_width, img_height = 2, 3
    rgb_pixel_matrix = [c1 c3 c5
                        c2 c4 c6]
    image = HdrImage(rgb_pixel_matrix)

    # test eltype
    @testset "eltype" begin
        @test eltype(image) == RGB{Float32}
    end

    # test fill!
    @testset "fill!" begin
        a = HdrImage(fill(RGB(NaN32,NaN32,NaN32), img_width, img_height))
        @test all(HdrImage(img_width, img_height) .== fill!(a, zero(RGB{Float32})))

        a = HdrImage(ones(RGB{Float32}, img_width, img_height))
        @test all(fill!(HdrImage(img_width, img_height), one(RGB{Float32})) .== a)
    end

    # test size
    @testset "size" begin
        @test size(image) == (img_width, img_height)
    end
end


text:
# Tests for distances
function t_dist()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

# Tests
D = [Distances.Euclidean(),
     Distances.Cityblock(),
     Distances.Chebyshev(),
     Distances.Jaccard()] # Add more just to be sure ...
for tr in [A01, A02, A03, vectordata]
	for ts in [A01, A02, A03, vectordata]
		for d in D
			Test.@test try
				ts |> (tr |> j4pr.dist(d));
				true
			catch
				false
			end
		end
	end
end

for tr in [A04, A05, A06, matrixdata]
	for ts in [A04, A05, A06, matrixdata]
		for d in D
			Test.@test try
				ts |> (tr |> j4pr.dist(d));
				true
			catch
				false
			end
		end
	end
end

# Test training with tuple and the generic dist(Array,Array)
Test.@test try
	w = j4pr.dist(j4pr.strip(A03),D[1])
	true
catch
	false
end

end

ext: jl


text:
struct GLEScript
    io::IOBuffer
end
GLEScript() = GLEScript(IOBuffer())

const GS = GLEScript

# write to buffer (or buffer encapsulating object) with form (s |> b)
|>(s::String, g::GS)      = write(g.io, s, " ")  # writing to g
|>(s::String, tio::Tuple) = @. |>(s, tio)        # writing to (g1, g2, ...)

# read buffer encapsulated by `b`
take!(g::GS) = take!(g.io)

|>(gi::GS, go::GS) = write(go.io, take!(gi))

gsave(g::GS)    = "\ngsave"    |> g
grestore(g::GS) = "\ngrestore" |> g


text:
using ITensors
include("dmrg_lbo.jl")

## STRUCTS ##

struct Parameters
    # Model parameters
    N::Int
    Nx::Int
    Ny::Int
    t::Real
    U::Real
    ω::Real
    g0::Real
    g1::Real
    λ::Real
    yperiodic::Bool
    doping::Real
    max_phonons::Int
    init_phonons::Int

    # DMRG parameters
    DMRG_numsweeps
    DMRG_noise
    DMRG_maxdim
    DMRG_cutoff
    DMRG_LBO::Bool
    max_LBO_dim::Int
    min_LBO_dim::Int

    # TEBD parameters
    mid::Int
    T::Int
    τ::Real
    TEBD_cutoff
    TEBD_maxdim
    TEBD_LBO::Bool
end

struct HubbardHolsteinModel
    # Lattice
    sites
    # Hamiltonian (in two forms)
    mpo::MPO
    gates
end

struct DMRGResults
    ground_state
    ground_state_energy
    ground_state_entropy
    optimized_basis
end

struct EquilibriumCorrelations
    spin
    charge
    sSC
    pSC
    dSC
end

struct TEBDResults
    entropy
    self_overlap
    corrs
    phi_t
    psi_t
end

## SETTING UP THE MODEL ##

function parameters(;Nx::Int, Ny::Int, t::Real, U::Real=nothing, ω::Real=nothing,
                    g0::Real=nothing, g1::Real=nothing, λ::Real=nothing, yperiodic=false,
                    doping::Real=0, max_phonons::Int=1, init_phonons::Int=0,
                    DMRG_numsweeps::Int=20, DMRG_noise=nothing,
                    DMRG_maxdim=nothing, DMRG_cutoff=nothing, DMRG_LBO=false,
                    max_LBO_dim=nothing, min_LBO_dim=4,
                    T::Int=25, τ::Real=0.1, TEBD_cutoff=1E-14, TEBD_maxdim=400, TEBD_LBO=false, )

    if isnothing(U)
        U = 8*t
    end
    if isnothing(ω)
        ω = 0.5*t
    end
    if isnothing(g0)
        g0 = 0.05*t
    end
    if isnothing(g1)
        g1 = 0
    end
    if isnothing(λ)
        λ = 0
    end
    if max_phonons==0
        @assert ω==g0==g1==0
    end
    if ω==0
        @assert λ==0
    end

    if isnothing(DMRG_noise)
        DMRG_noise = [1E-6,1E-6,1E-8,0]
    end
    if isnothing(DMRG_maxdim)
        DMRG_maxdim = [20,40,100,200,400]
    end
    if isnothing(DMRG_cutoff)
        DMRG_cutoff = 1E-10
    end
    if isnothing(max_LBO_dim)
        max_LBO_dim = 12
    end

    Parameters(Nx*Ny,Nx,Ny,t,U,ω,g0,g1,λ,yperiodic,doping,max_phonons,init_phonons,
                DMRG_numsweeps,DMRG_noise,DMRG_maxdim,DMRG_cutoff,DMRG_LBO,
                max_LBO_dim,min_LBO_dim,
                ceil(Int,Nx*Ny/2),T,τ,TEBD_cutoff,TEBD_maxdim,TEBD_LBO)
end

"""
Use this function for reconstructing from scratch given site indices
"""
function HubbardHolsteinModel(p::Parameters, sites::Vector{Index{Vector{Pair{QN, Int64}}}})
    N, Nx, Ny, t, U, ω, g0, g1, λ, τ, max_phonons, yperiodic = p.N, p.Nx, p.Ny, p.t, p.U, p.ω, p.g0, p.g1, p.λ, p.τ, p.max_phonons, p.yperiodic

    lattice = square_lattice(Nx, Ny; yperiodic=yperiodic)

    # make the hamiltonian
    ampo = OpSum()

    # bond terms
    for b in lattice
        ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
        ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
        ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
        ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
    end

    # on-site terms
    for j=1:N-1
        # ∑_j U * n_j↑ n_j↓
        ampo += U,"Nupdn",j,"I",j

        # ∑_j ω * nb_j
        if max_phonons >= 1
            ampo += ω,"Nb",j

            # # ∑_j g0 * nf_j (b^†_j + b_j)
            ampo += g0,"Ntot(Bd+B)",j

            # # ∑_⟨ij⟩ g1 * nf_j (b^†_i + b_i)
            ampo += g1,"Ntot",j,"Bdag+B",j+1
            ampo += g1,"Ntot",j+1,"Bdag+B",j

            # quartic term
            ampo += λ,"Nb^2",j
        end
    end
    # Edge site
    ampo += U,"Nupdn",N
    if max_phonons >= 1
        ampo += ω,"Nb",N
        ampo += g0,"Ntot(Bd+B)",N
        ampo += λ,"Nb^2",N
    end
    H = MPO(ampo,sites)

    ## Helper functions for making the Trotter gates ##
    function make_twosite(s1, s2)
        hj_twosite = -t*(op("Cdagup*F",s1) * op("Cup",s2)  # t * (c^†_jσ c_{j+1}σ + h.c.)
                 -op("Cup*F",s1) * op("Cdagup",s2)
                 +op("Cdagdn*F",s1) * op("Cdn",s2)
                 -op("Cdn*F",s1) * op("Cdagdn",s2))
        if max_phonons >= 1
            hj_twosite  = hj_twosite
                + g1*(op("Ntot",s1) * op("Bdag+B",s2))
                + g1*(op("Bdag+B",s1) * op("Ntot",s2))
        end

        return exp(-1.0im * τ/2 * hj_twosite)
    end

    function make_onesite(s1)
        hj_onesite = U * op("Nupdn",s1)
        if max_phonons >=1
            hj_onesite = hj_onesite
                    + ω * op("Nb",s1)
                    + g0 * op("Ntot(Bd+B)",s1)
                    + λ * op("Nb^2",s1)
        end

        return exp(-1.0im * τ/2 * hj_onesite)
    end

    function make_endsite(sn)
        hn = U*op("Nupdn",sn)
        if max_phonons >= 1
            hn = hn
            + ω*op("Nb",sn)
            + g0*op("Ntot(Bd+B)",sn)
            + λ*op("Nb^2",sn)
        end
        return exp(-1.0im * τ/2 * hn)
    end

    # make the trotter gates
    gates = ITensor[]
    for b in lattice
        s1_idx = min(b.s1, b.s2)
        s2_idx = max(b.s1, b.s2)
        s1 = sites[s1_idx]
        s2 = sites[s2_idx]

        # DEBUGGING !! ##
        # Check if we need to apply any swaps (sites are not adjacent)
        if abs(s2_idx-s1_idx) > 1
            # Apply SWAP to get to new configuration
            SWAP_forward = []
            s_last = s1_idx
            for s_next in s1_idx+1:(s2_idx-1)
                push!(SWAP_forward, op("fSWAP", sites[s_last], sites[s_next]))
                s_last = s_next
            end
            append!(gates, SWAP_forward)

            # Apply gate
            hj_twosite = make_twosite(sites[s2_idx-1], sites[s2_idx])
            push!(gates,hj_twosite)

            # Apply SWAP to get back
            SWAP_backward = []
            for s_next in s_last-1:s1_idx
                push!(SWAP_backward, op("fSWAP", sites[s_last], sites[s_next]))
                s_last = s_next
            end
            append!(gates, SWAP_backward)
            #append!(gates, reverse(SWAP_forward))

        else # no need to apply swaps
            hj_twosite = make_twosite(s1, s2)
            push!(gates,hj_twosite)
        end

    end

    # One-site terms (no swap needed)
    for j=1:N-1
        hj_onesite = make_onesite(sites[j])
        push!(gates,hj_onesite)
    end

    # End site (no swap needed)
    hn = make_endsite(sites[N])
    push!(gates,hn)

    # Add all the gates in reverse
    append!(gates,reverse(gates))

    # Return the struct
    HubbardHolsteinModel(sites, H, gates)
end

function HubbardHolsteinModel(p::Parameters)
    # make the sites
    sites = siteinds("HubHolst",p. N; dim=p.max_phonons+1)

    return HubbardHolsteinModel(p, sites)
end

## GENERIC STATE FUNCTIONS ##

function compute_overlap(ψ1::MPS, ψ2::MPS)
    LinearAlgebra.norm(inner(ψ1, ψ2))
end

function compute_phonon_number(ψ::MPS)
    expect(ψ,"Nb")
end

function compute_electron_number(ψ::MPS)
    expect(ψ,"Ntot")
end

function compute_entropy(ψ::MPS, b::Int)
    orthogonalize!(ψ, b)
    U,S,V = svd(ψ[b], (linkind(ψ, b-1), siteind(ψ,b)))
    SvN = 0.0
    for n=1:dim(S, 1)
        p = S[n,n]^2
        SvN -= p * log(p)
    end
    return SvN
end

## DMRG ##

function initialize_wavefcn(HH::HubbardHolsteinModel, p::Parameters)
    # Extract the information
    ups = "Up,"*string(p.init_phonons)
    downs = "Dn,"*string(p.init_phonons)
    emps = "Emp,"*string(p.init_phonons)

    state = [isodd(n) ? ups : downs for n=1:p.N]

    # Account for doping
    if p.doping > 0
        spacing = floor(Int,1/p.doping)
        state[1:spacing:end] .= emps
    end

    # NOTE: the QN of this state is preserved during DMRG
    productMPS(HH.sites,state)
end

function run_DMRG(HH::HubbardHolsteinModel, p::Parameters; alg="divide_and_conquer")
    # Set DMRG params
    sweeps = Sweeps(p.DMRG_numsweeps)
    setnoise!(sweeps, p.DMRG_noise...) # Very important to use noise for this model
    setmaxdim!(sweeps, p.DMRG_maxdim...)
    setcutoff!(sweeps, p.DMRG_cutoff...)

    ϕ0 = initialize_wavefcn(HH,p)
    @show flux(ϕ0)
    if p.DMRG_LBO # If performing local basis optimization
        energy, ϕ, Rs = dmrg_lbo(HH.mpo, ϕ0, sweeps, alg=alg, LBO=true,
                                    max_LBO_dim=p.max_LBO_dim, min_LBO_dim=p.min_LBO_dim)
    else
        energy, ϕ = dmrg(HH.mpo, ϕ0, sweeps, alg=alg)
        Rs = nothing
    end
    entropy = compute_entropy(ϕ, p.mid)
    return DMRGResults(ϕ, energy, entropy, Rs)
end

## CORRELATION FUNCTIONS ##

function apply_onesite_operator(ϕ::MPS, opname::String, sites, siteidx::Int)
    ϕ = copy(ϕ)

    ## Account for fermion sign using Jordan-Wigner strings ##
    if opname == "Cup" || opname == "Cdn"
        ϕ = apply_op(ϕ, op(opname,sites[siteidx]), siteidx)
        for i in reverse(1:(siteidx-1)) # Don't act with string on-site
            ϕ = apply_op(ϕ, op("F",sites[i]), i)
        end
        return ϕ

    elseif opname == "Cdagup" || opname == "Cdagdn"
        for i in 1:(siteidx-1) # Don't act with string on-site
            ϕ = apply_op(ϕ, op("F",sites[i]), i)
        end
        ϕ = apply_op(ϕ, op(opname,sites[siteidx]), siteidx)
        return ϕ
    end

    # Otherwise, just apply the operator as usual
    return apply_op(ϕ, op(opname,sites[siteidx]), siteidx)
end

function apply_twosite_operator(ϕ::MPS, opname::String, sites, siteidx::Int)
    ϕ = copy(ϕ)
    if opname == "pSC"
        # Apply F string first
        ϕ = apply_op(ϕ, op("F",sites[siteidx]), siteidx)

        # Then apply the operator
        Π = 1/sqrt(2)*(op("Cup",sites[siteidx])*op("Cdn",sites[siteidx+1])
                            + op("Cdn",sites[siteidx])*op("Cup",sites[siteidx+1]))
        return apply_op_twosite(ϕ,Π,siteidx)

    elseif opname == "pSCdag"
        # Apply the operator first
        Π_dag = 1/sqrt(2)*(op("Cdagdn",sites[siteidx+1])*op("Cdagup",sites[siteidx])
                            + op("Cdagup",sites[siteidx+1])*op("Cdagdn",sites[siteidx]))
        ϕ = apply_op_twosite(ϕ,Π_dag,siteidx)

        # Apply the F string
        return apply_op(ϕ, op("F",sites[siteidx]), siteidx)

    elseif opname == "dSC"
        # Apply F string first
        ϕ = apply_op(ϕ, op("F",sites[siteidx]), siteidx)

        # Then apply the operator
        Δ = 1/sqrt(2)*(op("Cup",sites[siteidx])*op("Cdn",sites[siteidx+1])
                            - op("Cdn",sites[siteidx])*op("Cup",sites[siteidx+1]))
        return apply_op_twosite(ϕ,Δ,siteidx)

    elseif opname == "dSCdag"
        # Apply the operator first
        Δ_dag = 1/sqrt(2)*(op("Cdagdn",sites[siteidx+1])*op("Cdagup",sites[siteidx])
                            - op("Cdagup",sites[siteidx+1])*op("Cdagdn",sites[siteidx]))
        ϕ = apply_op_twosite(ϕ,Δ_dag,siteidx)

        # Apply the F string
        return apply_op(ϕ, op("F",sites[siteidx]), siteidx)
    end
    @error "No recognized two-site operator"
end

function apply_op(ϕ::MPS, op::ITensor, siteidx::Int)
    ϕ = copy(ϕ) # Make a copy of the original state
    orthogonalize!(ϕ, siteidx)
    new_ϕj = op * ϕ[siteidx] # Apply the local operator
    noprime!(new_ϕj)
    ϕ[siteidx] = new_ϕj
    return ϕ
end

function apply_op_twosite(ϕ::MPS, G::ITensor, siteidx::Int; cutoff=1E-8)
    ϕ = copy(ϕ)
    # Note: siteidx corresponds to the leftermost site
    orthogonalize!(ϕ,siteidx)
    wf = (ϕ[siteidx] * ϕ[siteidx+1]) * G
    noprime!(wf)

    inds_site = uniqueinds(ϕ[siteidx],ϕ[siteidx+1])
    U,S,V = svd(wf,inds_site,cutoff=cutoff)
    ϕ[siteidx] = U
    ϕ[siteidx+1] = S*V
    return ϕ
end

function compute_all_equilibrium_correlations(dmrg_results::DMRGResults,
                                            HH::HubbardHolsteinModel;
                                            start=nothing, stop=nothing)

    N = length(HH.sites)
    if isnothing(start)
        start = floor(Int,0.25*N)
    end
    if isnothing(stop)
        stop = ceil(Int,0.75*N)
    end

    corrtypes = ["spin","charge","sSC","pSC","dSC"]
    corrs = []
    for corrtype in corrtypes
        corr = equilibrium_correlations(dmrg_results,corrtype,HH,start,stop)
        push!(corrs,corr)
    end
    EquilibriumCorrelations(corrs...)
end

function equilibrium_correlations(dmrg_results::DMRGResults, corrtype::String,
                                HH::HubbardHolsteinModel, start::Int, stop::Int,)

    ϕ = copy(dmrg_results.ground_state)
    j = start
    sites = HH.sites
    if corrtype=="spin"
        return correlation_matrix(ϕ, "Sz", "Sz")[start:stop,j]
    elseif corrtype=="charge"
        ninj = correlation_matrix(ϕ, "Ntot", "Ntot")[start:stop,j]
        ni = expect(ϕ, "Ntot")
        nj = ni[j]
        return ninj - nj .* (ni[start:stop])
    elseif corrtype=="sSC"
        ψ = apply_onesite_operator(ϕ, "Cupdn", sites, j)
        function compute_corr_sSC(i::Int)
            Σ_iψ = apply_onesite_operator(ψ, "Cdagupdn", sites, i)
            return inner(ϕ,Σ_iψ)
        end
        return compute_corr_sSC.(collect(start:stop))
    elseif corrtype=="pSC"
        ψ = apply_twosite_operator(ϕ, "pSC", sites, j)
        function compute_corr_pSC(i::Int)
            Π_iψ = apply_twosite_operator(ψ, "pSCdag", sites, i)
            return inner(ϕ,Π_iψ)
        end
        return compute_corr_pSC.(collect(start:stop))
    elseif corrtype=="dSC"
        ψ = apply_twosite_operator(ϕ, "dSC", sites, j)
        function compute_corr_dSC(i::Int)
            Δ_iψ = apply_twosite_operator(ψ, "dSCdag", sites, i)
            return inner(ϕ,Δ_iψ)
        end
        return compute_corr_dSC.(collect(start:stop))
    end
    @error "Unknown spin correlation type"
end

function compute_correlations(dmrg_results::DMRGResults,
    A_t0::String, A_t::String,
    HH::HubbardHolsteinModel, p::Parameters;
    interim_save=false, savepath=nothing)

    if p.TEBD_LBO
        return tebd_lbo_corrs(dmrg_results,
        A_t0, A_t,
        HH, p;
        interim_save=interim_save, savepath=savepath)
    else
        return tebd_corrs(dmrg_results,
        A_t0, A_t,
        HH, p;
        interim_save=interim_save, savepath=savepath)
    end
end

function tebd_lbo_corrs(dmrg_results::DMRGResults,
                        A_t0::String, A_t::String,
                        HH::HubbardHolsteinModel, p::Parameters;
                        interim_save=false, savepath=nothing)

    # The wavefunction being acted upon at t=0, |ψ⟩ = A_t0|ϕ⟩
    ϕ = copy(dmrg_results.ground_state)
    ψ = copy(ϕ)
    Rs = copy(dmrg_results.optimized_basis)
    gates = copy(HH.gates)

    # Apply A_t0 to middle site
    ψ = apply_onesite_operator(ψ, A_t0, HH.sites, p.mid)

    # Measure correlation before time evolution
    A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, p.mid)
    # Measure the correlation fcn
    @show inner(ϕ, A_tψ)

    # Parameters for time evolution
    nsteps = floor(p.T/p.τ) # Number of time steps for time evolution
    t = 0.0

    # Results
    corrs = []
    entropy = []
    self_overlap = Float64[]

    for step in 1:nsteps
        print(floor(Int,step),"-")

        ## TO DO: IMPLEMENT LBO FOR THE TEBD STEP ##
        # 1. Take the gate and contract (on only ONE side) with the optimized basis
        # 2. Put ψ in the optimized basis
        # 3. Contract ψ and gates in the optimized basis. Will be in the bare basis on the uncontracted side
        # 4. Optimize the basis of ψ again using the local reduced density matrix
        # 5. Update the rotation matrix

        ϕ = apply(HH.gates, ϕ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff)
        ψ = apply(HH.gates, ψ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff)

        t += p.τ

        ### SANITY CHECKS
        if step%1==0
            # Compute entropy
            push!(entropy, [compute_entropy(ϕ, p.mid),compute_entropy(ψ, p.mid)])
            println("Entropy of ϕ, ψ: ", entropy[end])

            # Compute overlap of ϕ with its original self
            push!(self_overlap, compute_overlap(ϕ,dmrg_results.ground_state))
            println("Overlap of ϕ with itself: ", self_overlap[end])
        end

        # Calculate ⟨ϕ(t)|c_j^† c_i|ϕ(0)⟩
        function measure_corr(j::Int)
            # Apply the second measurement operator
            A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, j)
            return dot(ϕ,A_tψ)
        end

        # Measure the correlation fcn
        push!(corrs,measure_corr.(collect(1:p.N)))

        if interim_save && step%10==0
            @assert !isnothing(savepath)
            tebd_results_interim = TEBDResults(hcat(entropy...), self_overlap,
                                hcat(corrs...), ϕ, ψ)
            save_structs(tebd_results_interim, save_path)
        end

    end

    return TEBDResults(hcat(entropy...), self_overlap,
                        hcat(corrs...), ϕ, ψ)

end

function tebd_corrs(dmrg_results::DMRGResults,
                            A_t0::String, A_t::String,
                            HH::HubbardHolsteinModel, p::Parameters;
                            interim_save=false, savepath=nothing)

    # The wavefunction being acted upon at t=0, |ψ⟩ = A_t0|ϕ⟩
    ϕ = copy(dmrg_results.ground_state)
    ψ = copy(ϕ)

    # Apply A_t0 to middle site
    ψ = apply_onesite_operator(ψ, A_t0, HH.sites, p.mid)

    # Measure correlation before time evolution
    A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, p.mid)
    # Measure the correlation fcn
    @show inner(ϕ, A_tψ)

    # Parameters for time evolution
    nsteps = floor(p.T/p.τ) # Number of time steps for time evolution
    t = 0.0

    # Results
    corrs = []
    entropy = []
    self_overlap = Float64[]

    for step in 1:nsteps
        print(floor(Int,step),"-")

        ϕ = apply(HH.gates, ϕ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff)
        ψ = apply(HH.gates, ψ; maxdim=p.TEBD_maxdim, cutoff=p.TEBD_cutoff)

        t += p.τ

        ### SANITY CHECKS
        if step%1==0
            # Compute entropy
            push!(entropy, [compute_entropy(ϕ, p.mid),compute_entropy(ψ, p.mid)])
            println("Entropy of ϕ, ψ: ", entropy[end])

            # Compute overlap of ϕ with its original self
            push!(self_overlap, compute_overlap(ϕ,dmrg_results.ground_state))
            println("Overlap of ϕ with itself: ", self_overlap[end])
        end

        # Calculate ⟨ϕ(t)|c_j^† c_i|ϕ(0)⟩
        function measure_corr(j::Int)
            # Apply the second measurement operator
            A_tψ = apply_onesite_operator(ψ, A_t, HH.sites, j)
            return dot(ϕ,A_tψ)
        end

        # Measure the correlation fcn
        push!(corrs,measure_corr.(collect(1:p.N)))

        if interim_save && step%10==0
            @assert !isnothing(savepath)
            tebd_results_interim = TEBDResults(hcat(entropy...), self_overlap,
                                hcat(corrs...), ϕ, ψ)
            save_structs(tebd_results_interim, save_path)
        end

    end

    return TEBDResults(hcat(entropy...), self_overlap,
                        hcat(corrs...), ϕ, ψ)
end


text:
module QSimulator

using Interpolations

import Base: getindex
import Base.+
using LinearAlgebra
using SparseArrays

eye(dim::Integer) = Matrix{Float64}(LinearAlgebra.I, dim, dim)
eye(::Type{T}, dim::Integer) where T = Matrix{T}(LinearAlgebra.I, dim, dim)
speye(dim::Integer) = sparse(one(Float64) * LinearAlgebra.I, dim, dim)
speye(::Type{T}, dim::Integer) where T = sparse(one(T) * LinearAlgebra.I, dim, dim)

include("base-types.jl")
include("systems.jl")
include("controls.jl")
include("interactions.jl")
include("composite-systems.jl")
include("dissipation.jl")
include("evolution.jl")
include("read_APS_file.jl")

end


text:
module Yolk

using Mixtape
import Mixtape: CompilationContext, allow, optimize!
using MacroTools: @capture, postwalk, rmlines, unblock
using CodeInfoTools
using CompilerPluginTools: inline_const!
using Core.Compiler: Const, is_pure_intrinsic_infer, intrinsic_nothrow, anymap, quoted, cfg_simplify!, compact!, adce_pass!

using Metatheory
using Metatheory.Library
using Metatheory.EGraphs

@metatheory_init()

include("methodtheory.jl")

struct YolkOptimizer <: CompilationContext
    theory::Vector{Rule}
    block_size_limit::Int
    opt::Bool
    YolkOptimizer(t) = new(t, 2, true)
    YolkOptimizer(t, b) = new(t, 2, b)
    YolkOptimizer() = new(Rule[])
end
allow(::YolkOptimizer, m::Module, args...) = true

shrink(ex::Expr, t) = t
shrink(ex::T, t) where T = Const(ex)

local_inline(g, stmts) = g
local_inline(g::GlobalRef, stmts) = CodeInfoTools.resolve(g)
function local_inline(g::Core.Compiler.SSAValue, stmts)
    g.id <= length(stmts) || return g
    local_inline(Core.Compiler.getindex(stmts, g.id)[:inst], stmts)
end
local_inline(ex::Expr, stmts) = Expr(ex.head, map(a -> local_inline(a, stmts), ex.args)...)

function extract!(ctx, ir)
    theory = ctx.theory
    for i in 1 : ctx.block_size_limit : length(ir.stmts)
        stmt = ir.stmts[i][:inst]
        type = ir.stmts[i][:type]
        if stmt isa Expr && stmt.head == :call
            if allow(ctx, CodeInfoTools.resolve(stmt.args[1]))
                r = local_inline(stmt, ir.stmts)
                G = EGraph(r)
                saturate!(G, theory)
                ex = Metatheory.EGraphs.extract!(G, astsize)
                if !(ex == r)
                    Core.Compiler.setindex!(ir.stmts[i], ex, :inst)
                    Core.Compiler.setindex!(ir.stmts[i], shrink(ex, type), :type)
                end
            end
        end
    end
    return ir
end

function optimize!(ctx::YolkOptimizer, b)
    ir = get_ir(b)
    if !ctx.opt
        ir = julia_passes!(ir, b.sv.src, b.sv)
    else
        ir = extract!(ctx, ir)
        ir = inline_const!(ir)
        ir = compact!(ir, true)
        ir = cfg_simplify!(ir)
        ir = julia_passes!(ir, b.sv.src, b.sv)
    end
    return ir
end

function opt(fn, tt::Type{T};
        emit_callable = false,
        ctx = YolkOptimizer(),
        opt = true) where T <: Tuple
    emit_callable && return jit(fn, tt; ctx = ctx, opt = opt)
    return emit(fn, tt; ctx = ctx, opt = opt)
end

export opt, YolkOptimizer, @methodtheory

end # module


text:
# ---
# title: 201. Bitwise AND of Numbers Range
# id: problem201
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/bitwise-and-of-numbers-range/description/>
# hidden: true
# ---
#
# Given a range [m, n] where 0 <= m <= n <= 2147483647, return the bitwise AND
# of all numbers in this range, inclusive.
#
# **Example 1:**
#
#
#
#     Input: [5,7]
#     Output: 4
#
#
# **Example 2:**
#
#
#
#     Input: [0,1]
#     Output: 0
#
#
## @lc code=start
using LeetCode

function range_bitwise_and(m::Int, n::Int)::Int
    i = 0
    while m != n
        m >>= 1
        n >>= 1
        i += 1
    end
    return n << i
end
## @lc code=end


text:
@testset "Analysis" begin
    @test true == true
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
using FeynmanDiagram
using Lehmann
using LinearAlgebra

const diagType = Ver4Diag
const Order = 3
const Circle = 100000

const kF = 1.919
const β = 25.0 / kF^2
const Λs = 1.0

function benchmark(tree, N, varK, varT)
    for i in 1:N
        ExprTree.evalNaive!(tree, varK, varT) #evaluate the expression tree
    end
end

function DiagTree.eval(id::BareGreenId, K, extT, varT)
    τin, τout = varT[extT[1]], varT[extT[2]]
    # ϵ = dot(K, K) - kF^2
    ϵ = K[1] * K[1] + K[2] * K[2] + K[3] * K[3] - kF^2
    τ = τout - τin
    if τ ≈ 0.0
        return Spectral.kernelFermiT(-1e-8, ϵ, β)
    else
        return Spectral.kernelFermiT(τ, ϵ, β)
    end
end

# DiagTree.eval(id::BareInteractionId, K, extT, varT) = 8π / (dot(K, K) + Λs)
DiagTree.eval(id::BareInteractionId, K, extT, varT) = 8π / (K[1] * K[1] + K[2] * K[2] + K[3] * K[3] + Λs)
# DiagTree.eval(id, K, extT, varT) = 1.0

diagPara(order) = GenericPara(diagType = diagType, innerLoopNum = order, hasTau = true,
    interaction = [FeynmanDiagram.Interaction(ChargeCharge, Instant),],  #instant charge-charge interaction
    # filter = [NoFock,])
    filter = [Girreducible,])

println("Build the diagrams into an experssion tree ...")
const para = [diagPara(o) for o in 1:Order]

#diagram of different orders
if diagType == SigmaDiag
    diags = [Parquet.sigma(para[i]) for i in 1:Order]
elseif diagType == PolarDiag
    diags = [Parquet.polarization(para[i]) for i in 1:Order]
elseif diagType == Ver4Diag
    diags = [Parquet.vertex4(para[i]) for i in 1:Order]
else
    error("not implemented!")
end
const extT = [diags[o].extT for o in 1:Order]                        #external tau of each diagram
println("Building tree")
const tree = [ExprTree.build(mergeby(diags[o]).diagram[1], verbose = 1) for o in 1:Order]     #experssion tree representation of diagrams
println("Done.")
# const tree = [ExprTree.build(DiagTree.optimize(mergeby(diags[o]).diagram[1])) for o in 1:Order]     #experssion tree representation of diagrams

varK = rand(3, 16)
varT = [rand() for i in 1:8]

benchmark(tree[1], 100, varK, varT)
println("Run evaluation for $Circle times.")
for i in 1:Order
    @time benchmark(tree[i], Circle, varK, varT)
end


text:
# Laplace distribution

export Laplace

@parameterized Laplace() ≪ (1/2) * Lebesgue(ℝ)

Laplace(nt::NamedTuple{(:μ,:σ)}) = Affine(nt, Laplace())
Laplace(nt::NamedTuple{(:μ,:ω)}) = Affine(nt, Laplace())
Laplace(nt::NamedTuple{(:σ,)}) = Affine(nt, Laplace())
Laplace(nt::NamedTuple{(:ω,)}) = Affine(nt, Laplace())
Laplace(nt::NamedTuple{(:μ,)}) = Affine(nt, Laplace())

@affinepars Laplace

function logdensity(d::Laplace{()} , x)
    return -abs(x)
end

Base.rand(rng::AbstractRNG, μ::Laplace{()}) = rand(rng, Dists.Laplace())

≪(::Laplace, ::Lebesgue{X}) where X <: Real = true

TV.as(::Laplace) = asℝ

distproxy(::Laplace{()}) = Dists.Laplace()


text:
using DemoCards
using DemoCards: democard, MarkdownDemoCard, DemoSection, DemoPage
using DemoCards: generate
using Test, ReferenceTests, Suppressor

# support both `include("runtests.jl")` and `include("test/runtests.jl")`
test_root = basename(pwd()) == "test" ? "." : "test"

cd(test_root) do
    include("types/card.jl")
    include("types/section.jl")
    include("types/page.jl")

    include("generate.jl")
    include("show.jl")
    include("utils.jl")
end

nothing


text:
using DataStructures # for OrderedDict
using JuMP

##################################################
# Variables
#mutable struct UCVariable
#  "Name of the variable."
#  name::Symbol
#  "What does the variable represent?"
#  description::String
#  "Global lower bound for the variable (may be adjusted later)."
#  lb::Float64
#  "Global upper bound for the variable (may be adjusted later)."
#  ub::Float64
#  "Is the variable integer-restricted?"
#  integer::Bool
#  "What are we indexing over?"*
#  " Recursive structure, e.g., [X,Y] means Y is a field in X,"*
#  " and [X,[Y1,Z],Y2] means Y1 and Y2 are fields in X and Z is a field in Y1.\n"*
#  " [ X, [Y,A,B], [Y,A,A], [Z,[D,E],F], T ]\n"*
#  " => [x, y1, y1.a, y1.b, y2, y2.a1, y2.a2, z, z.d, z.d.e, z.f, t]."
#  indices::Vector
#  "Function to add the variable; if this is missing, we will attempt to add the variable automatically using the `indices`. Signature should be (variable, model.vars.familyname, mip, instance)."
#  add_variable::Union{Function,Nothing}
#end # UCVariable

# TODO Above did not work for some reason
mutable struct UCVariable
  name::Symbol
  description::String
  lb::Float64
  ub::Float64
  integer::Bool
  indices::Vector
  add_variable::Union{Function,Nothing}
end

"""
It holds that x(t,t') = 0 if t' does not belong to 𝒢 = [t+DT, t+TC-1].
This is because DT is the minimum downtime, so there is no way x(t,t')=1 for t'<t+DT
and TC is the "time until cold" => if the generator starts afterwards, always has max cost.
"""
function add_downtime_arcs(var::UCVariable,
                           x::OrderedDict,
                           mip::JuMP.Model,
                           instance::UnitCommitmentInstance)
  T = instance.time
  for g in instance.units
    S = length(g.startup_categories)
    if S == 0
      continue
    end

    DT = g.min_downtime # minimum time offline
    TC = g.startup_categories[S].delay # time offline until totally cold

    for t1 = 1:T-1
      for t2 = t1+1:T
        # It holds that x(t,t') = 0 if t' does not belong to 𝒢 = [t+DT, t+TC-1]
        # This is because DT is the minimum downtime, so there is no way x(t,t')=1 for t'<t+DT
        # and TC is the "time until cold" => if the generator starts afterwards, always has max cost
        if (t2 < t1 + DT) || (t2 >= t1 + TC)
          continue
        end

        name = string(var.name, "[", g.name, ",", t1, ",", t2, "]")
        x[g.name, t1, t2] = @variable(mip,
                                      lower_bound=var.lb,
                                      upper_bound=var.ub,
                                      integer=var.integer,
                                      base_name=name)
      end # loop over time 2
    end # loop over time 1
  end # loop over units
end # add_downtime_arcs

"""
If there is a penalty specified for not meeting the reserve, then we add a reserve shortfall variable.
"""
function add_reserve_shortfall(var::UCVariable,
                               x::OrderedDict,
                               mip::JuMP.Model,
                               instance::UnitCommitmentInstance)
  T = instance.time
  for t = 1:T
    if instance.shortfall_penalty[t] > 1e-7
      name = string(var.name, "[", t, "]")
      x[t] = @variable(mip,
                       lower_bound=var.lb,
                       upper_bound=var.ub,
                       integer=var.integer,
                       base_name=name)
    end
  end # loop over time
end # add_reserve_shortfall

"""
Variables that the model may (or may not) use.

Note the relationship
  r_g(t) = bar{p}_g(t) - p_g(t)
         = bar{p}'_g(t) - p'_g(t)
"""
var_list = OrderedDict{Symbol,UCVariable}(
  :prod
    => UCVariable(:prod,
                  "[gen, t]; power from generator gen at time t; p_g(t) = p'_g(t) + g.min_power[t] * u_g(t)",
                  0., Inf, false,
                  [Unit, Time], nothing),
  :prod_above
    => UCVariable(:prod_above,
                  "[gen, t]; production above minimum required level; p'_g(t)",
                  0., Inf, false,
                  [Unit, Time], nothing ),
  :max_power_avail
    => UCVariable(:max_power_avail,
                  "[gen, t]; maximum power available from generator gen at time t; bar{p}_g(t) = p_g(t) + r_g(t)",
                  0., Inf, false,
                  [Unit, Time], nothing),
  :max_power_avail_above
    => UCVariable(:max_power_avail_above,
                  "[gen, t]; maximum power available above minimum from generator gen at time t; bar{p}'_g(t)",
                  0., Inf, false,
                  [Unit, Time], nothing),
  :segprod
    => UCVariable(:segprod,
                  "[gen, seg, t]; how much generator gen produces on segment seg in time t; p_g^l(t)",
                  0., Inf, false,
                  [ [Unit, CostSegment], Time], nothing),
  :reserve
    => UCVariable(:reserve,
                  "[gen, t]; reserves provided by gen at t; r_g(t)",
                  0., Inf, false,
                  [Unit, Time], nothing),
  :reserve_shortfall
    => UCVariable(:reserve_shortfall,
                  "[t]; reserve shortfall at gen at t; s_R(t)",
                  0., Inf, false,
                  [Time], add_reserve_shortfall),
  :is_on
    => UCVariable(:is_on,
                  "[gen, t]; is gen on at t; u_g(t)",
                  0., 1., true,
                  [Unit, Time], nothing),
  :switch_on
    => UCVariable(:switch_on,
                  "[gen, t]; indicator that gen will be turned on at t; v_g(t)",
                  0., 1., true,
                  [Unit, Time], nothing),
  :switch_off
    => UCVariable(:switch_off,
                  "[gen, t]; indicator that gen will be turned off at t; w_g(t)",
                  0., 1., true,
                  [Unit, Time], nothing),
  :net_injection
    => UCVariable(:net_injection,
                  "[bus.name, t]",
                  -1e100, Inf, false,
                  [Bus, Time], nothing),
  :curtail
    => UCVariable(:curtail,
                  "[bus.name, t]; upper bound is max load at the bus at time t",
                  0., Inf, false,
                  [Bus, Time], nothing),
  :flow
    => UCVariable(:flow,
                  "[violation.monitored_line.name, t]",
                  -1e100, Inf, false,
                  [Violation, Time], nothing),
  :overflow
    => UCVariable(:overflow,
                  "[transmission_line.name, t]; how much flow above the transmission limits (in MW) is allowed",
                  0., Inf, false,
                  [TransmissionLine, Time], nothing),
  :loads
    => UCVariable(:loads,
                  "[price_sensitive_load.name, t]; production to meet demand at a set price, if it is economically sensible, independent of the rest of the demand; upper bound is demand at this price at time t",
                  0., Inf, false,
                  [PriceSensitiveLoad, Time], nothing),
  :startup
    => UCVariable(:startup,
                  "[gen, startup_category, t]; indicator that generator g starts up in startup_category at time t; 𝛿_g^s(t)",
                  0., 1., true,
                  [ [Unit, StartupCategory], Time], nothing),
  :downtime_arc
    => UCVariable(:downtime_arc,
                  "[gen, t, t']; indicator for shutdown at t and starting at t'",
                  0., 1., true,
                  [Unit, Time, Time], add_downtime_arcs),
) # var_list

#var_symbol_list =
#  [
#   :prod_above,     # [gen, t], ≥ 0
#   :segprod,        # [gen, t, segment], ≥ 0
#   :reserve,        # [gen, t], ≥ 0
#   :is_on,          # [gen, t], binary
#   :switch_on,      # [gen, t], binary
#   :switch_off,     # [gen, t], binary
#   :net_injection,  # [bus.name, t], urs?
#   :curtail,        # [bus.name, t], domain [0, b.load[t]]
#   :overflow,       # [transmission_line.name, t], ≥ 0
#   :loads,          # [price_sensitive_load.name, t], domain [0, ps.demand[t]]
#   :startup         # [gen, t, startup_category], binary
#  ]

"""
For a particular UCElement, which is the field in UnitCommitmentInstance that this corresponds to?
This is used to determine indexing and ranges, e.g., `is_on` is indexed over Unit and Time,
so the variable `is_on` will range in the first index from 1 to length(instance.units)
and on the second index from 1 to instance.time.
"""
ind_to_field_dict = OrderedDict{Type{<:UCElement},Symbol}(
  Time                => :time,
  Bus                 => :buses,
  Unit                => :units,
  TransmissionLine    => :lines,
  PriceSensitiveLoad  => :price_sensitive_loads,
  CostSegment         => :cost_segments,
  StartupCategory     => :startup_categories,
) # ind_to_field_dict

"""
Take indices and convert them to fields of UnitCommitmentInstance.
"""
function ind_to_field(index::Union{Vector,Type{<:UCElement}}) :: Union{Vector,Symbol}
  if isa(index, Type{<:UCElement})
    return ind_to_field_dict[index]
  else
    return [ ind_to_field(t) for t in index ]
  end
end # ind_to_field

function num_indices(v) :: Int64
  if !isa(v, Array)
    return 1
  else
    return sum(num_indices(v[i]) for i in 1:length(v))
  end
end # num_indices

"""
Can return
  * UnitRange -> iterate over this range
  * Array{UnitRange} -> cross product of the ranges in the array
  * Tuple(UnitRange, Array{UnitRange}) -> the array length should be the same as the range of the UnitRange
"""
function get_indices_tuple(obj::Any, fields::Union{Symbol,Vector,Nothing} = nothing)
  if isa(fields, Symbol)
    return get_indices_tuple(getfield(obj,fields))
  end
  if fields == nothing || (isa(fields,Array) && length(fields) == 0)
    if isa(obj, Array)
      return UnitRange(1,length(obj))
    elseif isa(obj, Int)
      return UnitRange(1,obj)
    else
      return UnitRange{Int64}(0:-1)
      #return UnitRange(1,1)
    end
  end

  if isa(obj,Array)
    indices = (
               UnitRange(1,length(obj)),
               ([
                 isa(f,Array) ? get_indices_tuple(getfield(x, f[1]), f[2:end]) : get_indices_tuple(getfield(x, f))
                 for x in obj
                ] for f in fields)...
              )
    #         more_indices = ([
    #                 isa(f,Array) ? get_indices_tuple(getfield(x, f[1]), f[2:end]) : get_indices_tuple(getfield(x, f))
    #                 for x in obj
    #               ] for f in fields
    #             )
    #         indices = (UnitRange(1,length(obj)),more_indices...)
  else
    indices = ()
    for f in fields
      if isa(f,Array)
        indices = (indices..., get_indices_tuple(getfield(obj, f[1]), f[2:end]))
      else
        indices = (indices..., get_indices_tuple(obj,f))
      end
    end
    #         indices = (
    #             isa(f,Array) ? get_indices_tuple(getfield(obj, f[1]), f[2:end]) : get_indices_tuple(getfield(obj, f))
    #                 for f in fields
    #             )
    #         (
    #             isa(f,Array) ? get_indices_tuple(getfield(obj, f[1]), f[2:end]) : get_indices_tuple(getfield(obj, f))
    #             for f in fields
    #         )
    #         indices = (indices...,)
  end # check if obj is Array or not

  return indices
end # get_indices_tuple

function loop_over_indices(indices::Any)
  loop = nothing
  should_print = false

  if isa(indices, UnitRange)
    loop = indices
  elseif isa(indices, Array{UnitRange{Int64}}) || isa(indices, Tuple{Int, UnitRange})
    loop = Base.product(Tuple(indices)...)
  elseif isa(indices, Tuple{UnitRange, Array})
    loop = ()
    for t in zip(indices...)
      loop = (loop..., loop_over_indices(t)...)
    end
  elseif isa(indices,Tuple)
    loop = ()
    for i in indices
      loop = (loop..., loop_over_indices(i))
    end
    loop = Base.product(loop...)
  else
    error("Why are we here?")
    #loop = Base.product(loop_over_indices(indices)...)
  end

  if should_print
    for i in loop
      @show i
    end
  end
  return loop
end # loop_over_indices

function expand_tuple(x::Tuple)
  y = ()
  for i in x
    if isa(i, Tuple)
      y = (y..., expand_tuple(i)...)
    else
      y = (y..., i)
    end
  end
  return y
end # expand_tuple

function expand_tuple(X::Array{<:Tuple})
  return [ expand_tuple(x) for x in X ]
end # expand_tuple

function get_indices(x::Array)
  return expand_tuple(x)
end
function get_indices(x::Base.Iterators.ProductIterator)
  return get_indices(collect(x))
end

"""
Access `t.f`, special terminal case of `get_nested_field`.
"""
function get_nested_field(t::Any, f::Symbol)
  return getfield(t,f)
end # get_nested_field

"""
Access `t.f`, where `f` could be a subfield.
"""
function get_nested_field(t::Any, f::Vector{Symbol})
  if length(f) > 1
    return get_nested_field(getfield(t,f[1]), f[2:end])
  else
    return getfield(t,f[1])
  end
end # get_nested_field

"""
Given a set of indices of UCVariable, e.g., [[X,Y],T],
and a UnitCommitmentInstance instance,
if we want to access the field corresponding to Y,
then we call get_nested_field(instance, [[X,Y],T], 2, (4,3)),
which will return instance.X[4].Y[3] if Y is a vector,
and just instance.X[4].Y otherwise.

===
Termination Conditions

If `i` <= 0, then we only care about instance, and not the field.
If `field` is a Symbol and `i` >= 1, then we want to explore instance.field (index t[i] or t).
Note that if `i` >= 1, then `field` must be a symbol.
If `i` == 1, then `t` can be an Int.

===
Parameters
  * instance::Any --> all fields will be from instance, or nested fields of fields of instance.
  * field::Union{Vector,Symbol,Nothing} --> either the field we want to access, or a vector of fields, and we will want field[i].
  * i::Int --> which field to access.
  * t::Tuple --> how to go through the fields of instance to get the right field, length needs to be at least `i`.
"""
function get_nested_field(instance::Any, field::Union{Vector,Symbol,Nothing}, i::Int, t::Union{Tuple, Int})
  # Check data
  if isa(field, Vector)
    if i >= 2 && (!isa(t,Tuple) || length(t) < i)
      error("Tuple of indices to get nested field needs to be at least the length of the index we want to get.")
    end
  end

  if isa(field, Symbol) || i <= 0
    # i = 0 can happen in the recursive call
    # What it means is that we do not want a field of the instance, but the instance itself
    # TODO handle other iterable types and empty arrays
    f = (isa(field, Symbol) && i >= 1) ? getfield(instance, field) : instance
    if isa(f,Vector)
      if length(f) == 0
        error("Trying to iterate over empty field!")
      else
        return isa(t,Int) ? f[t] : f[t[i]]
      end
    else
      return f
    end
  end # check termination conditions (f is field or i <= 0)

  # Loop over the fields until we find where index i is located
  # It may be nested inside an array, so that is why we recurse
  start_ind = 0
  for f in field
    curr_len = isa(f, Vector) ? length(f) : 1
    if start_ind + curr_len >= i
      if isa(f, Vector)
        new_field_is_iterable = isa(getfield(instance, f[1]), Vector)
        if new_field_is_iterable
          return get_nested_field(getfield(instance, f[1])[t[start_ind+1]], f[2:end], i - start_ind - 1, isa(t,Tuple) ? t[start_ind+2:end] : t)
        else
          return get_nested_field(getfield(instance, f[1]), f[2:end], i - start_ind - 1, isa(t,Tuple) ? t[start_ind+2:end] : t)
        end
      else
        # f is hopefully a symbol...
        return get_nested_field(instance, f, 1, isa(t,Tuple) ? t[start_ind+1] : t)
      end
    end
    start_ind += curr_len
  end

  return nothing
end # get_nested_field

#"""
#Get ranges for the indices of a UCVariable along dimension `i`,
#making sure that the right fields ranges are calculated via `get_nested_field` and `ind_to_field`.
#"""
#function get_range(arr::UCVariable, instance::UnitCommitmentInstance, i::Int) :: UnitRange
#  arr = ind_to_field[var.indices[i]]
#  f = get_nested_field(instance, arr)
#  if isa(f, Array)
#    return 1:length(f)
#  elseif isa(f, Int)
#    return 1:f
#  else
#    error("Unknown type to generate UnitRange from: ", typeof(f))
#  end
#end # get_range

export UCVariable


text:
"""
    gaussian_kernel_hyperparameters
This contains the hyperparameters for the gaussian kernel. The functional form for correlation between x_p and x_q   is   w_0 exp( -0.5 * \\sum_{d=1}^D ((x_{p,d} - x_{q,d})/w_d)^2)
    where D is the number of dimensions. and each term of the summation is a different dimension of x_p and x_q.
"""
struct GaussianKernelHyperparameters{T<:Real}
    w_0::T
    w_i::Array{T,1}
    function GaussianKernelHyperparameters(w_0::S,w_i::Array{G,1}) where S<:Real where G<:Real
        promo = promote_type(S,G)
        return new{promo}(promo(w_0), promo.(w_i))
    end
    function GaussianKernelHyperparameters{S}(aa::GaussianKernelHyperparameters{G}) where S<:Real where G<:Real
        return new{S}(S(aa.w_0), S.(aa.w_i))
    end
end

"""
    KrigingModel
    Our kriging model is defined by a level \\mu, a set of X coordinates and corresponding y values, a set of hyperparameters and a noise level.
In addition we save a number of matrices in order to faciliate computations. One is
This contains the hyperparameters for the gaussian kernel.
"""

struct KrigingModel{T<:Real}
    X::AbstractArray{T,2}
    f::AbstractArray{T,1}
    hyper::GaussianKernelHyperparameters{T}
    cov_func::Function
    noise::T
    # Now we have stuff that is calculated from the above. Saved in the model purely for efficiency reasons.
    K::Hermitian{T}
    invK::Hermitian{T}
    mu::T
    sigma2::T
    y_m_mu::Array{T,1}
    function KrigingModel(X::AbstractArray{T,2}, f::AbstractArray{R,1}, hyper::GaussianKernelHyperparameters{S};
                          cov_func::Function = gaussian_kernel, noise::U = 1000*eps(),
                          K::Hermitian{W} = K_matrix(X, cov_func, hyper, noise), invK::Hermitian{Z} = LinearAlgebra.inv(K)) where T<:Real where R<:Real where S<:Real where U<:Real where W<:Real where Z<:Real
            tt = promote_type(T, R, S, U, W, Z)

            dim = size(invK)[1]
            mu = tt(sum(invK * f) / sum(invK * ones(dim)))
            y_m_mu = f .- mu
            sigma2 = tt((y_m_mu' * (invK * y_m_mu)) / dim)
            return new{tt}(Array{tt,2}(X), Array{tt,1}(f), GaussianKernelHyperparameters{tt}(hyper) , cov_func, tt(noise), LinearAlgebra.Hermitian(Array{tt,2}(K)), LinearAlgebra.Hermitian(Array{tt,2}(invK)), mu, sigma2, y_m_mu)
    end
end


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

mutable struct TopologyResource <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    associations::Any # spec type: Union{ Nothing, Vector{TopologyAssociation} } # spec name: associations

    function TopologyResource(;name=nothing, id=nothing, location=nothing, associations=nothing)
        o = new()
        validate_property(TopologyResource, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(TopologyResource, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(TopologyResource, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(TopologyResource, Symbol("associations"), associations)
        setfield!(o, Symbol("associations"), associations)
        o
    end
end # type TopologyResource

const _property_map_TopologyResource = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("id")=>Symbol("id"), Symbol("location")=>Symbol("location"), Symbol("associations")=>Symbol("associations"))
const _property_types_TopologyResource = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("id")=>"String", Symbol("location")=>"String", Symbol("associations")=>"Vector{TopologyAssociation}")
Base.propertynames(::Type{ TopologyResource }) = collect(keys(_property_map_TopologyResource))
Swagger.property_type(::Type{ TopologyResource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_TopologyResource[name]))}
Swagger.field_name(::Type{ TopologyResource }, property_name::Symbol) =  _property_map_TopologyResource[property_name]

function check_required(o::TopologyResource)
    true
end

function validate_property(::Type{ TopologyResource }, name::Symbol, val)
end


text:
using Test

@testset "KnetMetrics" begin
    include("classification-metrics.jl")
    include("regression-metrics.jl")
    include("pairwise-metrics.jl")
end