text:
using Flux
using ConditionalDists
using GenerativeModels
import GenerativeModels: VAE
using ValueHistories
using MLDataPattern: RandomBatches
using Distributions
using DistributionsAD
using StatsBase
using Random

"""
	safe_softplus(x::T)

Safe version of softplus.
"""
safe_softplus(x::T) where T  = softplus(x) + T(0.000001)

"""
	init_vamp(pseudoinput_mean,k::Int)

Initializes the VAMP prior from a mean vector and number of components.
"""
function init_vamp(pseudoinput_mean, k::Int)
	T = eltype(pseudoinput_mean)
	pseudoinputs = T(1) .* randn(T, size(pseudoinput_mean)[1:end-1]..., k) .+ pseudoinput_mean
	VAMP(pseudoinputs)
end

"""
	vae_constructor(;idim::Int=1, zdim::Int=1, activation = "relu", hdim=128, nlayers::Int=3,
		init_seed=nothing, prior="normal", pseudoinput_mean=nothing, k=1, kwargs...)

Constructs a classical variational autoencoder.

# Arguments
	- `idim::Int`: input dimension.
	- `zdim::Int`: latent space dimension.
	- `activation::String="relu"`: activation function.
	- `hdim::Int=128`: size of hidden dimension.
	- `nlayers::Int=3`: number of decoder/encoder layers, must be >= 3.
	- `init_seed=nothing`: seed to initialize weights.
	- `prior="normal"`: one of ["normal", "vamp"].
	- `pseudoinput_mean=nothing`: mean of data used to initialize the VAMP prior.
	- `k::Int=1`: number of VAMP components.
	- `var="scalar"`: decoder covariance computation, one of ["scalar", "diag"].
"""
function vae_constructor(;idim::Int=1, zdim::Int=1, activation="relu", hdim=128, nlayers::Int=3,
	init_seed=nothing, prior="normal", pseudoinput_mean=nothing, k=1, var="scalar", kwargs...)
	(nlayers < 3) ? error("Less than 3 layers are not supported") : nothing

	# if seed is given, set it
	(init_seed != nothing) ? Random.seed!(init_seed) : nothing

	# construct the model
	# encoder - diagonal covariance
	encoder_map = Chain(
		build_mlp(idim, hdim, hdim, nlayers-1, activation=activation)...,
		ConditionalDists.SplitLayer(hdim, [zdim, zdim], [identity, safe_softplus])
		)
	encoder = ConditionalMvNormal(encoder_map)

	# decoder - we will optimize only a shared scalar variance for all dimensions
	if var=="scalar"
		decoder_map = Chain(
			build_mlp(zdim, hdim, hdim, nlayers-1, activation=activation)...,
			ConditionalDists.SplitLayer(hdim, [idim, 1], [identity, safe_softplus])
			)
	else
		decoder_map = Chain(
				build_mlp(zdim, hdim, hdim, nlayers-1, activation=activation)...,
				ConditionalDists.SplitLayer(hdim, [idim, idim], [identity, safe_softplus])
				)
	end
	decoder = ConditionalMvNormal(decoder_map)

	# prior
	if prior == "normal"
		prior_arg = zdim
	elseif prior == "vamp"
		(pseudoinput_mean === nothing) ? error("if `prior=vamp`, supply pseudoinput array") : nothing
		prior_arg = init_vamp(pseudoinput_mean, k)
	end

	# reset seed
	(init_seed !== nothing) ? Random.seed!() : nothing

	# get the vanilla VAE
	model = VAE(prior_arg, encoder, decoder)
end

"""
	conv_vae_constructor(idim=(2,2,1), zdim::Int=1, activation="relu", hdim=1024, kernelsizes=(1,1),
		channels=(1,1), scalings=(1,1), init_seed=nothing, prior="normal", pseudoinput_mean=nothing,
		k=1, batchnorm=false, kwargs...)

Constructs a classical variational autoencoder.

# Arguments
	- `idim::Int`: input dimension.
	- `zdim::Int`: latent space dimension.
	- `activation::String="relu"`: activation function.
	- `hdim::Int=1024`: size of hidden dimension.
	- `kernelsizes=(1,1)`: kernelsizes in consequent layers.
	- `channels=(1,1)`: number of channels.
	- `scalings=(1,1)`: scalings in subsewuent layers.
	- `init_seed=nothing`: seed to initialize weights.
	- `prior="normal"`: one of ["normal", "vamp"].
	- `pseudoinput_mean=nothing`: mean of data used to initialize the VAMP prior.
	- `k::Int=1`: number of VAMP components.
	- `batchnorm=false`: use batchnorm (discouraged).
"""
function conv_vae_constructor(;idim=(2,2,1), zdim::Int=1, activation="relu", hdim=1024,
	kernelsizes=(1,1), channels=(1,1), scalings=(1,1),
	init_seed=nothing, prior="normal", pseudoinput_mean=nothing, k=1,
	batchnorm=false, kwargs...)
	# if seed is given, set it
	(init_seed != nothing) ? Random.seed!(init_seed) : nothing

	# construct the model
	# encoder - diagonal covariance
	encoder_map = Chain(
		conv_encoder(idim, hdim, kernelsizes, channels, scalings; activation=activation,
			batchnorm=batchnorm)...,
		ConditionalDists.SplitLayer(hdim, [zdim, zdim], [identity, safe_softplus])
		)
	encoder = ConditionalMvNormal(encoder_map)

	# decoder - we will optimize only a shared scalar variance for all dimensions
	# also, the decoder output will be vectorized so the usual logpdfs vcan be used
	vecdim = reduce(*,idim[1:3]) # size of vectorized data
	decoder_map = Chain(
		conv_decoder(idim, zdim, reverse(kernelsizes), reverse(channels), reverse(scalings),
			activation=activation, vec_output=true, vec_output_dim=nothing, batchnorm=batchnorm)...,
		ConditionalDists.SplitLayer(vecdim, [vecdim, 1], [identity, safe_softplus])
		)
	decoder = ConditionalMvNormal(decoder_map)

	# prior
	if prior == "normal"
		prior_arg = zdim
	elseif prior == "vamp"
		(pseudoinput_mean == nothing) ? error("if `prior=vamp`, supply pseudoinput array") : nothing
		prior_arg = init_vamp(pseudoinput_mean, k)
	end

	# reset seed
	(init_seed != nothing) ? Random.seed!() : nothing

	# get the vanilla VAE
	model = VAE(prior_arg, encoder, decoder)
end

"""
	StatsBase.fit!(model::VAE, data::Tuple, loss::Function; max_train_time=82800, lr=0.001,
		batchsize=64, patience=30, check_interval::Int=10, kwargs...)
"""
function StatsBase.fit!(model::VAE, data::Tuple, loss::Function; max_iters=10000, max_train_time=82800,
	lr=0.001, batchsize=64, patience=30, check_interval::Int=10, kwargs...)
	history = MVHistory()
	opt = ADAM(lr)

	tr_model = deepcopy(model)
	ps = Flux.params(tr_model)
	_patience = patience

	tr_x = data[1][1]
	# i know this could be done generally for all sizes but it is very ugly afaik
	if ndims(tr_x) == 2
		val_x = data[2][1][:,data[2][2] .== 0]
	elseif ndims(tr_x) == 4
		val_x = data[2][1][:,:,:,data[2][2] .== 0]
	else
		error("not implemented for other than 2D and 4D data")
	end
	val_N = size(val_x,ndims(val_x))

	# on large datasets, batching loss is faster
	best_val_loss = Inf
	i = 1
	start_time = time() # end the training loop after 23hrs
	for batch in RandomBatches(tr_x, batchsize)
		# batch loss
		batch_loss = 0f0
		gs = gradient(() -> begin
			batch_loss = loss(tr_model,batch)
		end, ps)
	 	Flux.update!(opt, ps, gs)

		push!(history, :training_loss, i, batch_loss)
		if mod(i, check_interval) == 0

			# validation/early stopping
			val_loss = (val_N > 5000) ? loss(tr_model, val_x, 256) : loss(tr_model, val_x)

			@info "$i - loss: $(batch_loss) (batch) | $(val_loss) (validation)"

			if isnan(val_loss) || isnan(batch_loss)
				error("Encountered invalid values in loss function.")
			end

			push!(history, :validation_likelihood, i, val_loss)

			if val_loss < best_val_loss
				best_val_loss = val_loss
				_patience = patience

				# this should save the model at least once
				# when the validation loss is decreasing
				model = deepcopy(tr_model)
			else # else stop if the model has not improved for `patience` iterations
				_patience -= 1
				if _patience == 0
					@info "Stopped training after $(i) iterations."
					break
				end
			end
		end
		if (time() - start_time > max_train_time) | (i > max_iters) # stop early if time is running out
			model = deepcopy(tr_model)
			@info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours."
			break
		end
		i += 1
	end
	# again, this is not optimal, the model should be passed by reference and only the reference should be edited
	(history=history, iterations=i, model=model, npars=sum(map(p->length(p), Flux.params(model))))
end


text:
include("../../src/EDKit.jl")
using .EDKit
using LinearAlgebra
using DelimitedFiles
using Test
using PyCall
np = pyimport("numpy")

@testset "PXP-L-26" begin
    L, k, p = 26, 0, 1
    mat = begin
        P = Diagonal([1, 1, 1, 0, 1, 1, 0, 0])
        X = [0 1; 1 0]
        P * kron(I(2), X, I(2)) * P
    end
    pxpf(v::Vector{<:Integer}) = all(v[i]==0 || v[mod(i, length(v))+1]==0 for i=1:length(v))
    print("Searching time :")
    @time basis = TranslationParityBasis(pxpf, k, p, L, threaded=true)
    H = trans_inv_operator(mat, 2, basis) |> Array
    E = np.linalg.eigvalsh(H)
    PXPE = readdlm("evals_periodic_N26_k0_p0.npy.txt")
    @test norm(PXPE-E) ≈ 0.0 atol=1e-10
end

@testset "PXP-L-28" begin
    L, k, p = 28, 0, 1
    mat = begin
        P = Diagonal([1, 1, 1, 0, 1, 1, 0, 0])
        X = [0 1; 1 0]
        P * kron(I(2), X, I(2)) * P
    end
    pxpf(v::Vector{<:Integer}) = all(v[i]==0 || v[mod(i, length(v))+1]==0 for i=1:length(v))
    print("Searching time :")
    @time basis = TranslationParityBasis(pxpf, k, p, L, threaded=true)
    H = trans_inv_operator(mat, 2, basis) |> Array
    E = np.linalg.eigvalsh(H)
    PXPE = readdlm("evals_periodic_N28_k0_p0.npy.txt")
    @test norm(PXPE-E) ≈ 0.0 atol=1e-10
end

ext: jl


text:
push!(LOAD_PATH, abspath("./src"))

using QMPS
using Test, Yao
import Random: seed!
import Statistics: mean

@testset "CircuitBuilder.jl + MPSC.jl" begin
    #####Test Functions#####
    function cr_test(reg::ArrayReg, cr::CompositeBlock, crt::CompositeBlock)
        n = nqubits(reg)
        seed!(seedNum)
        ExpReg1x = expect(opx(n), copy(reg) |> cr) |> mean |> real
        seed!(seedNum)
        ExpReg1y = expect(opy(n), copy(reg) |> cr) |> mean |> real
        seed!(seedNum)
        ExpReg1z = expect(opz(n), copy(reg) |> cr) |> mean |> real
        seed!(seedNum)
        ExpReg2x = expect(opx(n), copy(reg) |> crt) |> mean |> real
        seed!(seedNum)
        ExpReg2y = expect(opy(n), copy(reg) |> crt) |> mean |> real
        seed!(seedNum)
        ExpReg2z = expect(opz(n), copy(reg) |> crt) |> mean |> real
        @test ExpReg1x ≈ ExpReg2x
        @test ExpReg1y ≈ ExpReg2y
        @test ExpReg1z ≈ ExpReg2z
    end

    function ce_test(reg::ArrayReg, ce::CompositeBlock, cet::CompositeBlock)
        reg3 = copy(reg) |> ce
        reg4 = copy(reg) |> cet
        @test reg3 ≈ reg4
    end

    function re_test(regr::ArrayReg, rege::ArrayReg, cr::CompositeBlock, ce::CompositeBlock, vBit::Int64)
        nr = nqubits(regr)
        ne = nqubits(rege)
        seed!(seedNum)
        ExpRegrx = expect(opVx(nr, vBit), copy(regr) |> cr) |> mean |> real
        seed!(seedNum)
        ExpRegry = expect(opVy(nr, vBit), copy(regr) |> cr) |> mean |> real
        seed!(seedNum)
        ExpRegrz = expect(opVz(nr, vBit), copy(regr) |> cr) |> mean |> real
        seed!(seedNum)
        ExpRegex = expect(opVx(ne, vBit), copy(rege) |> ce) |> mean |> real
        seed!(seedNum)
        ExpRegey = expect(opVy(ne, vBit), copy(rege) |> ce) |> mean |> real
        seed!(seedNum)
        ExpRegez = expect(opVz(ne, vBit), copy(rege) |> ce) |> mean |> real
        @test isapprox(ExpRegrx, ExpRegex, atol = 0.5e-2)
        @test isapprox(ExpRegry, ExpRegey, atol = 0.5e-2)
        @test isapprox(ExpRegrz, ExpRegez, atol = 0.5e-2)
    end

    opx(n::Int64) = chain(n, [put(n, i=>X) for i=1:n])
    opy(n::Int64) = chain(n, [put(n, i=>Y) for i=1:n])
    opz(n::Int64) = chain(n, [put(n, i=>Z) for i=1:n])
    opVx(n::Int64, V::Int64) = chain(n, [put(n, i=>X) for i=1:V])
    opVy(n::Int64, V::Int64) = chain(n, [put(n, i=>Y) for i=1:V])
    opVz(n::Int64, V::Int64) = chain(n, [put(n, i=>Z) for i=1:V])

    #####CircuitBuilder.jl#####
    seedNum = 1234
    # DCbuilder(nBit::Int64, depth::Int64)
    ## Preparing section
    dc = DCbuilder(4,2)
    head  = chain(4, chain(4, put(4,4=>Rx(0)), put(4,3=>Rx(0)), put(4,2=>Rx(0)), put(4,1=>Rx(0))),
                     chain(4, put(4,4=>Rz(0)), put(4,3=>Rz(0)), put(4,2=>Rz(0)), put(4,1=>Rz(0))),
                     chain(4, control(4,4,3=>X), control(4,3,2=>X), control(4,2,1=>X), control(4,1,4=>X)))
    block = chain(4, chain(4, put(4,4=>Rz(0)), put(4,3=>Rz(0)), put(4,2=>Rz(0)), put(4,1=>Rz(0))),
                     chain(4, put(4,4=>Rx(0)), put(4,3=>Rx(0)), put(4,2=>Rx(0)), put(4,1=>Rx(0))),
                     chain(4, put(4,4=>Rz(0)), put(4,3=>Rz(0)), put(4,2=>Rz(0)), put(4,1=>Rz(0))),
                     chain(4, control(4,4,3=>X), control(4,3,2=>X), control(4,2,1=>X), control(4,1,4=>X)))
    tail  = chain(4, chain(4, put(4,4=>Rz(0)), put(4,3=>Rz(0)), put(4,2=>Rz(0)), put(4,1=>Rz(0))),
                     chain(4, put(4,4=>Rx(0)), put(4,3=>Rx(0)), put(4,2=>Rx(0)), put(4,1=>Rx(0))))
    Cblock = chain(4, block, block)
    body = chain(4, head, block, tail)
    ## Testing section
    @test head == dc.head
    @test block == dc.block
    @test tail == dc.tail
    @test Cblock == dc.circuit
    @test body == dc.fullbody

    # MPSbuilder(nBitA::Int64, vBit::Int64, rBit::Int64, blockT::String)
    ## Preparing section
    regCSr = rand_state(2, nbatch=5000)
    regCSe = rand_state(4)
    regCSr0 = zero_state(2, nbatch=10000)
    regCSe0 = zero_state(4, nbatch=10000)
    mpsCS = MPSbuilder(4, 1, 1, "CS")
    CScr = mpsCS.circuit
    CSce = mpsCS.cExtend
    CScrt = chain(2, chain(2, repeat(2, H, (2,1)), control(2, 1, 2=>Z), Measure(2, locs=2, resetto=0)),
                     chain(2, put(2, (2,1)=>SWAP), put(2, 1=>H), control(2, 1, 2=>Z), Measure(2, locs=2, resetto=0)),
                     chain(2, put(2, (2,1)=>SWAP), put(2, 1=>H), control(2, 1, 2=>Z)))
    CScet = chain(4, repeat(4, H, (4,3,2,1)), control(4, 3, 4=>Z), control(4, 2, 3=>Z), control(4, 1, 2=>Z))
    ## Testing section
    print(IOBuffer(), mpsCS)
    re_test(regCSr0, regCSe0, CScr, CSce, 1)
    re_test(regCSr0, regCSe0, CScrt, CScet, 1)
    cr_test(regCSr, CScr, CScrt)
    ce_test(regCSe, CSce, CScet)

    # MPSbuilder(nBitA::Int64, vBit::Int64, rBit::Int64, blockT::Tuple{String, Int64})
    ## Preparing section
    regDCr = rand_state(4, nbatch=5000)
    regDCe = rand_state(6)
    regDCr0 = zero_state(4, nbatch=10000)
    regDCe0 = zero_state(6, nbatch=10000)
    mpsDC = MPSbuilder(6, 2, 2, ("DC", 2))
    DCcr = mpsDC.circuit
    DCce = mpsDC.cExtend
    Cb1 = deepcopy(Cblock) |> markDiff
    Cb2 = deepcopy(Cblock) |> markDiff
    dispatch!(Cb1, :random)
    dispatch!(Cb2, :random)
    pars1 = getiparams.(content.(collect_blocks(QMPS.QDiff, Cb1)))
    pars2 = getiparams.(content.(collect_blocks(QMPS.QDiff, Cb2)))
    dispatch!(DCcr[1], pars1)
    dispatch!(DCcr[2], pars2)
    swap1 = chain(4, put(4, (2,3)=>SWAP), put(4, (3,4)=>SWAP), put(4, (1,2)=>SWAP), put(4, (2,3)=>SWAP))
    swap2 = chain(4, put(4, (3,2)=>SWAP), put(4, (2,1)=>SWAP), put(4, (4,3)=>SWAP), put(4, (3,2)=>SWAP))
    DCcrt = chain(4, Cb1, swap1, Measure(4, locs=(3,4), resetto=0), swap2, Cb2, swap1)
    DCcet = chain(6, subroutine(6, chain(4, Cb1, swap1), 3:6), subroutine(6, chain(4, Cb2, swap1), 1:4))
    ## Testing section
    print(IOBuffer(), mpsDC)
    re_test(regDCr0, regDCe0, DCcr, DCce, 2)
    re_test(regDCr0, regDCe0, DCcrt, DCcet, 2)
    cr_test(regDCr, DCcr, DCcrt)
    ce_test(regDCe, DCce, DCcet)

    #####MPSC.jl#####
    # MPSpar(nBitA::Int64, vBit::Int64, rBit::Int64)
    ## Preparing section
    nA = 13
    v = 4
    r = 3
    d = 2
    par = MPSpar(nA, v, r)
    ## Testing section
    @test par.nBitA == nA
    @test par.vBit == v
    @test par.rBit == r
    @test par.nBit == v+r
    @test par.nBlock == Int((nA - v) / r)

    # MPSDCpar(circuit::ChainBlock)
    ## Preparing section
    mpsDC2 = MPSbuilder(nA, v, r, ("DC", d))
    par2 = MPSDCpar(mpsDC2.circuit)
    par3 = MPSDCpar(mpsDC2.cExtend)
    ## Testing section
    @test par2.nBitA == nA
    @test par2.vBit == v
    @test par2.rBit == r
    @test par2.nBit == v+r
    @test par2.nBlock == length(mpsDC2.circuit)
    @test par2.depth == d
    @test par3.nBitA == nA
    @test par3.vBit == v
    @test par3.rBit == r
    @test par3.nBit == v+r
    @test par3.nBlock == length(mpsDC2.cExtend)
    @test par3.depth == d

    # MPSC(blockT, nBitA::Int64, vBit::Int64, rBit::Int64=1; dBlocksPar=0)
    ## blockT = "CS"
    ### Preparing section
    mpssCS = MPSC("CS", 4, 1, 1)
    ### Testing section
    re_test(regCSr0, regCSe0, mpssCS.circuit, mpssCS.cExtend, 1)
    cr_test(regCSr, mpssCS.circuit, CScrt)
    ce_test(regCSe, mpssCS.cExtend, CScet)
    @test mpssCS.mpsBlocks == [chain(2, repeat(2, H, (2,1)), control(2, 1, 2=>Z)),
                               chain(2, put(2, (2,1)=>SWAP), put(2, 1=>H), control(2, 1, 2=>Z)),
                               chain(2, put(2, (2,1)=>SWAP), put(2, 1=>H), control(2, 1, 2=>Z))]
    @test mpssCS.cEBlocks == CScet.blocks
    @test mpssCS.nBit == nqubits(CScrt)
    @test mpssCS.nBlock == length(mpssCS.circuit)
    ## blockT = ("DC", 2)
    ### Preparing section
    mpssDC = MPSC(("DC", 2), 6, 2, 2)
    reg2DCe = rand_state(6)
    seed!(seedNum)
    reg2DCe_1 = copy(reg2DCe)|> mpssDC.cExtend
    seed!(seedNum)
    reg2DCe_2 = copy(reg2DCe)|> DCcet
    seed!(seedNum)
    regDCr0_1 = copy(regDCr0)|> mpssDC.circuit
    seed!(seedNum)
    regDCr0_2 = copy(regDCr0)|> DCcrt
    dG2 = collect_blocks(QMPS.QDiff, mpssDC.circuit)
    dispatch!.(dG2, [pars1;pars2])
    ### Testing section
    @test (reg2DCe_1 ≈ reg2DCe_2) == false
    @test (regDCr0_1 ≈ regDCr0_2) == false
    cr_test(regDCr0, mpssDC.circuit, DCcrt)
    ce_test(reg2DCe, mpssDC.cExtend, DCcet)
    re_test(regDCr0, regDCe0, mpssDC.circuit, mpssDC.cExtend, 2)
    @test mpssDC.mpsBlocks == [chain(4, Cb1, swap1), chain(4, swap2, Cb2, swap1)]
    @test mpssDC.cEBlocks == DCce.blocks
    @test mpssDC.dGates == collect_blocks(QMPS.QDiff, DCce)
    @test mpssDC.nBit == nqubits(DCcr)
    @test mpssDC.nBlock == length(mpssDC.circuit)
    @test parameters(MPSC(("DC", 1), 3, 1, 1, dBlocksPar=[1.0:12.0;]).circuit) == [1.0:12.0;]

    #####Error Exception Tests#####
    @test_throws ErrorException MPSbuilder(6, 3, 2, ("DC", 2))
    @test_throws ErrorException MPSpar(5, 1, 3)
    @test_throws ErrorException mpssDC2 = MPSC(("DC", 1), 5, 2, 2)

    @test_throws ErrorException MPSbuilder(4, 2, 1, "CS")
    @test_throws ErrorException MPSbuilder(5, 1, 2, "CS")

    @test_throws ErrorException MPSDCpar(chain(chain(4, Cb1, put(4, 1=>X))))
    @test_throws ErrorException mpssDC2 = MPSC("CS", 3, 1, 1, dBlocksPar=[1:0.5:13;])
    @test_throws ErrorException mpssDC2 = MPSC(("DC", 1), 3, 1, 1, dBlocksPar=[1:0.5:13;])
end

@testset "Diff.jl" begin
    seedNum = 1234
    # Testing basic Yao-compatible functions for QDiff{GT, N} and markDiff(block::AbstractBlock).
    seed!(seedNum)
    n = 3
    c = chain(n, put(n, 1=>Rx(pi)), put(n, 1=>shift(0)), put(n, 2=>X), control(n, 2, 1=>Ry(0)), control(n, 2, 1=>shift(pi/2)))
    c = markDiff(c)
    DB_Rx = collect_blocks(QMPS.QDiff, c)[1]
    DB_CS = collect_blocks(QMPS.QDiff, c)[end]
    applyTestReg1 = rand_state(1)
    applyTestReg2 = rand_state(n)
    @test (copy(applyTestReg1) |> DB_Rx) == (apply!(copy(applyTestReg1), DB_Rx))
    @test (copy(applyTestReg1) |> DB_Rx.block) == (apply!(copy(applyTestReg1), DB_Rx.block))
    @test (copy(applyTestReg1) |> DB_Rx) == (copy(applyTestReg1) |> DB_Rx.block)
    @test (copy(applyTestReg2) |> DB_CS) == (apply!(copy(applyTestReg2), DB_CS))
    @test (copy(applyTestReg2) |> DB_CS.block) == (apply!(copy(applyTestReg2), DB_CS.block))
    @test (copy(applyTestReg2) |> DB_CS) == (copy(applyTestReg2) |> DB_CS.block)
    dG = collect_blocks(QMPS.QDiff, c)
    @test dG == [ QMPS.QDiff(Rx(pi)), QMPS.QDiff(control(n, 2, 1=>shift(pi/2))) ]
    @test dG[1].mat == mat(dG[1])
    @test (dG[1])' == adjoint(dG[1]) == QMPS.QDiff( adjoint(dG[1].block) )

    # Testing functions getQdiff!() and getNdiff!().
    ## Define functions.
    function test_diff!(reg::ArrayReg, C::CompositeBlock, op::AbstractBlock, del::Float64, tolerantSetting::Symbol)
        c = deepcopy(C)
        dGates = collect_blocks(QMPS.QDiff, c)
        ldG = length(dGates)
        seed!(seedNum)
        qg = getQdiff!.(()->(copy(reg) |> c), dGates, Ref(op))
        @test qg == [dGates[i].grad for i=1:ldG]
        seed!(seedNum)
        ng = getNdiff!.(()->(copy(reg) |> c), dGates, Ref(op), δ=del)
        seed!(seedNum)
        qg_m2 = getQdiff!(()->(copy(reg) |> c), dGates, op)
        @test qg_m2 == [dGates[i].grad for i=1:ldG]
        seed!(seedNum)
        ng_m2 = getNdiff!(()->(copy(reg) |> c), dGates, op, δ=del)
        tg = zeros(ldG)
        for i = 1:ldG
            dispatch!(-, dGates[i], (del,))
            psi1 = expect(op, copy(reg) |> c) |> mean |> real
            dispatch!(+, dGates[i], (2del,))
            psi2 = expect(op, copy(reg) |> c) |> mean |> real
            dispatch!(-, dGates[i], (del,))
            tg[i] = (psi2 - psi1) / (2del)
        end
        if tolerantSetting == :dc
            @test qg ≈ qg_m2
            @test ng ≈ ng_m2
            @test tg ≈ ng
            @test ng ≈ qg
        end
        if tolerantSetting == :cExtend
            @test qg ≈ qg_m2
            @test ng ≈ ng_m2
            @test tg ≈ ng
            for i=1:ldG @test isapprox(ng[i], qg[i], atol = 2e-4) end
        end
        if tolerantSetting == :cQMPS
             for i=1:ldG @test isapprox(ng[i], ng_m2[i], atol = 2e-3) end
             for i=1:ldG @test isapprox(ng[i], ng_m2[i], atol = 2e-2) end
             for i=1:ldG @test isapprox(tg[i], ng[i], atol = 2e-2) end
             for i=1:ldG @test isapprox(ng[i], qg[i], atol = 1e-1) end
        end
        qg
    end
    function OP(n::Int64)
        seed!(seedNum)
        res = rand_hermitian(2^n) |> matblock
    end
    ## Testing diffs of a differentiable circuit.
    seed!(seedNum)
    del1 = 1e-6
    n1 = 4
    d1 = 3
    reg1 = rand_state(n1)
    c1 = DCbuilder(n1,d1).circuit |> markDiff
    dispatch!(c1, :random)
    test_diff!(reg1, c1, OP(n1), del1, :dc)
    ## Testing diffs of a QMPS-DC and its extended circuit.
    d2 = 2
    n2 = 4
    v2 = 2
    r2 = 1
    seed!(seedNum)
    qmps = MPSC(("DC", d2), n2, v2, r2)
    c2 = deepcopy(qmps.circuit)
    c3 = deepcopy(qmps.cExtend)
    CommonOPs(n::Int64, V::Int64) = [chain(n, [put(n, i=>X) for i=1:V]);
                                     chain(n, [put(n, i=>Y) for i=1:V]);
                                     chain(n, [put(n, i=>Z) for i=1:V])]
    ### Diff of a QMPS when iput state is zero state.
    qg2 = test_diff!(zero_state(v2+r2, nbatch=100000), c2, OP(v2+r2), 5e-1, :cQMPS)
    ### Diff of a MPS extended circuit when iput state is zero state.
    seed!(seedNum)
    test_diff!(rand_state(n2), c3, OP(n2), 1e-2, :cExtend)
    ### Compare the Diffs bewteen the QMPS and corresponding MPS extended circuit.
    for i=1:3
        qg2 = test_diff!(zero_state(v2+r2, nbatch=100000), c2, CommonOPs(v2+r2, v2)[i], 5e-1, :cQMPS)
        qg3 = test_diff!(zero_state(n2), c3, CommonOPs(n2, v2)[i], 1e-2, :cExtend)
        for i=1:length(qg2) @test isapprox(qg2[i], qg3[i], atol = 2e-3) end
    end
end


text:
@testset "fsm_active_close.jl" begin
    base_seq = WrappingInt32(1 << 31)
    DEFAULT_CAPACITY = 64000
    TIMEOUT_DFLT = 1000
    @testset "start in TIME_WAIT, timeout" begin
        conn = TCPConnection()
        #Listen will do nothing
        expect_state(conn, JLSponge.LISTEN)
        tick!(conn, 1)
        expect_state(conn, JLSponge.LISTEN)

        send_syn!(conn, WrappingInt32(0))
        tick!(conn, 1)
        seg = expect_one_seg(conn; ack=true, syn=true, ackno=WrappingInt32(1))
        expect_state(conn, JLSponge.SYN_RCVD)
        send_ack!(conn, WrappingInt32(1), seg.header.seqno + 1)
        tick!(conn, 1)
        expect_no_seg(conn)

        expect_state(conn, JLSponge.ESTABLISHED)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
module Utilities

using Test
using MLJBase
using StableRNGs
import Random
using ComputationalResources

@test MLJBase.finaltypes(Union{Missing,Int}) == [Union{Missing,Int64}]
@test MLJBase.finaltypes(Float64) == [Float64]

abstract type Foo end
struct Bar <: Foo end
struct Baz <: Foo end

@test MLJBase.finaltypes(Foo) == [Bar, Baz]

@testset "flat_values" begin
    t = (X = (x = 1, y = 2), Y = 3)
    @test flat_values(t) == (1, 2, 3)
end

mutable struct M
    a1
    a2
end
mutable struct A1
    a11
    a12
end
mutable struct A2
    a21
end
mutable struct A21
    a211
    a212
end

@testset "recursive getproperty, setproperty!" begin
    m = (a1 = (a11 = 10, a12 = 20), a2 = (a21 = (a211 = 30, a212 = 40),))

    @test MLJBase.recursive_getproperty(m, :(a1.a12)) == 20
    @test MLJBase.recursive_getproperty(m, :a1) == (a11 = 10, a12 = 20)
    @test MLJBase.recursive_getproperty(m, :(a2.a21.a212)) == 40

    m = M(A1(10, 20), A2(A21(30, 40)))
    MLJBase.recursive_setproperty!(m, :(a2.a21.a212), 42)
    @test MLJBase.recursive_getproperty(m, :(a1.a11)) == 10
    @test MLJBase.recursive_getproperty(m, :(a1.a12)) == 20
    @test MLJBase.recursive_getproperty(m, :(a2.a21.a211)) == 30
    @test MLJBase.recursive_getproperty(m, :(a2.a21.a212)) == 42
    @test MLJBase.recursive_getproperty(
        MLJBase.recursive_getproperty(m, :(a2.a21)), :a212) == 42
end

@testset "shuffle rows" begin
    rng = StableRNG(5996661)
    # check dims
    x = randn(rng, 5)
    y = randn(rng, 5, 5)
    z = randn(rng, 5, 5)
    @test MLJBase.check_dimensions(x, y) === nothing
    @test MLJBase.check_dimensions(z, x) === nothing
    @test MLJBase.check_dimensions(y, z) === nothing
    @test_throws DimensionMismatch MLJBase.check_dimensions(x, randn(4))

    x = 1:5 |> collect
    y = 1:5 |> collect
    rng = 555
    Random.seed!(rng)
    perm = Random.randperm(5)
    @test MLJBase.shuffle_rows(x, y; rng=rng) == (x[perm], y[perm])
    y = randn(5,5)
    @test MLJBase.shuffle_rows(x, y; rng=rng) == (x[perm], y[perm,:])
    @test MLJBase.shuffle_rows(z, y; rng=rng) == (z[perm,:], y[perm,:])
    @test MLJBase.shuffle_rows(x, x; rng=rng) == (x[perm], x[perm])
end

@testset "unwind" begin
    iterators = ([1, 2], ["a","b"], ["x", "y", "z"])
    @test unwind(iterators...) ==
        [1  "a"  "x";
         2  "a"  "x";
         1  "b"  "x";
         2  "b"  "x";
         1  "a"  "y";
         2  "a"  "y";
         1  "b"  "y";
         2  "b"  "y";
         1  "a"  "z";
         2  "a"  "z";
         1  "b"  "z";
         2  "b"  "z"]
end

@testset "comp resources" begin
    @test default_resource() == CPU1()
    default_resource(ComputationalResources.CPUProcesses())
    @test default_resource() == ComputationalResources.CPUProcesses()
    default_resource(CPU1())
end

@static if VERSION >= v"1.3.0-DEV.573"
    @testset "Chunks" begin
        nthreads = Threads.nthreads()
        #test for cases with exactly same work as threads
        @test length(MLJBase.chunks(1:nthreads, nthreads)) == nthreads
        #test for cases with more work than threads
        @test length(MLJBase.chunks(1:nthreads + rand(1:nthreads), nthreads)) == nthreads
        #test for cases with less work than threads
        @test length(MLJBase.chunks(1:nthreads-1, nthreads)) == nthreads - 1
    end
end

@testset "available_name" begin
    @test MLJBase.available_name(Utilities, :pizza) == :pizza
    Utilities.eval(:(orange = 5))
    @test MLJBase.available_name(Utilities, :orange) == :orange2
    Utilities.eval(:(orange2 = 6))
    @test MLJBase.available_name(Utilities, :orange) == :orange3
end

end # module
true


text:
import Pkg
# Tell Julia to use the environment specified by Manifest.toml and Project.toml
Pkg.activate(".")
# Install all specified packages
Pkg.instantiate()

# Ensure all packages are precompiled before the parallel workers are started
import JLD
import Flux
import Distributions
import Distributed
import GLM
import StatsModels
import DataFrames
import BoltzmannMachines
import SplitApplyCombine
import Statistics

using Distributed
# Use all processes on this machine:
addprocs()
# Or use multiple machines:
# addprocs([("lenz@imbi2",:auto), ("lenz@imbi10",:auto), ("lenz@imbi12",:auto)];
#        topology = :master_worker)

# Load the required packages and modules on all workers
@everywhere import Pkg
@everywhere Pkg.activate(".")
@everywhere include("./evaluate_odds.jl")

# Define the hyperparameters
seeds = 1:15
modeldefinitions = [
      Dict("model" => "mice");
      map(i -> Dict("model" => "dbm", "epochs" => i), 50:20:2000);
      Dict("model" => "vae", "epochs" => 50:20:2000);
      Dict("model" => "gan", "epochs" => 50:20:2000);
      Dict("model" => "im");
      ]

# Run the experiment
@time result = EvaluateOdds.evaluateoddsratios(modeldefinitions, seeds);

# Write output to tab-separated file
using DelimitedFiles
resultmat = permutedims(hcat(map(collect, result)...))
resultmat = vcat(permutedims(string.(collect(keys(result[1])))), resultmat)
writedlm("results_odds.tsv", resultmat, "\t")

max_stars_repo_name: stefan-m-lenz/dist-gen-comp

max_issues_repo_name: stefan-m-lenz/dist-gen-comp

max_forks_repo_name: stefan-m-lenz/dist-gen-comp


text:
export Simplicial_PPR3_decomposed,
    Simplicial_PPR3_combined,
    SimplicialPROperator,
    grad_and_curl

function hodge_normalization(grad::SpIntMat, curl::SpIntMat,
                             edge_map::Dict{NTuple{2,Int64},Int64})
    d1 = convert(Vector{Float64}, vec(sum(abs.(grad), dims=1)))
    d2 = convert(Vector{Float64}, vec(sum(abs.(curl), dims=1)))
    for ((i, j), ind) in edge_map; d2[ind] += d1[i] + d1[j]; end
    nonzero_inds = findall(d2 .> 0)
    d2[nonzero_inds] = 1.0 ./ d2[nonzero_inds]
    return d2
end

"""
SimplicialPROperator
-------------

Construct the Simplicial PageRank Operator

SimplicialPROperator(grad::SpIntMat, curl::SpIntMat,
                     edge_map::Dict{NTuple{2,Int64},Int64}, α::Float64)

Input parameters:
- grad::SpIntMat: gradient operator (as a matrix)
- curl::SpIntMat: curl operator (as a matrix)
- edge_map::Dict{NTuple{2,Int64}, Int64}: maps an a sorted edge tuple to an index for the matrices
- α::Float64: teleportation parameter

returns simplicial PageRank operator
"""
function SimplicialPROperator(grad::SpIntMat, curl::SpIntMat,
                              edge_map::Dict{NTuple{2,Int64},Int64}, α::Float64)
    G = LinearOperator(convert(SpFltMat, grad))
    C = LinearOperator(convert(SpFltMat, curl))
    Dinv = opDiagonal(hodge_normalization(grad, curl, edge_map))
    L0 = (G * transpose(G) + transpose(C) * C) * Dinv
    β0 = 1 / α - 1
    return (β0 * opEye(size(L0, 1)) + L0)
end

"""
grad_and_curl
-------------

Construct the gradient and curl operators.

grad_and_curl(A::SpIntMat, At::SpIntMat, B::SpIntMat)

Input parameters:
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- At::SpIntMat: the transpose of A
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.

returns tuple (grad, curl, edge_map):
- grad::SpIntMat: gradient operator (as a matrix)
- curl::SpIntMat: curl operator (as a matrix)
- edge_map::Dict{NTuple{2,Int64}, Int64}: maps an a sorted edge tuple to an index for the matrices
"""
function grad_and_curl(A::SpIntMat, At::SpIntMat, B::SpIntMat)
    # Gradient
    edge_map = Dict{NTuple{2, Int64}, Int64}()
    curr_ind = 1
    I, J, V = Int64[], Int64[], Int64[]
    nnodes = size(B, 2)
    for j in 1:nnodes, i in nz_row_inds(B, j)
        if i < j
            edge_map[(i, j)] = curr_ind
            push!(I, curr_ind, curr_ind)
            curr_ind += 1
            push!(J, i, j)
            push!(V, -1, 1)
        end
    end
    nedges = length(edge_map)
    println("$nedges edges...")
    grad = convert(SpIntMat, sparse(I, J, V, nedges, nnodes))

    # Curl
    nthreads = Threads.nthreads()
    Js = Vector{Vector{Int64}}(undef, nthreads)
    for tid in 1:nthreads
        Js[tid] = Vector{Int64}()
    end

    simplex_order = simplex_degree_order(At)
    triangle_order = proj_graph_degree_order(B)
    triangle_id = 1
    shuffled_inds = shuffle(collect(1:nnodes))
    Threads.@threads for ii = 1:length(shuffled_inds)
        i = shuffled_inds[ii]
        for (j, k) in neighbor_pairs(B, triangle_order, i)
            if i == j || i == k || j == k; continue; end
            if B[j, k] > 0 && triangle_closed(A, At, simplex_order, i, j, k)
                # Triangle (i, j, k)
                a, b, c = sort([i, j, k])
                ab, bc, ac = edge_map[(a, b)], edge_map[(b, c)], edge_map[(a, c)]
                push!(Js[Threads.threadid()], ab, bc, ac)
            end
        end
    end

    # combine arrays
    total = sum([length(J) for J in Js])
    cJ = Vector{Int64}(undef, total)
    curr_ind = 1
    for t in 1:nthreads
        size = length(Js[t])
        cJ[curr_ind:(curr_ind + size - 1)] = Js[t][:]
        curr_ind += size
    end
    # triangle indices increment by 1 every third element
    cI = zeros(Int64, length(cJ))
    curr_tri_id = 1
    curr_ind = 1
    while curr_ind <= length(cI)
        for _ in 1:3
            cI[curr_ind] = curr_tri_id
            curr_ind += 1
        end
        curr_tri_id += 1
    end
    ntriangles = curr_tri_id - 1
    # values are (1, 1, -1) for each triangle
    cV = ones(Int64, length(cJ))
    for j = 3:3:length(cV); cV[j] = -1; end
    curl = convert(SpIntMat, sparse(cI, cJ, cV, ntriangles, nedges))

    return (grad, curl, edge_map)
end

"""
Simplicial_PPR3_decomposed
--------------------------

Compute 3-way simplicial personalized PageRank scores for triangles, decomposed
into the harmonic, gradient, and curl components.

Simplicial_PPR3_decomposed(triangles::Vector{NTuple{3,Int64}}, A::SpIntMat, dense_solve::Bool=false, α::Float64=0.85)

Example usage:
```
(scores_comb, scores_curl, scores_grad, scores_harm, S_comb, S_curl, S_grad, S_harm, vec_edge_map) =
    Simplicial_PPR3_decomposed(T, A)
```

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- dense_solve::Bool=false: whether or not to use a dense solver. If the network is tiny, then it might be worth setting this option to true.
- α::Float64=0.85: teleportation parameter for PageRank

returns a tuple (scores_comb, scores_curl, scores_grad, scores_harm, S_comb, S_curl, S_grad, S_harm, vec_edge_map)
- scores_comb::Vector{Float64}: scores for "combined" simplicial personalized PageRank
- scores_curl::Vector{Float64}: scores for curl component
- scores_grad::Vector{Float64}: scores for gradient component
- scores_harm::Vector{Float64}: scores for harmonic component
- S_comb::SpFltMat: a sparse matrix, where S[x, y] is the combined simplicial PPR score for edge x with respect to edge y
- S_curl::SpFltMat: same as S_comb but for the curl component
- S_grad::SpFltMat: same as S_comb but for the gradient component
- S_harm::SpFltMat: same as S_comb but for the harmonic component
- vec_edge_map::Array{Int64,2}: 2 x (# edges) map of indices for S_comb, S_curl, S_grad, S_harm, such that the vec_edge_map[:, i] is the edge for ith index in the matrix
"""
function Simplicial_PPR3_decomposed(triangles::Vector{NTuple{3,Int64}},
                                    A::SpIntMat, dense_solve::Bool=false,
                                    α::Float64=0.85)
    At = convert(SpIntMat, A')
    B = A * At
    B -= Diagonal(B)
    grad, curl, edge_map = grad_and_curl(A, At, B)
    nedges = length(edge_map)

    # find all edges that appear in an open triangle
    in_open_tri = zeros(Int64, nedges)
    I, J = Int64[], Int64[]
    for (i, j, k) in triangles
        a, b, c = sort([i, j, k], alg=InsertionSort)
        ind1, ind2, ind3 = edge_map[(a, b)], edge_map[(a, c)], edge_map[(b, c)]
        in_open_tri[[ind1, ind2, ind3]] .= 1
        push!(I, ind1, ind1, ind2)
        push!(J, ind2, ind3, ind3)
    end
    S = sparse(I, J, ones(Float64, length(I)), nedges, nedges)
    S = make_sparse_ones(S + S')
    S_comb = convert(SpFltMat, S)
    S_grad = copy(S_comb)
    S_curl = copy(S_comb)
    S_harm = copy(S_comb)
    β0 = 1 / α - 1

    if dense_solve
        grad_adj = grad'
        curl_adj = curl'
        Dinv = hodge_normalization(grad, curl, edge_map)
        L1 = (grad * grad_adj + curl_adj * curl) * sparse(Diagonal(Dinv))
        M = Matrix(β0 * sparse(Diagonal(ones(Float64, size(L1, 1)))) + L1)
        println("edge flow solve...")
        Full_S = β0 * inv(M)
        for edge = 1:nedges, i in nz_row_inds(S, edge)
            S_comb[i, edge] = Full_S[i, edge]
        end
        # Gradient component
        println("gradient solve...")
        Sol = grad * (pinv(Matrix(grad)) * Full_S)
        for edge = 1:nedges, i in nz_row_inds(S, edge)
            S_grad[i, edge] = Sol[i, edge]
        end
        # Curl component
        println("curl solve...")
        Sol = curl_adj * (pinv(Matrix(curl_adj)) * Full_S)
        for edge = 1:nedges, i in nz_row_inds(S, edge)
            S_curl[i, edge] = Sol[i, edge]
        end
        # Harmonic component
        for edge = 1:nedges, i in nz_row_inds(S, edge)
            S_harm[i, edge] = S[i, edge] - S_grad[i, edge] - S_curl[i, edge]
        end
    else
        M = SimplicialPROperator(grad, curl, edge_map, α)
        dim = size(M, 1)
        curl_adj = curl'
        edges_in_open_tri = findall(in_open_tri .> 0)
        num_edges_in_open_tri = length(edges_in_open_tri)
        println("$num_edges_in_open_tri edges in open triangles...")

        Threads.@threads for ind = 1:num_edges_in_open_tri
            tid = Threads.threadid()
            if tid == 1
                print(stdout, "$ind of $num_edges_in_open_tri \r")
                flush(stdout)
            end
            edge = edges_in_open_tri[ind]
            b = zeros(Float64, dim)
            b[edge] = 1.0
            sol_comb = dqgmres(M, b, atol=1e-4, rtol=1e-4)[1] * β0
            # split into gradient, harmonic, and curl components
            sol_grad = grad * lsqr(grad, sol_comb, atol=1e-3, btol=1e-3)
            sol_curl = curl_adj * lsqr(curl_adj, sol_comb, atol=1e-3, btol=1e-3)
            sol_harm = sol_comb - sol_grad - sol_curl
            for i in nz_row_inds(S, edge)
                S_comb[i, edge] = sol_comb[i]
                S_grad[i, edge] = sol_grad[i]
                S_curl[i, edge] = sol_curl[i]
                S_harm[i, edge] = sol_harm[i]
            end
        end
    end

    ntriangles = length(triangles)
    scores_comb = zeros(Float64, ntriangles)
    scores_grad = zeros(Float64, ntriangles)
    scores_curl = zeros(Float64, ntriangles)
    scores_harm = zeros(Float64, ntriangles)
    function mat_score(X::SpFltMat, r1::Int64, r2::Int64, r3::Int64)
        s  = abs(X[r1, r2])
        s += abs(X[r2, r1])
        s += abs(X[r1, r3])
        s += abs(X[r3, r1])
        s += abs(X[r2, r3])
        s += abs(X[r3, r2])
        return s
    end
    Threads.@threads for ind = 1:ntriangles
        i, j, k = triangles[ind]
        a, b, c = sort([i, j, k], alg=InsertionSort)
        r1 = edge_map[(a, b)]
        r2 = edge_map[(a, c)]
        r3 = edge_map[(b, c)]
        scores_comb[ind] = mat_score(S_comb, r1, r2, r3)
        scores_grad[ind] = mat_score(S_grad, r1, r2, r3)
        scores_curl[ind] = mat_score(S_curl, r1, r2, r3)
        scores_harm[ind] = mat_score(S_harm, r1, r2, r3)
    end

    vec_edge_map = zeros(Int64, 2, length(edge_map))
    for ((i, j), ind) in edge_map
        vec_edge_map[1, ind] = i
        vec_edge_map[2, ind] = j
    end

    return (scores_comb, scores_curl, scores_grad, scores_harm,
            S_comb,      S_curl,      S_grad,      S_harm,
            vec_edge_map)
end

"""
Simplicial_PPR3_combined
--------------------------

Compute the "combined" 3-way simplicial personalized PageRank scores for
triangles. This is less computationally expensive that the related function
Simplicial_PPR3_decomposed(), which decomposes the scores into the curl,
gradient, and harmonic components.

function Simplicial_PPR3_combined(triangles::Vector{NTuple{3,Int64}}, A::SpIntMat,
                                  At::SpIntMat, B::SpIntMat, α::Float64=0.85)

Example usage:
```
(scores_comb, S_comb, vec_edge_map) = Simplicial_PPR3_combined(T, A)
```

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- A::SpIntMat: (# nodes) x (# simplices) adjacency matrix
- α::Float64=0.85: teleportation parameter for PageRank

returns a tuple (scores_comb, S_comb, vec_edge_map)
- scores_comb::Vector{Float64}: scores for "combined" simplicial personalized PageRank
- S_comb::SpFltMat: a sparse matrix, where S[x, y] is the combined simplicial PPR score for edge x with respect to edge y
- vec_edge_map::Array{Int64,2}: 2 x (# edges) map of indices for S_comb such that the vec_edge_map[:, i] is the edge for ith index in the matrix
"""
function Simplicial_PPR3_combined(triangles::Vector{NTuple{3,Int64}},
                                  A::SpIntMat, α::Float64=0.85)
    At = convert(SpIntMat, A')
    B = A * At
    B -= Diagonal(B)
    grad, curl, edge_map = grad_and_curl(A, At, B)
    nedges = length(edge_map)

    # find all edges that appear in an open triangle
    in_open_tri = zeros(Int64, nedges)
    I, J = Int64[], Int64[]
    for (i, j, k) in triangles
        a, b, c = sort([i, j, k], alg=InsertionSort)
        ind1, ind2, ind3 = edge_map[(a, b)], edge_map[(a, c)], edge_map[(b, c)]
        in_open_tri[[ind1, ind2, ind3]] .= 1
        push!(I, ind1, ind1, ind2)
        push!(J, ind2, ind3, ind3)
    end
    S = sparse(I, J, ones(Float64, length(I)), nedges, nedges)
    S = make_sparse_ones(S + S')
    S_comb = convert(SpFltMat, S)

    M = SimplicialPROperator(grad, curl, edge_map, α)
    dim = size(M, 1)
    curl_adj = curl'
    edges_in_open_tri = findall(in_open_tri .> 0)
    num_edges_in_open_tri = length(edges_in_open_tri)
    println("$num_edges_in_open_tri edges in open triangles...")

    Threads.@threads for ind = 1:num_edges_in_open_tri
        tid = Threads.threadid()
        if tid == 1
            print(stdout, "$ind of $num_edges_in_open_tri \r")
            flush(stdout)
        end
        edge = edges_in_open_tri[ind]
        b = zeros(Float64, dim)
        b[edge] = 1.0
        sol_comb = dqgmres(M, b, atol=1e-4, rtol=1e-4)[1] * (1 - α)
        for i in nz_row_inds(S, edge)
            S_comb[i, edge] = sol_comb[i]
        end
    end

    ntriangles = length(triangles)
    scores_comb = zeros(Float64, ntriangles)
    Threads.@threads for ind = 1:ntriangles
        i, j, k = triangles[ind]
        a, b, c = sort([i, j, k], alg=InsertionSort)
        r1 = edge_map[(a, b)]
        r2 = edge_map[(a, c)]
        r3 = edge_map[(b, c)]
        score  = abs(S_comb[r1, r2])
        score += abs(S_comb[r2, r1])
        score += abs(S_comb[r1, r3])
        score += abs(S_comb[r3, r1])
        score += abs(S_comb[r2, r3])
        score += abs(S_comb[r3, r2])
        scores_comb[ind] = score
    end

    vec_edge_map = zeros(Int64, 2, length(edge_map))
    for ((i, j), ind) in edge_map
        vec_edge_map[1, ind] = i
        vec_edge_map[2, ind] = j
    end

    return (scores_comb, S_comb, vec_edge_map)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: synthetics
using AWS.Compat
using AWS.UUIDs

"""
    create_canary(artifact_s3_location, code, execution_role_arn, name, runtime_version, schedule)
    create_canary(artifact_s3_location, code, execution_role_arn, name, runtime_version, schedule, params::Dict{String,<:Any})

Creates a canary. Canaries are scripts that monitor your endpoints and APIs from the
outside-in. Canaries help you check the availability and latency of your web services and
troubleshoot anomalies by investigating load time data, screenshots of the UI, logs, and
metrics. You can set up a canary to run continuously or just once.  Do not use CreateCanary
to modify an existing canary. Use UpdateCanary instead. To create canaries, you must have
the CloudWatchSyntheticsFullAccess policy. If you are creating a new IAM role for the
canary, you also need the the iam:CreateRole, iam:CreatePolicy and iam:AttachRolePolicy
permissions. For more information, see Necessary Roles and Permissions. Do not include
secrets or proprietary information in your canary names. The canary name makes up part of
the Amazon Resource Name (ARN) for the canary, and the ARN is included in outbound calls
over the internet. For more information, see Security Considerations for Synthetics
Canaries.

# Arguments
- `artifact_s3_location`: The location in Amazon S3 where Synthetics stores artifacts from
  the test runs of this canary. Artifacts include the log file, screenshots, and HAR files.
  The name of the S3 bucket can't include a period (.).
- `code`: A structure that includes the entry point from which the canary should start
  running your script. If the script is stored in an S3 bucket, the bucket name, key, and
  version are also included.
- `execution_role_arn`: The ARN of the IAM role to be used to run the canary. This role
  must already exist, and must include lambda.amazonaws.com as a principal in the trust
  policy. The role must also have the following permissions:    s3:PutObject
  s3:GetBucketLocation     s3:ListAllMyBuckets     cloudwatch:PutMetricData
  logs:CreateLogGroup     logs:CreateLogStream     logs:PutLogEvents
- `name`: The name for this canary. Be sure to give it a descriptive name that
  distinguishes it from other canaries in your account. Do not include secrets or proprietary
  information in your canary names. The canary name makes up part of the canary ARN, and the
  ARN is included in outbound calls over the internet. For more information, see Security
  Considerations for Synthetics Canaries.
- `runtime_version`: Specifies the runtime version to use for the canary. For a list of
  valid runtime versions and more information about runtime versions, see  Canary Runtime
  Versions.
- `schedule`: A structure that contains information about how often the canary is to run
  and when these test runs are to stop.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ArtifactConfig"`: A structure that contains the configuration for canary artifacts,
  including the encryption-at-rest settings for artifacts that the canary uploads to Amazon
  S3.
- `"FailureRetentionPeriodInDays"`: The number of days to retain data about failed runs of
  this canary. If you omit this field, the default of 31 days is used. The valid range is 1
  to 455 days.
- `"RunConfig"`: A structure that contains the configuration for individual canary runs,
  such as timeout value.
- `"SuccessRetentionPeriodInDays"`: The number of days to retain data about successful runs
  of this canary. If you omit this field, the default of 31 days is used. The valid range is
  1 to 455 days.
- `"Tags"`: A list of key-value pairs to associate with the canary. You can associate as
  many as 50 tags with a canary. Tags can help you organize and categorize your resources.
  You can also use them to scope user permissions, by granting a user permission to access or
  change only the resources that have certain tag values.
- `"VpcConfig"`: If this canary is to test an endpoint in a VPC, this structure contains
  information about the subnet and security groups of the VPC endpoint. For more information,
  see  Running a Canary in a VPC.
"""
function create_canary(
    ArtifactS3Location,
    Code,
    ExecutionRoleArn,
    Name,
    RuntimeVersion,
    Schedule;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return synthetics(
        "POST",
        "/canary",
        Dict{String,Any}(
            "ArtifactS3Location" => ArtifactS3Location,
            "Code" => Code,
            "ExecutionRoleArn" => ExecutionRoleArn,
            "Name" => Name,
            "RuntimeVersion" => RuntimeVersion,
            "Schedule" => Schedule,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_canary(
    ArtifactS3Location,
    Code,
    ExecutionRoleArn,
    Name,
    RuntimeVersion,
    Schedule,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return synthetics(
        "POST",
        "/canary",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ArtifactS3Location" => ArtifactS3Location,
                    "Code" => Code,
                    "ExecutionRoleArn" => ExecutionRoleArn,
                    "Name" => Name,
                    "RuntimeVersion" => RuntimeVersion,
                    "Schedule" => Schedule,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_canary(name)
    delete_canary(name, params::Dict{String,<:Any})

Permanently deletes the specified canary. When you delete a canary, resources used and
created by the canary are not automatically deleted. After you delete a canary that you do
not intend to use again, you should also delete the following:   The Lambda functions and
layers used by this canary. These have the prefix cwsyn-MyCanaryName .   The CloudWatch
alarms created for this canary. These alarms have a name of
Synthetics-SharpDrop-Alarm-MyCanaryName .   Amazon S3 objects and buckets, such as the
canary's artifact location.   IAM roles created for the canary. If they were created in the
console, these roles have the name  role/service-role/CloudWatchSyntheticsRole-MyCanaryName
.   CloudWatch Logs log groups created for the canary. These logs groups have the name
/aws/lambda/cwsyn-MyCanaryName .    Before you delete a canary, you might want to use
GetCanary to display the information about this canary. Make note of the information
returned by this operation so that you can delete these resources after you delete the
canary.

# Arguments
- `name`: The name of the canary that you want to delete. To find the names of your
  canaries, use DescribeCanaries.

"""
function delete_canary(name; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "DELETE", "/canary/$(name)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function delete_canary(
    name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "DELETE",
        "/canary/$(name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_canaries()
    describe_canaries(params::Dict{String,<:Any})

This operation returns a list of the canaries in your account, along with full details
about each canary. This operation does not have resource-level authorization, so if a user
is able to use DescribeCanaries, the user can see all of the canaries in the account. A
deny policy can only be used to restrict access to all canaries. It cannot be used on
specific resources.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: Specify this parameter to limit how many canaries are returned each time
  you use the DescribeCanaries operation. If you omit this parameter, the default of 100 is
  used.
- `"NextToken"`: A token that indicates that there is more data available. You can use this
  token in a subsequent operation to retrieve the next set of results.
"""
function describe_canaries(; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST", "/canaries"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_canaries(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "POST", "/canaries", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    describe_canaries_last_run()
    describe_canaries_last_run(params::Dict{String,<:Any})

Use this operation to see information from the most recent run of each canary that you have
created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: Specify this parameter to limit how many runs are returned each time you
  use the DescribeLastRun operation. If you omit this parameter, the default of 100 is used.
- `"NextToken"`: A token that indicates that there is more data available. You can use this
  token in a subsequent DescribeCanaries operation to retrieve the next set of results.
"""
function describe_canaries_last_run(; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST", "/canaries/last-run"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_canaries_last_run(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "POST",
        "/canaries/last-run",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_runtime_versions()
    describe_runtime_versions(params::Dict{String,<:Any})

Returns a list of Synthetics canary runtime versions. For more information, see  Canary
Runtime Versions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: Specify this parameter to limit how many runs are returned each time you
  use the DescribeRuntimeVersions operation. If you omit this parameter, the default of 100
  is used.
- `"NextToken"`: A token that indicates that there is more data available. You can use this
  token in a subsequent DescribeRuntimeVersions operation to retrieve the next set of results.
"""
function describe_runtime_versions(; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST", "/runtime-versions"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_runtime_versions(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "POST",
        "/runtime-versions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_canary(name)
    get_canary(name, params::Dict{String,<:Any})

Retrieves complete information about one canary. You must specify the name of the canary
that you want. To get a list of canaries and their names, use DescribeCanaries.

# Arguments
- `name`: The name of the canary that you want details for.

"""
function get_canary(name; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "GET", "/canary/$(name)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_canary(
    name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "GET",
        "/canary/$(name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_canary_runs(name)
    get_canary_runs(name, params::Dict{String,<:Any})

Retrieves a list of runs for a specified canary.

# Arguments
- `name`: The name of the canary that you want to see runs for.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: Specify this parameter to limit how many runs are returned each time you
  use the GetCanaryRuns operation. If you omit this parameter, the default of 100 is used.
- `"NextToken"`: A token that indicates that there is more data available. You can use this
  token in a subsequent GetCanaryRuns operation to retrieve the next set of results.
"""
function get_canary_runs(name; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST",
        "/canary/$(name)/runs";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_canary_runs(
    name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "POST",
        "/canary/$(name)/runs",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Displays the tags associated with a canary.

# Arguments
- `resource_arn`: The ARN of the canary that you want to view tags for. The ARN format of a
  canary is arn:aws:synthetics:Region:account-id:canary:canary-name .

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return synthetics(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_canary(name)
    start_canary(name, params::Dict{String,<:Any})

Use this operation to run a canary that has already been created. The frequency of the
canary runs is determined by the value of the canary's Schedule. To see a canary's
schedule, use GetCanary.

# Arguments
- `name`: The name of the canary that you want to run. To find canary names, use
  DescribeCanaries.

"""
function start_canary(name; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST",
        "/canary/$(name)/start";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_canary(
    name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "POST",
        "/canary/$(name)/start",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    stop_canary(name)
    stop_canary(name, params::Dict{String,<:Any})

Stops the canary to prevent all future runs. If the canary is currently running, Synthetics
stops waiting for the current run of the specified canary to complete. The run that is in
progress completes on its own, publishes metrics, and uploads artifacts, but it is not
recorded in Synthetics as a completed run. You can use StartCanary to start it running
again with the canary’s current schedule at any point in the future.

# Arguments
- `name`: The name of the canary that you want to stop. To find the names of your canaries,
  use DescribeCanaries.

"""
function stop_canary(name; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST",
        "/canary/$(name)/stop";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function stop_canary(
    name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "POST",
        "/canary/$(name)/stop",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(tags, resource_arn)
    tag_resource(tags, resource_arn, params::Dict{String,<:Any})

Assigns one or more tags (key-value pairs) to the specified canary.  Tags can help you
organize and categorize your resources. You can also use them to scope user permissions, by
granting a user permission to access or change only resources with certain tag values. Tags
don't have any semantic meaning to Amazon Web Services and are interpreted strictly as
strings of characters. You can use the TagResource action with a canary that already has
tags. If you specify a new tag key for the alarm, this tag is appended to the list of tags
associated with the alarm. If you specify a tag key that is already associated with the
alarm, the new tag value that you specify replaces the previous value for that tag. You can
associate as many as 50 tags with a canary.

# Arguments
- `tags`: The list of key-value pairs to associate with the canary.
- `resource_arn`: The ARN of the canary that you're adding tags to. The ARN format of a
  canary is arn:aws:synthetics:Region:account-id:canary:canary-name .

"""
function tag_resource(Tags, resourceArn; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    Tags,
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return synthetics(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Tags" => Tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes one or more tags from the specified canary.

# Arguments
- `resource_arn`: The ARN of the canary that you're removing tags from. The ARN format of a
  canary is arn:aws:synthetics:Region:account-id:canary:canary-name .
- `tag_keys`: The list of tag keys to remove from the resource.

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return synthetics(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_canary(name)
    update_canary(name, params::Dict{String,<:Any})

Use this operation to change the settings of a canary that has already been created. You
can't use this operation to update the tags of an existing canary. To change the tags of an
existing canary, use TagResource.

# Arguments
- `name`: The name of the canary that you want to update. To find the names of your
  canaries, use DescribeCanaries. You cannot change the name of a canary that has already
  been created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ArtifactConfig"`: A structure that contains the configuration for canary artifacts,
  including the encryption-at-rest settings for artifacts that the canary uploads to Amazon
  S3.
- `"ArtifactS3Location"`: The location in Amazon S3 where Synthetics stores artifacts from
  the test runs of this canary. Artifacts include the log file, screenshots, and HAR files.
  The name of the S3 bucket can't include a period (.).
- `"Code"`: A structure that includes the entry point from which the canary should start
  running your script. If the script is stored in an S3 bucket, the bucket name, key, and
  version are also included.
- `"ExecutionRoleArn"`: The ARN of the IAM role to be used to run the canary. This role
  must already exist, and must include lambda.amazonaws.com as a principal in the trust
  policy. The role must also have the following permissions:    s3:PutObject
  s3:GetBucketLocation     s3:ListAllMyBuckets     cloudwatch:PutMetricData
  logs:CreateLogGroup     logs:CreateLogStream     logs:CreateLogStream
- `"FailureRetentionPeriodInDays"`: The number of days to retain data about failed runs of
  this canary.
- `"RunConfig"`: A structure that contains the timeout value that is used for each
  individual run of the canary.
- `"RuntimeVersion"`: Specifies the runtime version to use for the canary. For a list of
  valid runtime versions and for more information about runtime versions, see  Canary Runtime
  Versions.
- `"Schedule"`: A structure that contains information about how often the canary is to run,
  and when these runs are to stop.
- `"SuccessRetentionPeriodInDays"`: The number of days to retain data about successful runs
  of this canary.
- `"VisualReference"`: Defines the screenshots to use as the baseline for comparisons
  during visual monitoring comparisons during future runs of this canary. If you omit this
  parameter, no changes are made to any baseline screenshots that the canary might be using
  already. Visual monitoring is supported only on canaries running the syn-puppeteer-node-3.2
  runtime or later. For more information, see  Visual monitoring and  Visual monitoring
  blueprint
- `"VpcConfig"`: If this canary is to test an endpoint in a VPC, this structure contains
  information about the subnet and security groups of the VPC endpoint. For more information,
  see  Running a Canary in a VPC.
"""
function update_canary(name; aws_config::AbstractAWSConfig=global_aws_config())
    return synthetics(
        "PATCH", "/canary/$(name)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function update_canary(
    name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return synthetics(
        "PATCH",
        "/canary/$(name)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
# Execution control
# copied from CUDAdrv/src/execution.jl

export ROCDim, ROCModule, ROCFunction, roccall

mutable struct ROCModule
    #data::String
    data::Vector{UInt8}
    options::Dict{Any,Any}
end
mutable struct ROCFunction
    mod::ROCModule
    entry::String
end

"""
    ROCDim3(x)

    ROCDim3((x,))
    ROCDim3((x, y))
    ROCDim3((x, y, x))

A type used to specify dimensions, consisting of 3 integers for respectively
the `x`, `y` and `z` dimension. Unspecified dimensions default to `1`.

Often accepted as argument through the `ROCDim` type alias, eg. in the case of
[`roccall`](@ref) or [`launch`](@ref), allowing to pass dimensions as a plain
integer or a tuple without having to construct an explicit `ROCDim3` object.
"""
struct ROCDim3
    x::Cuint
    y::Cuint
    z::Cuint
end

ROCDim3(dims::Integer)             = ROCDim3(dims,    Cuint(1), Cuint(1))
ROCDim3(dims::NTuple{1,<:Integer}) = ROCDim3(dims[1], Cuint(1), Cuint(1))
ROCDim3(dims::NTuple{2,<:Integer}) = ROCDim3(dims[1], dims[2],  Cuint(1))
ROCDim3(dims::NTuple{3,<:Integer}) = ROCDim3(dims[1], dims[2],  dims[3])

# Type alias for conveniently specifying the dimensions
# (e.g. `(len, 2)` instead of `ROCDim3((len, 2))`)
const ROCDim = Union{Integer,
                    Tuple{Integer},
                    Tuple{Integer, Integer},
                    Tuple{Integer, Integer, Integer}}

function Base.getindex(dims::ROCDim3, idx::Int)
    return idx == 1 ? dims.x :
           idx == 2 ? dims.y :
           idx == 3 ? dims.z :
           error("Invalid dimension: $idx")
end

"""
    launch(queue::HSAQueue, signal::HSASignal, f::ROCFunction,
           groupsize::ROCDim, gridsize::ROCDim, args...)

Low-level call to launch a ROC function `f` on the GPU, using `groupsize` and
`gridsize` as respectively the grid and block configuration. The kernel is
launched on queue `queue` and is waited on by signal `signal`.

Arguments to a kernel should either be bitstype, in which case they will be
copied to the internal kernel parameter buffer, or a pointer to device memory.

This is a low-level call, prefer to use [`roccall`](@ref) instead.
"""
@inline function launch(queue::HSAQueue, signal::HSASignal, f::ROCFunction,
                        groupsize::ROCDim, gridsize::ROCDim, args...)
    groupsize = ROCDim3(groupsize)
    gridsize = ROCDim3(gridsize)
    (groupsize.x>0 && groupsize.y>0 && groupsize.z>0) ||
        throw(ArgumentError("Group dimensions should be non-null"))
    (gridsize.x>0 && gridsize.y>0 && gridsize.z>0) ||
        throw(ArgumentError("Grid dimensions should be non-null"))

    _launch(queue, signal, f, groupsize, gridsize, args...)
end

# we need a generated function to get an args array,
# without having to inspect the types at runtime
@generated function _launch(queue::HSAQueue, signal::HSASignal, f::ROCFunction,
                            groupsize::ROCDim3, gridsize::ROCDim3,
                            args::NTuple{N,Any}) where N
    all(isbitstype, args.parameters) ||
        throw(ArgumentError("Arguments to kernel should be bitstype."))

    ex = Expr(:block)
    push!(ex.args, :(Base.@_inline_meta))

    # If f has N parameters, then kernelParams needs to be an array of N
    # pointers.  Each of kernelParams[0] through kernelParams[N-1] must point
    # to a region of memory from which the actual kernel parameter will be
    # copied.

    # put arguments in Ref boxes so that we can get a pointers to them
    arg_refs = Vector{Symbol}(undef, N)
    for i in 1:N
        arg_refs[i] = gensym()
        push!(ex.args, :($(arg_refs[i]) = Base.RefValue(args[$i])))
    end

    # generate an array with pointers
    arg_ptrs = [:(Base.unsafe_convert(Ptr{Cvoid}, $(arg_refs[i]))) for i in 1:N]

    append!(ex.args, (quote
        GC.@preserve $(arg_refs...) begin
            kernelParams = [$(arg_ptrs...)]

            # link with ld.lld
            ld_path = HSARuntime.ld_lld_path
            @assert ld_path != "" "ld.lld was not found; cannot link kernel"
            # TODO: Do this more idiomatically
            io = open("/tmp/amdgpu-dump.o", "w")
            write(io, f.mod.data)
            close(io)
            run(`$ld_path -shared -o /tmp/amdgpu.exe /tmp/amdgpu-dump.o`)
            io = open("/tmp/amdgpu.exe", "r")
            data = read(io)
            close(io)

            # generate executable and kernel instance
            exe = HSAExecutable(queue.agent, data, f.entry)
            kern = HSAKernelInstance(queue.agent, exe, f.entry, args)
            HSARuntime.launch!(queue, kern, signal;
                workgroup_size=groupsize, grid_size=gridsize)
        end
    end).args)

    return ex
end

"""
    roccall(queue::HSAQueue, signal::HSASignal, f::ROCFunction, types, values...;
             groupsize::ROCDim, gridsize::ROCDim)

`ccall`-like interface for launching a ROC function `f` on a GPU.

For example:

    vadd = ROCFunction(md, "vadd")
    a = rand(Float32, 10)
    b = rand(Float32, 10)
    ad = Mem.upload(a)
    bd = Mem.upload(b)
    c = zeros(Float32, 10)
    cd = Mem.alloc(c)

    roccall(vadd, (Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat}), ad, bd, cd;
             gridsize=10)
    Mem.download!(c, cd)

The `groupsize` and `gridsize` arguments control the launch configuration, and should both
consist of either an integer, or a tuple of 1 to 3 integers (omitted dimensions default to
1). The `types` argument can contain both a tuple of types, and a tuple type, the latter
being slightly faster.
"""
roccall

@inline function roccall(queue::HSAQueue, signal::HSASignal, f::ROCFunction, types::NTuple{N,DataType}, values::Vararg{Any,N};
                          kwargs...) where N
    # this cannot be inferred properly (because types only contains `DataType`s),
    # which results in the call `@generated _roccall` getting expanded upon first use
    _roccall(queue, signal, f, Tuple{types...}, values; kwargs...)
end

@inline function roccall(queue::HSAQueue, signal::HSASignal, f::ROCFunction, tt::Type, values::Vararg{Any,N};
                          kwargs...) where N
    # in this case, the type of `tt` is `Tuple{<:DataType,...}`,
    # which means the generated function can be expanded earlier
    _roccall(queue, signal, f, tt, values; kwargs...)
end

# we need a generated function to get a tuple of converted arguments (using unsafe_convert),
# without having to inspect the types at runtime
@generated function _roccall(queue::HSAQueue, signal::HSASignal, f::ROCFunction, tt::Type, args::NTuple{N,Any};
                              groupsize::ROCDim=1, gridsize::ROCDim=1) where N
    types = tt.parameters[1].parameters     # the type of `tt` is Type{Tuple{<:DataType...}}

    ex = Expr(:block)
    push!(ex.args, :(Base.@_inline_meta))

    # convert the argument values to match the kernel's signature (specified by the user)
    # (this mimics `lower-ccall` in julia-syntax.scm)
    converted_args = Vector{Symbol}(undef, N)
    arg_ptrs = Vector{Symbol}(undef, N)
    for i in 1:N
        converted_args[i] = gensym()
        arg_ptrs[i] = gensym()
        push!(ex.args, :($(converted_args[i]) = Base.cconvert($(types[i]), args[$i])))
        push!(ex.args, :($(arg_ptrs[i]) = Base.unsafe_convert($(types[i]), $(converted_args[i]))))
    end

    append!(ex.args, (quote
        GC.@preserve $(converted_args...) begin
            launch(queue, signal, f, groupsize, gridsize, ($(arg_ptrs...),))
        end
    end).args)

    return ex
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
Base.zeros(::Type{ANode}, sz::Int...) = ANode(:zeros, (), sz)
Base.zeros(::Type{ANode}, sz::Tuple{Vararg{Int}}) = ANode(:zeros, (), sz)
do_forward!(::ANode{:zeros}, out) = fill!(out, zero(eltype(out)))
elementwise(::ANode{:zeros}) = true
@llvm llvm_val(n::ANode{:zeros}, ::ModuleRef, b::BuilderRef, t::Type) = llvm_const(t, 0)

Base.ones(::Type{ANode}, sz::Int...) = ANode(:ones, (), sz)
Base.ones(::Type{ANode}, sz::Tuple{Vararg{Int}}) = ANode(:ones, (), sz)
do_forward!(::ANode{:ones}, out) = fill!(out, one(eltype(out)))
elementwise(::ANode{:ones}) = true
@llvm llvm_val(n::ANode{:ones}, ::ModuleRef, b::BuilderRef, t::Type) = llvm_const(t, 1)

Base.fill!(a::ANode, value) = ANode(:fill!, (a,), size(a), convert(eltype(a), value))
mutates(::ANode{:fill!}) = true
do_forward!(a::ANode{:fill!}, out) = fill!(out, a.arg)
gradient_node!(a::ANode{:fill!}, wrt::ANode, b) = zeros(size(wrt))
#elementwise(::ANode{:fill!}) = true
#@llvm llvm_val(n::ANode{:fill!}, ::ModuleRef, b::BuilderRef, t::Type) = llvm_const(t, n.arg)

Base.rand!(a::ANode) = ANode(:rand!, (a,), size(a))
mutates(::ANode{:rand!}) = true
do_forward!(a::ANode{:rand!}, out) = rand!(out)
gradient_node!(a::ANode{:rand!}, wrt::ANode, b) = zeros(size(wrt))

Base.randn!(a::ANode) = ANode(:randn!, (a,), size(a))
mutates(::ANode{:randn!}) = true
do_forward!(a::ANode{:randn!}, out) = randn!(out)
gradient_node!(a::ANode{:randn!}, wrt::ANode, b) = zeros(size(wrt))


text:
# This file is not an active part of the package. This is the code
# that uses the Clang.jl package to wrap Libui using the headers.

# Find all headers
# To run the script from Julia console:
push!(Libdl.DL_LOAD_PATH, "/usr/local/Cellar/llvm/6.0.1/lib")
using Clang.wrap_c
const wdir = dirname(@__FILE__)
const pkgbasedir = joinpath(wdir, "..")

# `outpath` specifies, where the julian wrappers would be generated.
# If the generated .jl files are ok, they have to be copied to the "src" folder
# overwriting the old ones
const outpath = joinpath(pkgbasedir, "new")
rm(outpath, recursive = true, force = true)
mkpath(outpath)

const incpath = joinpath(pkgbasedir, "deps", "src")
if !isdir(incpath)
  error("Run Pkg.build(\"Libui\") before trying to wrap C headers.")
end

info("Scanning Libui headers in $incpath...")
const Libui_header_file = joinpath(incpath, "ui.h")

const clang_path = "/usr/local/Cellar/llvm/6.0.1/lib/clang/6.0.1/" # change to your clang location
const includes = [
	joinpath(clang_path, "include"),
	incpath,
]

function find_outfile(s)
	joinpath(outpath, string(first(splitext(basename(s))), ".jl"))
end
find_libfile(s) = "libui"

const context = wrap_c.init(
	headers = [Libui_header_file,],
	common_file = joinpath(outpath, "ui_h.jl"),
	clang_args = [
		"-D", "__STDC_LIMIT_MACROS",
		"-D", "__STDC_CONSTANT_MACROS",
		# "-v"
	],
	# clang_diagnostics = true,
	header_library = find_libfile,
	header_outputfile = find_outfile,
	clang_includes = includes,
)

info("Generating .jl wrappers for Libui in $outpath...")
run(context)
info("Done generating .jl wrappers for Libui in $outpath")


text:
# Header type with restricted choice of keys and value types

module Headers

export Header, headerval

abstract type Header
	#=
	struct
		D<:Dict{Symbol,Any}
	end
	=#
end

@inline Base.getproperty(header::Header, name::Symbol) =
		getproperty(header, Val(name))

@inline Base.getproperty(header::Header, ::Val{:D}) = getfield(header, :D)

@inline Base.getproperty(header::Header, ::Val{name}) where {name} =
		getindex(header, name)

@inline Base.setproperty!(header::Header, name::Symbol, v) =
		setindex!(header, v, name)

Base.propertynames(header::Header, private::Bool = false) =
		private ? (:D, keys(header)...) : (keys(header)...,)

Base.show(io::IO, header::Header) = print(io, "$(typeof(header))()")

function Base.show(io::IO, ::MIME"text/plain", header::Header)
	print(io, "$(length(header))-entry $(typeof(header))")
	if length(header) > 0
		print(io, ":")
		for key in keys(header)
			print(io, "\n $(key)")
		end
	end
end

Base.getindex(header::Header, key::Symbol) = header.D[key]

Base.setindex!(header::Header, v, key::Symbol) =
		(header.D[key] = headerval(header, Val(key), v))

Base.haskey(header::Header, key::Symbol) = haskey(header.D, key)

Base.getkey(header::Header, key::Symbol, default) =
		getkey(header.D, key, default)

Base.get(header::Header, key::Symbol, default) = get(header.D, key, default)

Base.get(f, header::Header, key::Symbol) = get(f, header.D, key)

function Base.get!(f, header::Header, key::Symbol)
	get!(header.D, key) do
		headerval(header, Val(key), f())
	end
end

Base.get!(header::Header, key::Symbol, default) =
		(get!(header, key) do; default; end)

function Base.delete!(header::Header, key::Symbol)
	delete!(header.D, key)
	header
end

function Base.empty!(header::Header)
	empty!(header.D)
	header
end

Base.iterate(header::Header) = iterate(header.D)

Base.iterate(header::Header, state) = iterate(header.D, state)

Base.IteratorSize(header::Header) = Base.IteratorSize(header.D)

Base.IteratorEltype(header::Header) = Base.IteratorEltype(header.D)

Base.eltype(header::Header) = Pair{Symbol,Any}

Base.length(header::Header) = length(header.D)

Base.keys(header::Header) = keys(header.D)

Base.values(header::Header) = values(header.D)

Base.pairs(header::Header) = pairs(header.D)

Base.:(==)(header1::Header, header2::Header) = (header1 === header2)

function Base.merge!(header::Header, src)
	for (key, value) in pairs(src)
		header[key] = value
	end
	header
end

function headerval end

end # module


text:
using VT100
import VT100: Line, Cursor, Cell
import Base: getindex, convert, show
using VT100.Flags
using VT100.Attributes

import VT100: get_image_cell, pos_for_image_cell

using Colors
using FixedPointNumbers

typealias RGB8 RGB{UFixed8}

import Base: ==, hash

immutable CellLoc
    row::Int
    col::Int
end
convert(::Type{CellLoc},t::Tuple{Int,Int}) = CellLoc(t...)
==(x::CellLoc,y::CellLoc) = x.row == y.row && x.col == y.col

immutable CellRect
    height::Int
    width::Int
end
convert(::Type{CellRect},t::Tuple{Int,Int}) = CellRect(t...)
==(x::CellRect,y::CellRect) = x.height == y.height && x.width == y.width

immutable CellBox
    topleft::CellLoc
    size::CellRect
end
==(x::CellBox, y::CellBox) = x.topleft == y.topleft && x.size == y.size

typealias LocOrRect Union{CellLoc,CellRect}

==(x::LocOrRect,y::Tuple{Int,Int}) = x == (typeof(x))(y...)
==(x::Tuple{Int,Int},y::LocOrRect) = (typeof(y))(x...) == y

getindex(x::LocOrRect,i) = getfield(x,i)

# Used to minimize the number of characters to be written.
# The written characters themselves still have to be buffered,
# which is done by an IOBuffer in the `render` method.
type DoubleBufferedTerminalScreen <: Screen
    # If this is empty, the screen state is unkown and should be cleared first
    have::Vector{Line}
    want::Vector{Line}
    wantcursor::Cursor
    # Technically these are available by looking
    # at the above, but to be explicit, let's save
    # it here.
    size::CellRect
    # For mouse tracking
    subscreens::Vector{TrackingSubscreen}
    topwidget::Nullable{Widget}
    # For iterm2 image support
    images::Vector{Tuple{Vector{UInt8},CellRect}}
    freelist::IntSet     # Which indicies to recycle
    # Which indicies were dropped by the widget and can
    # be recycled if they are no longer in the image
    dropped::IntSet
    # For embedding support
    offset::CellLoc
    fullsize::CellRect
end
istracking(s::DoubleBufferedTerminalScreen) = true
topwidget(s::DoubleBufferedTerminalScreen) = get(s.topwidget)
subscreens(s::DoubleBufferedTerminalScreen) = s.subscreens
show(io::IO,s::DoubleBufferedTerminalScreen) =
    print(io,"DoubleBufferedTerminalScreen(",s.size[2],"x",s.size[1],")")

function resized(s::DoubleBufferedTerminalScreen,size)
    s.have = Vector{Line}()
    s.want = Vector{Line}()
    s.wantcursor = Cursor(1,1)
    if isfullscreen(s)
        s.size = size
    end
    s.fullsize = size
    empty!(s.subscreens)
    s.topwidget = Nullable{Widget}()
    s.want = emptybuffer(s)
    s
end

function clear(s::DoubleBufferedTerminalScreen)
    s.have = Vector{Line}()
    s.want = emptybuffer(s)
end

function invalidate(s::DoubleBufferedTerminalScreen)
    s.have = Vector{Line}()
    push!(invalidated2,s)
end

function allocate_char_for_image(s::DoubleBufferedTerminalScreen, data::Vector{UInt8}, size)
    if !isempty(s.freelist)
        idx = pop!(s.freelist)
        s.images[idx] = (data, size)
    else
        push!(s.images, (data, size))
        idx = length(s.images)
    end
    return Cell(Cell(Char(0xF0000+idx)), flags = CELL_IS_IMG)
end

function subscreen_for_pos(s::Screen,subscreens::Vector{TrackingSubscreen}, pos)
    for ss in subscreens
        if pos[1] in simple(ss).rows && pos[2] in simple(ss).cols
            return subscreen_for_pos(ss, ss.children, translate⁻¹(ss,pos))
        end
    end
    return (s, pos)
end

function add_child(s::DoubleBufferedTerminalScreen, ss::TrackingSubscreen)
    push!(s.subscreens, ss)
end

function redraw(s::DoubleBufferedTerminalScreen, w::Widget)
    empty!(s.subscreens)
    s.topwidget = w
    draw(s,w)
end

width(s::DoubleBufferedTerminalScreen) = s.size[2]
height(s::DoubleBufferedTerminalScreen) = s.size[1]
size(s::DoubleBufferedTerminalScreen) = s.size

isfullscreen(s::DoubleBufferedTerminalScreen) = s.offset == (1,1) && s.size == s.fullsize

function DoubleBufferedTerminalScreen(size; offset = (1,1), fullsize = size)
    s = DoubleBufferedTerminalScreen(
        Vector{Line}(),
        Vector{Line}(),
        Cursor(1,1),
        size,
        Vector{TrackingSubscreen}(),
        Nullable{Widget}(),
        Vector{Tuple{Vector{UInt8},CellRect}}(),
        IntSet(), IntSet(),
        offset,fullsize)
    s.want = emptybuffer(s)
    s
end

function invalidbuffer(s::DoubleBufferedTerminalScreen)
    Line[Line(Cell[Cell(0) for _ in 1:s.size[2]]) for _ in 1:s.size[1]]
end

function emptybuffer(s::DoubleBufferedTerminalScreen)
    Line[Line(Cell[Cell(' ') for _ in 1:s.size[2]]) for _ in 1:s.size[1]]
end

for f in (:cmove, :cmove_col, :cmove_line_up, :cmove_line_down,
    :cmove_up, :cmove_down)
    @eval function $(f)(s::DoubleBufferedTerminalScreen,args...)
        s.wantcursor = $(f)(s.wantcursor,args...)
    end
end

fresh_screen(s::DoubleBufferedTerminalScreen) = empty!(s.have)

function resize!(s::DoubleBufferedTerminalScreen,size)
    error("DoubleBufferedTerminalScreen cannot be resized. Create a new screen instead")
end

# Will be filled in in swrite write previous color
ascellparams(color::Void) = (0,RGB8(0,0,0),false)
ascellparams(symbol::Symbol) = (VT100.colorlist[symbol],RGB8(0,0,0),false)
ascellparams(c::Color)  = (0,convert(RGB8,c),true)

function dummycell(fg,bg,attrs = 0)
    fg, fg_rgb, fg_is_rgb = ascellparams(fg)
    bg, bg_rgb, bg_is_rgb = ascellparams(bg)
    Cell(Char(0),0 | (fg_is_rgb ? FG_IS_RGB : 0) |
        (bg_is_rgb ? BG_IS_RGB : 0),fg,bg,attrs,fg_rgb,bg_rgb)
end

function update_cell!(line,idx,dcell,char,fg,bg)
    oldcell = line[idx]
    cell = Cell(dcell; content = char)
    fg === nothing && (cell = Cell(cell; fg = oldcell.fg, fg_rgb = oldcell.fg_rgb))
    bg === nothing && (cell = Cell(cell; bg = oldcell.bg, bg_rgb = oldcell.bg_rgb))
    line[idx] = cell
end

function swrite(s::DoubleBufferedTerminalScreen, rows, cols, cell::Cell)
    for r in rows
        line = s.want[r]
        for c in cols
            line[c] = cell
        end
    end
end

function swrite(s::DoubleBufferedTerminalScreen, rows, cols, char::Char; fg = :default, bg = :default, attrs = 0)
    if attrs & IsACS != 0
        i = findfirst(VT100.LineDrawing,char)
        @assert i != -1
        char = Char(i-1)
    end
    dcell = dummycell(fg,bg,attrs)
    for r in rows
        line = s.want[r]
        for c in cols
            update_cell!(line,c,dcell,char,fg,bg)
        end
    end
end

function swrite_image(s::DoubleBufferedTerminalScreen, rows, cols, data::Vector{UInt8})
    c = allocate_char_for_image(s,data,CellRect(length(rows),length(cols)))
    swrite_image(s,rows,cols,c)
    c
end

function swrite_image(s::DoubleBufferedTerminalScreen, rows, cols, c::Cell)
    @assert (c.flags & CELL_IS_IMG) != 0
    for row in rows
        line = s.want[row]
        for col in cols
            line[col] = get_image_cell(c,row-first(rows)+1,col-first(cols)+1)
        end
    end
    c
end

function clear(s::DoubleBufferedTerminalScreen, rows, cols)
    swrite(s, rows, cols, ' ')
end

function swrite(s::DoubleBufferedTerminalScreen, rows, cols, string::AbstractString;
        fg = :default, bg = :default, attrs = 0)
    @assert attrs & IsACS == 0
    dcell = dummycell(fg,bg,attrs)
    for r in rows
        line = s.want[r]
        for c in cols
            for (i,char) in enumerate(string)
                update_cell!(line,c+i-1,dcell,char,fg,bg)
            end
        end
    end
end

function do_erase(buf,n)
    if n <= 5
        for i = 1:n
            write(buf,' ')
        end
    else
        write(buf,CSI,string(n),'X') # ECH
        do_move(buf,n)
    end
end

function do_clear_rect(buf,box::CellBox)
    tl = box.topleft
    top = tl.row
    left = tl.col
    bot = top + box.size.height
    right = left + box.size.width

    # For terminals that support DECERA
    # write(buf,CSI,top,';',left,';',bot,';',right,"\$z")
    for i in top:bot
        do_absmove(buf,i,left)
        write(buf,CSI,string(right-left),'X') # ECH
    end
end

function do_absmove(buf,row,col)
    write(buf, CSI, string(row), ';', string(col), 'H') # CUP
end

do_absmove(buf,s,row,col) = do_absmove(buf,s.offset[1]+row-1,s.offset[2]+col-1)

function do_move(buf,n)
    write(buf, CSI, string(n), 'C') # CUF
end

function change_color(buf, flags, color, rgb, is_256, is_rgb)
    if (flags & is_256) != 0
        write(buf,"8;5;",string(color))
    elseif (flags & is_rgb) != 0
        write(buf,"8;2;",string(rgb.r.i),';',
                          string(rgb.g.i),';',
                          string(rgb.b.i))
    else
        write(buf,string(color))
    end
    write(buf,'m')
end

function change_fg_color(buf, cell)
    write(buf,CSI,'3')
    change_color(buf, cell.flags, cell.fg, cell.fg_rgb, FG_IS_256, FG_IS_RGB)
end

function change_bg_color(buf, cell)
    write(buf,CSI,'4')
    change_color(buf, cell.flags, cell.bg, cell.bg_rgb, BG_IS_256, BG_IS_RGB)
end

function change_attrs(buf, want, have)
    # If we are clearing any bits
    if (want & IsACS) != (have & IsACS)
        write(buf,(want & IsACS) == 0 ? "\e(B" : "\e(0")
    end
    if (have & ~want) != 0
        write(buf,CSI,"0m")
        have = 0
    end
    if (want & Bright) != 0
        write(buf,CSI,"1m")
    end
end

function render_cell(buf,wantc,attrs)
    (wantc.bg != attrs.bg || wantc.bg_rgb != attrs.bg_rgb) && change_bg_color(buf,wantc)
    (wantc.fg != attrs.fg || wantc.fg_rgb != attrs.fg_rgb) && change_fg_color(buf,wantc)
    (wantc.attrs != attrs.attrs) && change_attrs(buf,wantc.attrs,attrs.attrs)
    write(buf,wantc.content)
    wantc
end

same_attributes(havec,wantc) = Cell(havec; content = '\0') == Cell(wantc; content = '\0')

include("iterm2.jl")

"""
Core screen rendering functionality. Loop through the screen and update all
the cells that have changed since the last time we wenth through (in the most
efficient way possible).

- We do images first because they need to be written out as a chunk, and it
  is not worth the incresed complexity to go back and figure out exactly
  which areas we'd have to reprocess.

"""
function render(s::DoubleBufferedTerminalScreen)
    buf = IOBuffer()
    # For now use a cell to hold the attributestate
    attributestate = dummycell(:default,:default)
    full = isfullscreen(s)
    if !full
        # If our subscreen extends off the screen we may have to scroll first
        if s.offset.row + s.size.height > s.fullsize.height
            difference = s.offset.row + s.size.height-s.fullsize.height
            write(buf,CSI,string(difference),'S')
            s.offset = CellLoc(s.offset.row-difference,s.offset.col)
        end
    end
    change_fg_color(buf,attributestate)
    change_bg_color(buf,attributestate)
    # Make sure everything gets reset
    change_attrs(buf,0,0xFF)
    if isempty(s.have)
        s.have = emptybuffer(s)
        if full
            write(buf,CSI,'H',CSI,"2J")
        else
            do_clear_rect(buf,CellBox(s.offset,s.size))
        end
    end
    # Stage 1: Render images
    images_to_draw = Dict{Tuple{Int,Int},Int}()
    for (row,(havel,wantl)) in enumerate(zip(s.have,s.want))
        for (j,(havec,wantc)) in enumerate(zip(havel,wantl))
            if (wantc.flags & CELL_IS_IMG) == 0 ||
                havec == wantc
                continue
            end
            cur_pos = pos_for_image_cell(wantc)
            start_pos = (row-cur_pos[1]+1,j-cur_pos[2]+1)
            if !haskey(images_to_draw,start_pos)
                images_to_draw[start_pos] = wantc.content - 0xF0000
            end
        end
    end

    todelete = copy(s.dropped)
    for (pos,image) in images_to_draw
        # Drop this image from the todelete list if it was on there
        pop!(todelete, image, image)
        data, size = s.images[image]
        do_absmove(buf,s,pos...)
        display_file(buf, data, height = size.height, width = size.width, inline = true)
        #@show (pos,size)
        for row in pos[1]:(pos[1]+size.height-1)
            line = s.have[row]
            for col in pos[2]:(pos[2]+size.width-1)
                line[col] = get_image_cell(Cell(Char(image+0xF0000)),row-pos[1]+1,col-pos[2]+1)
            end
        end
    end
    for image in todelete
        s.images[image] = (Vector{UInt8}(),CellRect(0,0))
        pop!(s.dropped,image)
        push!(s.freelist,image)
    end

    # Stage 2: Every thing else
    for (row,(havel,wantl)) in enumerate(zip(s.have,s.want))
        didmoverow = false
        nerase = 0
        nmove = 0
        # Look through the cells
        for (j,(havec,wantc)) in enumerate(zip(havel,wantl))
            # Update the have buffer in place
            havel[j] = wantc
            # Now do the rendering
            if havec == wantc
                if nerase != 0
                    do_erase(buf, nerase)
                    nerase = 0
                end
                nmove += 1
                continue
            end
            # First we need to move to the right spot
            if nmove != 0
                if !didmoverow
                    do_absmove(buf,s,row,nmove+1)
                    didmoverow = true
                else
                    do_move(buf, nmove)
                end
                nmove = 0
            end
            if wantc.content == ' ' && same_attributes(wantc,attributestate)
                # This can be handled by an erase
                nerase += 1
                continue
            else
                @assert (wantc.flags & CELL_IS_IMG) == 0
                if !didmoverow
                    do_absmove(buf,s,row,1)
                    didmoverow = true
                end
                if nerase != 0
                    do_erase(buf, nerase)
                    nerase = 0
                end
                attributestate = render_cell(buf,wantc,attributestate)
            end
        end
    end
    change_attrs(buf,0,0xFF)
    buf
end

function afterembed(io::IO,s::DoubleBufferedTerminalScreen)
    # Go to the line after the embedding
    # However, if we are at the end, we might have to scroll a line
    if s.offset.row + s.size.height == s.fullsize.height
        write(io,CSI,"1S")
    end
    write(io,CSI,"$(s.offset.row + s.size.height);1H")
end

function render(io::IO,s::DoubleBufferedTerminalScreen)
    write(io,takebuf_array(render(s)))
end


text:
#!/usr/bin/env julia

using Pkg
haskey(Pkg.installed(), "Documenter") || Pkg.add("Documenter")

using Documenter
using NarrativeTest

# Highlight indented code blocks as Julia code.
using Markdown
Markdown.Code(code) = Markdown.Code("julia", code)

makedocs(
    sitename = "NarrativeTest.jl",
    pages = [
        "Home" => "index.md",
        "guide.md",
        "reference.md",
        "test.md",
    ],
    modules = [NarrativeTest])

deploydocs(
    repo = "github.com/rbt-lang/NarrativeTest.jl.git",
)