text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 95eaf04a-08e5-11eb-043f-27646f1bf207
using DelimitedFiles, Dictionaries, IdentityRanges, LsqFit, NPZ, OffsetArrays, Plots, Plots.Measures, SparseArrays, Statistics

# ╔═╡ bfbc2e08-0c3b-11eb-10c0-7bd1d3a723b9
# Code to load CSC matrix from NPZ file
# based on description here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html
# see also:
# + https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.save_npz.html
# + https://docs.julialang.org/en/v1/stdlib/SparseArrays/#man-csc

module NPZ_CSC
using PyCall: pyimport, pyrepr
using SparseArrays: SparseMatrixCSC

np = pyimport("numpy")

function load(filename)
	contents = np.load(filename, allow_pickle=true)
	format = get(contents, "format")
	pyrepr(format) == "array(b'csc', dtype='|S3')" ||
		throw("Array format seems wrong. Got $(pyrepr(format))")

	shape = get(contents, "shape")
	data = get(contents, "data")
	indptr = get(contents, "indptr") .+ 1
	indices = get(contents, "indices") .+ 1

	SparseMatrixCSC(shape[1], shape[2], indptr, indices, data)
end

end

# ╔═╡ 849adbca-08e5-11eb-1956-cd02f2a6be7a
md"""
# Plots for two-wave epidemic
"""

# ╔═╡ 34f9cdda-0901-11eb-1117-099323e644ef
gr(linewidth=2,markerstrokecolor=:auto,markersize=1.5,legend=nothing,
	dpi=200,size=(900,400),titlefontsize=9,guidefontsize=9,tickfontsize=7)

# ╔═╡ 289fbed4-21ff-11eb-261a-25eace4cbf4e
# polynomial fit
function polyfit(xx,yy,order)
	axes(xx) == axes(yy) || throw("Arrays must have same axes. Got: $(axes(xx)) and $(axes(yy))")
	x, y = OffsetArrays.no_offset_view(xx), OffsetArrays.no_offset_view(yy)

	fit_func = (x,p) -> evalpoly.(x,Ref(p))
	popt = curve_fit(fit_func,x,y,ones(order)).param
	return fit_func(xx,popt)
end

# ╔═╡ 6539e014-543f-11eb-17c3-991de18e923c
# plot a rectangle
rectangle(xlims,ylims) = Shape([xlims; reverse(xlims)],repeat(ylims,inner=2))

# ╔═╡ eb7bec9e-5453-11eb-0876-c93037e57988
hypername((n,m,q)) = "HYPER (m=$m, q=$q)"

# ╔═╡ 210d331a-09c3-11eb-06ed-e1f4ddc023fc
md"""
## Load data
"""

# ╔═╡ fd643c2c-08e5-11eb-0d25-014798a505d1
md"""
### Utility functions
"""

# ╔═╡ fe7cf842-08e5-11eb-0d07-41ed2902921a
function loadresults(dir, pattern)
	regmatches = filter(!isnothing, match.(pattern, readdir(dir)))
	pairs = map(regmatches) do regmatch
		design = (n=parse(Int,regmatch[:n]), m=parse(Int,regmatch[:m]), q=parse(Int,regmatch[:q]))
		result = OffsetVector(npzread(joinpath(dir,regmatch.match)),-1)
		design => result
	end
	return dictionary(pairs)
end

# ╔═╡ 3f1b753e-0901-11eb-0060-1f1269588bef
loadeff(dir; pattern=r"^Eff_avg.n-(?<n>\d+)_m-(?<m>\d+)_q-(?<q>\d+).npy$") = loadresults(dir, pattern)

# ╔═╡ 3f1c5288-0901-11eb-391f-dd7dd02d3d9b
loadsens(dir; pattern=r"^Recall_combined.n-(?<n>\d+)_m-(?<m>\d+)_q-(?<q>\d+).npy$") = loadresults(dir, pattern)

# ╔═╡ 05c6923e-08e6-11eb-1d9b-21afc86a6b04
md"""
### Viral load distributions and average prevalence
"""

# ╔═╡ c6c3f3ca-0c3b-11eb-2593-439620bd9e78
ViralLoad = NPZ_CSC.load("Simulated_populations_two_wave/seir_viral_loads_swab_switch_SEIR_new_1.viral_loads.npz")

# ╔═╡ cf1b2b1a-0c3b-11eb-2159-6ffeb682d595
prevalence = OffsetArray(vec(mean(>(0), ViralLoad, dims=1)), -1)

# ╔═╡ 122151e0-08e6-11eb-18f5-3b277145bc9d
md"""
### HYPER
"""

# ╔═╡ 47d9acce-2fa6-11eb-0260-7767c0252b8c
hypereff = loadeff("designs/hypergraph/")

# ╔═╡ f8274a32-2fa6-11eb-3836-9bbe0aaf63cf
hypersens = loadsens("designs/hypergraph/")

# ╔═╡ 2e0cc554-08e6-11eb-0186-bba8c2aee7ec
md"""
### Array
"""

# ╔═╡ 339ca8ce-08e6-11eb-19c3-5f0608ee2794
arrayeff = (dictionary∘Dict)(
	(8,12) => only(loadeff("designs/array,8x12/")),
	(16,24) => only(loadeff("designs/array,16x24/"))
)

# ╔═╡ 33a20596-0902-11eb-2d50-cb33f008ffd9
arraysens = (dictionary∘Dict)(
	(8,12) => only(loadsens("designs/array,8x12/")),
	(16,24) => only(loadsens("designs/array,16x24/"))
)

# ╔═╡ 3706cb8e-08e6-11eb-3c27-f1142c9d9159
md"""
### P-BEST
"""

# ╔═╡ 6f20ef1a-2216-11eb-37d8-cb470acc99e4
pbesteff = only(loadeff("designs/pbest/"))

# ╔═╡ 744ec278-2216-11eb-2aee-6138f60f7722
pbestsens = only(loadsens("designs/pbest/"))

# ╔═╡ a92db3c4-5436-11eb-181b-3f1a71204b12
md"""
### Individual testing
"""

# ╔═╡ 936c39ce-368e-11eb-0cb2-938376a6c18f
LOD = 100.0

# ╔═╡ 72f1ae0a-544c-11eb-0a55-f18f192d64c2
indiveff = [1.0 for _ in eachindex(prevalence)]

# ╔═╡ 21073ca8-3684-11eb-0a25-8d16a9799cf1
indivsens = mapslices(ViralLoad, dims=1) do vl
	posloads = filter(>(0),vl)
	isempty(posloads) ? 1.0 : mean(>(LOD),posloads)
end |> x->OffsetArray(vec(x),-1)

# ╔═╡ 0fd98c50-2203-11eb-0f58-d946dcf484b4
md"""
## Figures
"""

# ╔═╡ 54b274ec-a780-11eb-3a85-7bc7a4f3d0cd
md"""
### Figs a,b
"""

# ╔═╡ 24d773d0-b9d2-11eb-009d-59905303c23f
histprev = 0.01

# ╔═╡ 5afac4cc-b9d1-11eb-0af1-3d9bb7a112dd
histdays = findall(i->prevalence[i] < histprev < prevalence[i+1] || prevalence[i] > histprev > prevalence[i+1],0:199)

# ╔═╡ 0045d946-a2ce-11eb-33f3-a30397432557
let
	prevplt = plot(prevalence;
		xlabel="day",xlims=(10,190),xticks=10:20:190,
		ylabel="prevalence",ylims=(0,0.02),
		yticks=(0:0.005:0.02,string.(0:0.5:2,"%")),
		color=palette(:default)[2])

	bins=-4.5:0.1:21
	distplt = plot(palette=PlotThemes.wong_palette)
	for day in histdays
		vl = filter(!iszero,ViralLoad[:,day]) |> collect .|> log10
		histogram!(vl;bins,opacity=0.5,linewidth=0,normalize=:pdf,
			label="Day $day: prevalence of $(round(100*prevalence[day],digits=1))%")
	end
	plot!(legend=:topright,bglegend=nothing,fglegend=nothing,legendfontsize=6,
		xlabel="log10(viral load)",yaxis=false,yticks=nothing,
		xlims=(-5,21),xticks=-4:2:20,
		left_margin=12.75mm)

	annotate!(prevplt,-11,0.0205,text("a",14))
	annotate!(distplt,-5,0.18,text("b",14))

	plot(prevplt,distplt;size=(900,150),bottom_margin=4mm)
	savefig("fig-a,b.png"); plot!()
end

# ╔═╡ 71950c84-543b-11eb-3c6b-4b4114ade7bf
md"""
### Fig c
"""

# ╔═╡ 6c66ebe0-543b-11eb-1619-294095a5ae7a
let compfigs = [96=>(n=96,m=16,q=2),384=>(n=384,m=32,q=2)]
	# Plot colors
	dcolors = distinguishable_colors(7, [RGB(1,1,1)], dropseed=true)
	indivcolor = dcolors[1]
	pbestcolor = dcolors[2]
	arraycolor = dcolors[3]
	hypercolor = Dict(
		16=>dcolors[4],
		32=>dcolors[5],
	)

	# Efficiency gain plots
	effplts = map(compfigs) do (n,hyperdesign)
		plot()

		# Array
		plot!(arrayeff[n == 96 ? (8,12) : (16,24)],linestyle=:solid,color=arraycolor)

		# P-BEST
		if n == 384
			plot!(pbesteff,linestyle=:solid,color=pbestcolor)
		end

		# HYPER
		plot!(hypereff[hyperdesign],linestyle=:dash,color=hypercolor[hyperdesign.m])

		# Individual
		plot!(indiveff,linestyle=:dot,color=indivcolor)
	end

	# Sensitivity plots
	sensplts = map(compfigs) do (n,hyperdesign)
		plot()

		# Array
		scatter!(arraysens[n == 96 ? (8,12) : (16,24)],opacity=0.6,color=arraycolor)

		# P-BEST
		if n == 384
			scatter!(pbestsens,opacity=0.6,color=pbestcolor)
		end

		# HYPER
		scatter!(hypersens[hyperdesign],opacity=0.6,color=hypercolor[hyperdesign.m])

		# Individual
		scatter!(indivsens,opacity=0.6,color=indivcolor)
	end

	# Figure (titles, axes, annotations)
	title!.(effplts,string.("n=",first.(compfigs)," individuals"))

	plot!.(sensplts,xlims=(10,190),xticks=(10:20:190,["$day\n($(round(100*prevalence[day],digits=2))%)" for day in 10:20:190]),xlabel="day (prevalence)")
	plot!.(effplts, xlims=(10,190),xticks=(10:20:190,""))

	plot!(effplts[1],ylims=(0, 6*1.14),yticks=[0,1,96/20,96/16])
	plot!(effplts[2],ylims=(0,12*1.14),yticks=[0,1,384/40,384/48,384/32])
	plot!.(sensplts,ylims=(0.43,0.98),yticks=(0.5:0.1:0.9,string.(50:10:90,"%")))
	plot!(effplts[1], ylabel="individuals/test")
	plot!(sensplts[1],ylabel="sensitivity")

	annotate!(effplts[1],[
		(11.75, 6.1, text(hypername(last(compfigs[1])), 8, :left, :bottom, hypercolor[last(compfigs[1]).m])),
		(11.75, 4.9, text("8x12 array",                 8, :left, :bottom, arraycolor)),
		(11.75, 1.1, text("Individual testing",         8, :left, :bottom, indivcolor))
	])

	annotate!(effplts[2],[
		(11.75, 12.2, text(hypername(last(compfigs[2])), 8, :left, :bottom, hypercolor[last(compfigs[2]).m])),
		(11.75,  9.8, text("16x24 array",                8, :left, :bottom, arraycolor)),
		(11.75,  7.8, text("P-BEST",                     8, :left, :top,    pbestcolor)),
		(11.75,  1.2, text("Individual testing",         8, :left, :bottom, indivcolor))
	])

	annotate!(effplts[1],-11,7.25,text("c",14))

	plot!.(effplts, bottom_margin=0mm)
	plot!.(sensplts,bottom_margin=5mm,top_margin=-2.5mm)
	plot!.([effplts; sensplts],right_margin=3mm)

	plot(effplts...,sensplts...,layout=(2,2),link=:x)
	savefig("fig-c.png"); plot!()
end

# ╔═╡ 5b53ff88-5458-11eb-1313-b74a1235ad08
md"""
### Figs d-e
"""

# ╔═╡ 7eea53ac-5457-11eb-39e2-19c11bcff57d
let compfigs = [384=>[(n=384,m=32,q=2),(n=384,m=16,q=2),(n=384,m=12,q=2)],384=>[(n=384,m=12,q=1),(n=384,m=12,q=2),(n=384,m=12,q=3)]]
	# Plot colors and linestyles
	dcolors = distinguishable_colors(7, [RGB(1,1,1)], dropseed=true)
	indivcolor = dcolors[1]
	hypercolor = Dict(
		(16,2)=>dcolors[4],
		(32,2)=>dcolors[5],
		(12,2)=>dcolors[7],
		(12,3)=>PlotThemes.wong_palette[2],
		(12,1)=>PlotThemes.wong_palette[1],
	)
	hyperstyle = Dict(1=>:dot, 2=>:dash, 3=>:solid)

	# Efficiency gain plots
	effplts = map(compfigs) do (n,hyperdesigns)
		plot()

		# HYPER
		for design in hyperdesigns
			plot!(hypereff[design],linestyle=hyperstyle[design.q],color=hypercolor[(design.m,design.q)])
		end

		# Individual
		plot!(indiveff,linestyle=:dot,color=indivcolor)
	end

	# Sensitivity plots
	sensplts = map(compfigs) do (n,hyperdesigns)
		plot()

		# HYPER
		for design in hyperdesigns
			scatter!(hypersens[design],opacity=0.6,color=hypercolor[(design.m,design.q)])
		end

		# Individual
		scatter!(indivsens,opacity=0.6,color=indivcolor)
	end

	# Figure (titles, axes, annotations)
	title!.(effplts,string.("n=",first.(compfigs)," individuals"))

	plot!.(sensplts,xlims=(10,190),xticks=(10:20:190,["$day\n($(round(100*prevalence[day],digits=2))%)" for day in 10:20:190]),xlabel="day (prevalence)")
	plot!.(effplts, xlims=(10,190),xticks=(10:20:190,""))

	plot!.(effplts, ylims=(0,32*1.02),yticks=[1,384/32,384/16,384/12])
	plot!.(sensplts,ylims=(0.37,0.98),yticks=(0.4:0.1:0.9,string.(40:10:90,"%")))
	plot!(effplts[1], ylabel="individuals/test")
	plot!(sensplts[1],ylabel="sensitivity")

	annotate!(effplts[1],[
		(11.75, 11.2, text(hypername(last(compfigs[1])[1]), 8, :left, :top, hypercolor[(last(compfigs[1])[1].m,last(compfigs[1])[1].q)])),
		(11.75, 20.0, text(hypername(last(compfigs[1])[2]), 8, :left, :top, hypercolor[(last(compfigs[1])[2].m,last(compfigs[1])[2].q)])),
		(50.75, 32.0, text(hypername(last(compfigs[1])[3]), 8, :left, :top, hypercolor[(last(compfigs[1])[3].m,last(compfigs[1])[3].q)])),
		(11.75,  1.4, text("Individual testing",            8, :left, :bottom, indivcolor))
	])
	annotate!(effplts[1],-11,35,text("d",14))

	annotate!(effplts[2],[
		(11.75, 12.2, text(hypername(last(compfigs[2])[1]), 8, :left, :top, hypercolor[(last(compfigs[2])[1].m,last(compfigs[2])[1].q)])),
		(12.75, 31.4, text(hypername(last(compfigs[2])[2]), 8, :left, :top, hypercolor[(last(compfigs[2])[2].m,last(compfigs[2])[2].q)], rotation=-20)),
		(55.75, 32.0, text(hypername(last(compfigs[2])[3]), 8, :left, :top, hypercolor[(last(compfigs[2])[3].m,last(compfigs[2])[3].q)])),
		(11.75,  1.4, text("Individual testing",            8, :left, :bottom, indivcolor))
	])
	annotate!(effplts[2],-7,35,text("e",14))

	plot!.(effplts, bottom_margin=0mm)
	plot!.(sensplts,bottom_margin=5mm,top_margin=-2.5mm)
	plot!.([effplts; sensplts],right_margin=3mm)

	plot(effplts...,sensplts...,layout=(2,2),link=:x)
	savefig("fig-d,e.png"); plot!()
end

# ╔═╡ Cell order:
# ╟─849adbca-08e5-11eb-1956-cd02f2a6be7a
# ╠═95eaf04a-08e5-11eb-043f-27646f1bf207
# ╠═34f9cdda-0901-11eb-1117-099323e644ef
# ╟─289fbed4-21ff-11eb-261a-25eace4cbf4e
# ╟─6539e014-543f-11eb-17c3-991de18e923c
# ╟─eb7bec9e-5453-11eb-0876-c93037e57988
# ╟─210d331a-09c3-11eb-06ed-e1f4ddc023fc
# ╟─fd643c2c-08e5-11eb-0d25-014798a505d1
# ╟─fe7cf842-08e5-11eb-0d07-41ed2902921a
# ╟─3f1b753e-0901-11eb-0060-1f1269588bef
# ╟─3f1c5288-0901-11eb-391f-dd7dd02d3d9b
# ╟─05c6923e-08e6-11eb-1d9b-21afc86a6b04
# ╟─bfbc2e08-0c3b-11eb-10c0-7bd1d3a723b9
# ╟─c6c3f3ca-0c3b-11eb-2593-439620bd9e78
# ╟─cf1b2b1a-0c3b-11eb-2159-6ffeb682d595
# ╟─122151e0-08e6-11eb-18f5-3b277145bc9d
# ╟─47d9acce-2fa6-11eb-0260-7767c0252b8c
# ╟─f8274a32-2fa6-11eb-3836-9bbe0aaf63cf
# ╟─2e0cc554-08e6-11eb-0186-bba8c2aee7ec
# ╟─339ca8ce-08e6-11eb-19c3-5f0608ee2794
# ╟─33a20596-0902-11eb-2d50-cb33f008ffd9
# ╟─3706cb8e-08e6-11eb-3c27-f1142c9d9159
# ╟─6f20ef1a-2216-11eb-37d8-cb470acc99e4
# ╟─744ec278-2216-11eb-2aee-6138f60f7722
# ╟─a92db3c4-5436-11eb-181b-3f1a71204b12
# ╟─936c39ce-368e-11eb-0cb2-938376a6c18f
# ╟─72f1ae0a-544c-11eb-0a55-f18f192d64c2
# ╟─21073ca8-3684-11eb-0a25-8d16a9799cf1
# ╟─0fd98c50-2203-11eb-0f58-d946dcf484b4
# ╟─54b274ec-a780-11eb-3a85-7bc7a4f3d0cd
# ╠═24d773d0-b9d2-11eb-009d-59905303c23f
# ╠═5afac4cc-b9d1-11eb-0af1-3d9bb7a112dd
# ╟─0045d946-a2ce-11eb-33f3-a30397432557
# ╟─71950c84-543b-11eb-3c6b-4b4114ade7bf
# ╟─6c66ebe0-543b-11eb-1619-294095a5ae7a
# ╟─5b53ff88-5458-11eb-1313-b74a1235ad08
# ╟─7eea53ac-5457-11eb-39e2-19c11bcff57d


text:
##################################################
# Check if the obtained transfer matrices are
# markov.
##################################################

"""
    is_markov(M::AbstractArray{T, 2}) -> Bool

Is the matrix Markov? """
function is_markov(M::AbstractArray{T, 2}) where T
    all(sum(M, dims = 2) .≈ 1)
end

"""
    is_markov(TO::AbstractTransferOperator) -> Bool

Is the transfer operator Markov? """
function is_markov(TO::AbstractTransferOperator)
    is_markov(TO.transfermatrix)
end

"""
    is_almostmarkov(M::AbstractArray{T, 2}; tol = 0.01) -> Bool

Is the matrix almost Markov?
"""
function is_almost_markov(M::AbstractArray{T, 2}; tol = 1e-3) where T
    all(sum(M, dims = 2) .> (1 - tol))
end

"""
    is_almostmarkov(TO::AbstractTransferOperator; tol = 0.01) -> Bool

Is the transfer operator almost Markov?
"""
function is_almost_markov(TO::AbstractTransferOperator; tol = 1e-3)
    is_almost_markov(TO.transfermatrix)
end

function zerocols(M::AbstractArray{T, 2}) where T
    findall(sum(M, dims = 1) .== 0)
end

function zerorows(M::AbstractArray{T, 2}) where T
    findall(sum(M, dims = 2) .== 0)
end

function zerocols(TO::AbstractTransferOperator)
    findall(sum(TO.transfermatrix, dims = 1) .== 0)
end

function zerorows(TO::AbstractTransferOperator)
    findall(sum(TO.transfermatrix, dims = 2) .== 0)
end

export
zerocols,
zerorows,
is_markov,
is_almost_markov

# #####################################################
# # How much do the exact and approximate triangulation
# # estimators differ?
# #####################################################
# """
#     max_discrep(to1::ExactSimplexTransferOperator,
#                 to2::ApproxSimplexTransferOperator)
#
# Measures the maximum discrepancy between an exact and an approximate estimate
# of the transfer operator.
#
# For this function to work, the transfer operators must have been generated from
# the same triangulation (otherwise, the number of simplices don't match up).
# """
# function max_discrep(exact::ExactSimplexTransferOperator,
#                      approx::ApproxSimplexTransferOperator)
#
# end

ext: jl


text:
"""
"""
abstract type AssemblyStrategy end

"""
"""
function row_map(a::AssemblyStrategy,row)
  @abstractmethod
end

"""
"""
function col_map(a::AssemblyStrategy,col)
  @abstractmethod
end

"""
"""
function row_mask(a::AssemblyStrategy,row)
  @abstractmethod
end

"""
"""
function col_mask(a::AssemblyStrategy,col)
  @abstractmethod
end

struct DefaultAssemblyStrategy <: AssemblyStrategy end

row_map(a::DefaultAssemblyStrategy,row) = row

col_map(a::DefaultAssemblyStrategy,col) = col

row_mask(a::DefaultAssemblyStrategy,row) = true

col_mask(a::DefaultAssemblyStrategy,col) = true

"""
"""
abstract type Assembler <: GridapType end

"""
"""
function get_test(a::Assembler)
  @abstractmethod
end

"""
"""
function get_trial(a::Assembler)
  @abstractmethod
end

"""
"""
function get_assembly_strategy(a::Assembler)
  @abstractmethod
end

"""
"""
function allocate_matrix(a::Assembler,matdata)
  @abstractmethod
end

"""
"""
function allocate_vector(a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function allocate_matrix_and_vector(a::Assembler,data)
  @abstractmethod
end

"""
"""
function assemble_matrix!(A,a::Assembler,matdata)
  @abstractmethod
end

"""
"""
function assemble_matrix_add!(mat,a::Assembler, matdata)
  @abstractmethod
end

"""
"""
function assemble_vector!(b,a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function assemble_vector_add!(b,a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function assemble_matrix_and_vector!(A,b,a::Assembler, data)
  @abstractmethod
end

function assemble_matrix_and_vector_add!(A,b,a::Assembler, data)
  @abstractmethod
end

"""
"""
function assemble_matrix(a::Assembler,matdata)
  A = allocate_matrix(a,matdata)
  assemble_matrix!(A,a,matdata)
  A
end

"""
"""
function assemble_vector(a::Assembler,vecdata)
  b = allocate_vector(a,vecdata)
  assemble_vector!(b,a,vecdata)
  b
end

"""
"""
function assemble_matrix_and_vector(a::Assembler,data)
  A, b = allocate_matrix_and_vector(a,data)
  assemble_matrix_and_vector!(A,b,a,data)
  (A, b)
end

"""
"""
function test_assembler(a::Assembler,matdata,vecdata,data)
  trial_fesp = get_trial(a)
  test_fesp = get_test(a)
  A = allocate_matrix(a,matdata)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  assemble_matrix!(A,a,matdata)
  assemble_matrix_add!(A,a,matdata)
  A = assemble_matrix(a,matdata)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  b = allocate_vector(a,vecdata)
  @test num_free_dofs(test_fesp) == length(b)
  assemble_vector!(b,a,vecdata)
  assemble_vector_add!(b,a,vecdata)
  b = assemble_vector(a,vecdata)
  @test num_free_dofs(test_fesp) == length(b)
  A, b = allocate_matrix_and_vector(a,data)
  assemble_matrix_and_vector!(A,b,a,data)
  assemble_matrix_and_vector_add!(A,b,a,data)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  @test num_free_dofs(test_fesp) == length(b)
  A, b = assemble_matrix_and_vector(a,data)
  @test num_free_dofs(trial_fesp) == size(A,2)
  @test num_free_dofs(test_fesp) == size(A,1)
  @test num_free_dofs(test_fesp) == length(b)
end


text:
Lar = LinearAlgebraicRepresentation

using IntervalTrees
using SparseArrays
using NearestNeighbors
using DataStructures

#-------------------------------------------------------------------------------
#   MODEL BOUNDING BOXES
#-------------------------------------------------------------------------------

"""
	getModelBoundingBoxes(
		model::CAGD.Model,
		deg::Int,
		cells = 1 : size(model, deg, 1);
		atol = 1e-10
	)::Array{Array{Float64,2},1}

Evaluates the bounding boxes related to each `deg`-cell of the `model`

The function builds a vector of matrices containing the two distincive vertices
of the axis-coherent hyper-cube that contains each `deg`-cell.

If `cells` is specified, then only those bounding boxes are evaluated

Keyword parameter `atol` specifies the tollerance added to the bounding box

See also: [`CAGD.spatialIndex`](@ref)
It uses:
---
# Examples
```jldoctest
julia> model = CAGD.Model(hcat([
		[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.5, 1.0], [0.5, 1.0, 0.5], [1.0, 1.0, 0.5]
	]...));
julia> CAGD.addModelCells!(model, 1, [[1, 2], [2, 5], [3, 4], [4, 5]]);

julia> CAGD.getModelBoundingBoxes(model,1)
4-element Array{Array{Float64,2},1}:
4-element Array{Array{Float64,2},1}:
 [0.0 1.0; 0.0 1.0; 0.0 0.0]
 [0.0 1.0; 1.0 1.0; 0.0 0.5]
 [0.0 0.5; 0.5 1.0; 0.5 1.0]
 [0.5 1.0; 1.0 1.0; 0.5 0.5]
```
"""

function getModelBoundingBoxes(
		model::CAGD.Model,
		deg::Int,
		cells = 1 : size(model, deg, 1);
		atol = 1e-10
	)::Array{Array{Float64,2},1}

	function findBBox(pts::Lar.Points)::Tuple{Array{Float64,2},Array{Float64,2}}
		minimum = mapslices(x->min(x...), pts, dims=2) .- 2*atol
		maximum = mapslices(x->max(x...), pts, dims=2) .+ 2*atol
		return minimum, maximum
	end

	cellpts = [ hcat(CAGD.getModelCellVertices(model, deg, c)...)  for c in cells ]
	return [hcat(findBBox(c)...) for c in cellpts]
end

#-------------------------------------------------------------------------------
#   MODEL BOUNDING BOXES INCLUSION
#-------------------------------------------------------------------------------

"""
	isBoundingBoxIncluded(bbox1::Array{Float64,2}, bbox2::Array{Float64,2})::Bool

Checks whether `bbox1` is totally included in `bbox2`
"""
function isBoundingBoxIncluded(bbox1::Array{Float64,2}, bbox2::Array{Float64,2})::Bool
	dim = size(bbox1, 1)
	for d = 1 : dim
		(bbox1[d, 1] >= bbox2[d, 1]) & (bbox1[d, 2] <= bbox2[d, 2]) || return false
	end
	return true
end

#-------------------------------------------------------------------------------
#   SPATIAL_INDEXING METHODS
#	 buildBBIntervalTree
#	 extractBoxCovering
#-------------------------------------------------------------------------------

"""
	buildBBIntervalTree(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::IntervalTrees.IntervalMap{Float64, Array}

Build the Interval Tree of a set of bounding boxes along the `dim` coordinate.
"""
function buildBBIntervalTree(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::IntervalTrees.IntervalMap{Float64, Array}

	# Build the dictionary related to the `dim`-coordinate
	boxDict = OrderedDict{Array{Float64,1}, Array{Int,1}}()
	for (h, box) in enumerate(bboxes)
		key = box[dim, :]
		if haskey(boxDict, key) == false
			boxDict[key] = [h]
		else
			push!(boxDict[key], h)
		end
	end

	# Build the related Interval Tree
	iT = IntervalTrees.IntervalMap{Float64, Array}()
	for (key, boxSet) in boxDict
		iT[tuple(key...)] = boxSet
	end

	return iT
end

"""
	extractBoxCovering(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::Array{Array{Int,1},1}

Evaluates possible `bboxes` intersections over the `dim`-th coordinate

The function evaluates the pairwise intersection between the `dim` coordinate
of each `bboxes` couple. The operation is performed via Interval Tree.
"""
function extractBoxCovering(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::Array{Array{Int,1},1}

	# Build the related Interval Tree
	tree = CAGD.buildBBIntervalTree(bboxes, dim)

	# Build the covers
	covers = [[] for k = 1 : length(bboxes)]
	for (i, bbox) in enumerate(bboxes)
		extent = bboxes[i][dim, :]
		iterator = IntervalTrees.intersect(tree, tuple(extent...))
		for x in iterator
			append!(covers[i], x.value)
		end
	end

	return covers
end

#-------------------------------------------------------------------------------
#   SPATIAL_INDEXING
#-------------------------------------------------------------------------------

"""
	spaceIndex(model::CAGD.Model, dim::Int)::Array{Array{Int,1},1}

Generation of *space indexes* for all ``(d-1)``-dim cell members of `model`.

*Spatial index* made by ``d`` *interval-trees* on
bounding boxes of ``sigma in S_{d−1}``. Spatial queries solved by
intersection of ``d`` queries on IntervalTrees generated by
bounding-boxes of geometric objects (LAR cells).

The return value is an array of arrays of `int`s, indexing cells whose
containment boxes are intersecting the containment box of the first cell.
According to Hoffmann, Hopcroft, and Karasick (1989) the worst-case complexity of
Boolean ops on such complexes equates the total sum of such numbers.

# Examples 2D

```
julia> model = CAGD.Model(hcat([[0.,0],[1,0],[1,1],[0,1],[2,1]]...));

julia> CAGD.addModelCells!(model, 1, [[1,2],[2,3],[3,4],[4,1],[1,5]]);

julia> Sigma = CAGD.spaceIndex(model, 1)
5-element Array{Array{Int64,1},1}:
 [4, 5, 2]
 [1, 3, 5]
 [4, 5, 2]
 [1, 3, 5]
 [4, 1, 3, 2]
```

# Example 3D

```julia
julia> model = CAGD.Model([
		0.0 1.0 0.0 0.0 0.0 0.0 1.0
		0.0 0.0 1.0 0.0 0.0 0.0 0.0
		0.0 0.0 0.0 1.0 2.0 3.0 2.0
	])
julia> CAGD.addModelCells!(model, 1, [
		[1,3],[1,2],[2,3],[3,4],[2,4],[1,4],[4,5],[5,7],[6,7],[5,6]
	])
julia> CAGD.addModelCells!(model, 2,[
		[1,2,3],[1,4,6],[2,5,6],[3,4,5],[8,9,10]
	])

julia> CAGD.spaceIndex(model, 1)
10-element Array{Array{Int64,1},1}:
 [4, 6, 2, 3, 5]
 [1, 4, 6, 3, 5]
 [1, 4, 6, 2, 5]
 [1, 6, 7, 2, 3, 5]
 [1, 4, 6, 7, 2, 3]
 [1, 4, 7, 2, 3, 5]
 [4, 6, 10, 5, 8, 9]
 [7, 10, 9]
 [7, 10, 8]
 [7, 8, 9]

CAGD.spaceIndex(model, 2)
5-element Array{Array{Int64,1},1}:
 [2, 3, 4]
 [1, 3, 4]
 [2, 1, 4]
 [2, 1, 3]
 []

```
"""
function spaceIndex(model::CAGD.Model, dim::Int)::Array{Array{Int,1},1}

	bboxes = CAGD.getModelBoundingBoxes(model, dim)
	# Build Box Coverings for each dimension and intersecting them
	covers = CAGD.extractBoxCovering(bboxes, 1)
	for d = 2 : length(model)
		dcovers = CAGD.extractBoxCovering(bboxes, d)
		covers  = [intersect(pair...) for pair in zip(covers, dcovers)]
	end

	# Remove each cell from its cover
	for k = 1 : length(covers)
		covers[k] = setdiff(covers[k],[k])
	end

	return covers
end


text:
__precompile__()
module Zeros

import Base: +, -, *, /, <, >, <=, >=, fma, muladd, mod, rem, modf,
     ldexp, copysign, flipsign, sign, round, floor, ceil, trunc,
     promote_rule, convert, show, significand, string
     AbstractFloat, Integer, Complex

export Zero, testzero, zero!

 "A type that stores no data, and holds the value zero."
struct Zero <: Integer
end

promote_rule(::Type{Zero}, ::Type{T}) where {T<:Number} = T
promote_rule(::Type{Zero}, ::Type{Zero}) = Float64

convert(::Type{Zero}, ::Zero) = Zero()
convert(::Type{T}, ::Zero) where {T<:Number} = zero(T)
convert(::Type{Zero}, x::T) where {T<:Number} = Zero(x)

if VERSION < v"0.7-"
    Zero(x::Number) = iszero(x) ? Zero() : throw(InexactError())

    # Disambiguation needed for Julia 0.6
    convert(::Type{T}, ::Zero) where {T<:Real} = zero(T)
    convert(::Type{BigInt}, ::Zero) = zero(BigInt)
    convert(::Type{BigFloat}, ::Zero) = zero(BigFloat)
    convert(::Type{Float16}, ::Zero) = zero(Float16)
    convert(::Type{Complex{T}}, ::Zero) where {T<:Real} = zero(Complex{T})
else
    Zero(x::Number) = iszero(x) ? Zero() : throw(InexactError(:Zero, Zero, x))
end

AbstractFloat(::Zero) = 0.0

Complex(x::Real, ::Zero) = x

# Loop over types in order to make methods specific enough to avoid ambiguities.
for T in (Number, Real, Integer, Complex, Complex{Bool})
    Base.:+(x::T,::Zero) = x
    Base.:+(::Zero,x::T) = x
    Base.:-(x::T,::Zero) = x
    Base.:-(::Zero,x::T) = -x
    Base.:*(::T,::Zero) = Zero()
    Base.:*(::Zero,::T) = Zero()
    Base.:/(::Zero, ::T) = Zero()
    Base.:/(::T, ::Zero) = throw(DivideError())
end

+(::Zero,::Zero) = Zero()
-(::Zero,::Zero) = Zero()
*(::Zero,::Zero) = Zero()
/(::Zero, ::Zero) = throw(DivideError())

<(::Zero,::Zero) = false
<=(::Zero,::Zero) = true

ldexp(::Zero, ::Integer) = Zero()
copysign(::Zero,::Real) = Zero()
flipsign(::Zero,::Real) = Zero()
modf(::Zero) = (Zero(), Zero())

# Alerady working due to default definitions:
# ==, !=, abs, isinf, isnan, isinteger, isreal, isimag,
# signed, unsigned, float, big, complex, isodd, iseven

for op in [:(-), :sign, :round, :floor, :ceil, :trunc, :significand]
  @eval $op(::Zero) = Zero()
end

# Avoid promotion of triplets
for op in [:fma :muladd]
    @eval $op(::Zero, ::Zero, ::Zero) = Zero()
    for T in (Real, Integer)
        @eval $op(::Zero, ::$T, ::Zero) = Zero()
        @eval $op(::$T, ::Zero, ::Zero) = Zero()
        @eval $op(::Zero, x::$T, y::$T) = convert(promote_type(typeof(x),typeof(y)),y)
        @eval $op(x::$T, ::Zero, y::$T) = convert(promote_type(typeof(x),typeof(y)),y)
        @eval $op(x::$T, y::$T, ::Zero) = x*y
    end
end

for op in [:mod, :rem]
  @eval $op(::Zero, ::Real) = Zero()
end

# Display Zero() as `0̸` (unicode "zero with combining long solidus overlay")
# so that it looks slightly different from `0`.
show(io::IO, ::Zero) = print(io, "0̸")
string(z::Zero) = Base.print_to_string(z)

# Display Array{Zero} without printing the actual zeros.
show(io::IO, x::Array{Zero,N}) where {N} = show(io, MIME"text/plain", x)
function show(io::IO, ::MIME"text/plain", x::Array{Zero,N}) where {N}
    join(io, map(dec, size(x)), "x")
    print(io, " Array{Zero,", N, "}")
end

# This function is intentionally not type-stable.
"Convert to Zero() if zero. (Use immediately before calling a function.)"
testzero(x::Number) = x==zero(x) ? Zero() : x

"Fill an array with zeros."
zero!(a::Array{T}) where T = fill!(a, Zero())

end # module


text:
import Adapt

using StaticArrays

dummy() = return

@testset "@veda" begin

@test_throws UndefVarError @veda undefined()
@test_throws MethodError @veda dummy(1)

@testset "low-level interface" begin
    k = vefunction(dummy)
    k()
    #k(; items=1)
end

#@testset "launch configuration" begin
#    @veda dummy()
#
#    @veda items=1 dummy()
#    @veda items=(1,1) dummy()
#    @veda items=(1,1,1) dummy()
#
#    @veda groups=1 dummy()
#    @veda groups=(1,1) dummy()
#    @veda groups=(1,1,1) dummy()
#end

@testset "launch=false" begin
    k = @veda launch=false dummy()
    k()
    #k(; items=1)
end

@testset "inference" begin
    foo() = @veda dummy()
    @inferred foo()

    # with arguments, we call kernel_convert
    kernel(a) = return
    bar(a) = @veda kernel(a)
    @inferred bar(VEArray([1]))
end

@testset "reflection" begin
    VectorEngine.code_lowered(dummy, Tuple{})
    VectorEngine.code_typed(dummy, Tuple{})
    VectorEngine.code_warntype(devnull, dummy, Tuple{})
    VectorEngine.code_llvm(devnull, dummy, Tuple{})
    VectorEngine.code_native(devnull, dummy, Tuple{})

    @device_code_lowered @veda dummy()
    @device_code_typed @veda dummy()
    @device_code_warntype io=devnull @veda dummy()
    @device_code_llvm io=devnull @veda dummy()
    @device_code_native io=devnull @veda dummy()

    mktempdir() do dir
        @device_code dir=dir @veda dummy()
    end

    @test_throws ErrorException @device_code_lowered nothing

    # make sure kernel name aliases are preserved in the generated code
    @test occursin("julia_dummy", sprint(io->(@device_code_llvm io=io optimize=false @veda dummy())))
    @test occursin("julia_dummy", sprint(io->(@device_code_llvm io=io @veda dummy())))
    @test occursin("julia_dummy", sprint(io->(@device_code_native io=io @veda dummy())))

    # make sure invalid kernels can be partially reflected upon
    let
        invalid_kernel() = throw()
        @test_throws VectorEngine.KernelError @veda invalid_kernel()
        @test_throws VectorEngine.KernelError @grab_output @device_code_warntype @veda invalid_kernel()
        out, err = @grab_output begin
            try
                @device_code_warntype @veda invalid_kernel()
            catch
            end
        end
        @test occursin("Body::Union{}", err)
    end

    let
        range_kernel() = (0.0:0.1:100.0; nothing)

        @test_throws VectorEngine.InvalidIRError @veda range_kernel()
    end

    # set name of kernel
    @test occursin("julia_mykernel", sprint(io->(@device_code_llvm io=io begin
        k = vefunction(dummy, name="mykernel")
        k()
    end)))
end

@testset "external kernels" begin
    @eval module KernelModule
        export external_dummy
        external_dummy() = return
    end
    import ...KernelModule
    @veda KernelModule.external_dummy()
    @eval begin
        using ...KernelModule
        @veda external_dummy()
    end

    @eval module WrapperModule
        using VectorEngine
        @eval dummy() = return
        wrapper() = @veda dummy()
    end
    WrapperModule.wrapper()
end

@testset "calling device function" begin
    @noinline child(i) = sink(i)
    function parent()
        child(1)
        return
    end

    VectorEngine.@veda parent()
end

@testset "varargs" begin
    function kernel(r, args...)
        r[] = args[2]
        return
    end

    r = Ref{Int}(-1)

    @veda kernel(r, 1, 2, 3)
    synchronize()
    @test r[] == 2
end

end

############################################################################################

@testset "argument passing" begin

dims = (16, 16)
len = prod(dims)

@testset "manually allocated" begin
    function kernel(input, output)
        for i = 1:length(input)
            output[i] = input[i]
        end
        return
    end

    input = round.(rand(Float32, dims) * 100)
    output = similar(input)

    input_dev = VEArray(input)
    output_dev = VEArray(output)

    @veda kernel(input_dev, output_dev)
    @test input ≈ Array(output_dev)
end

@testset "scalar through single-value array" begin
    function kernel(a, x)
        max = length(a)
        for i = 1:max
            if i == max
                x[] = a[i]
            end
        end
        return
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = VEArray(arr)
    val_dev = VEArray(val)

    @veda kernel(arr_dev, val_dev)
    @test arr[dims...] ≈ Array(val_dev)[1]
end

@testset "scalar through single-value array, using device function" begin
    @noinline child(a, i) = a[i]
    function parent(a, x)
        max = length(a)
        for i = 1:max
            if i == max
                x[] = child(a, i)
            end
        end
        return
    end

    arr = round.(rand(Float32, dims) * 100)
    val = [0f0]

    arr_dev = VEArray(arr)
    val_dev = VEArray(val)

    @veda parent(arr_dev, val_dev)
    @test arr[dims...] ≈ Array(val_dev)[1]
end

@testset "tuples" begin
    # issue #7: tuples not passed by pointer

    function kernel(keeps, out)
        if keeps[1]
            out[] = 1
        else
            out[] = 2
        end
        return
    end

    keeps = (true,)
    d_out = VEArray(zeros(Int))

    @veda kernel(keeps, d_out)
    @test Array(d_out)[] == 1
end

@testset "ghost function parameters" begin
    # bug: ghost type function parameters are elided by the compiler

    len = 60
    a = rand(Float32, len)
    b = rand(Float32, len)
    c = similar(a)

    d_a = VEArray(a)
    d_b = VEArray(b)
    d_c = VEArray(c)

    @eval struct ExecGhost end

    function kernel(ghost, a, b, c)
        for i = 1:length(a)
            c[i] = a[i] + b[i]
        end
        return
    end
    @veda kernel(ExecGhost(), d_a, d_b, d_c)
    @test a+b == Array(d_c)

    # bug: ghost type function parameters confused aggregate type rewriting
    function kernel(ghost, out, aggregate)
        for i = 1:length(out)
            out[i] = aggregate[1]
        end
        return
    end
    @veda kernel(ExecGhost(), d_c, (42,))

    @test all(val->val==42, Array(d_c))
end

@testset "immutables" begin
    # issue #15: immutables not passed by pointer

    function kernel(ptr, b)
        ptr[] = imag(b)
        return
    end

    arr = VEArray(zeros(Float32))
    x = ComplexF32(2,2)

    @veda kernel(arr, x)
    @test Array(arr)[] == imag(x)
end

@testset "automatic recompilation" begin
    arr = VEArray(zeros(Int))

    function kernel(ptr)
        ptr[] = 1
        return
    end

    @veda kernel(arr)
    @test Array(arr)[] == 1

    function kernel(ptr)
        ptr[] = 2
        return
    end

    @veda kernel(arr)
    @test Array(arr)[] == 2
end

@testset "automatic recompilation (bis)" begin
    arr = VEArray(zeros(Int))

    @eval doit(ptr) = ptr[] = 1

    function kernel(ptr)
        doit(ptr)
        return
    end

    @veda kernel(arr)
    @test Array(arr)[] == 1

    @eval doit(ptr) = ptr[] = 2

    @veda kernel(arr)
    @test Array(arr)[] == 2
end

@testset "non-isbits arguments" begin
    function kernel1(T, i)
        sink(i)
        return
    end
    @veda kernel1(Int, 1)

    function kernel2(T, i)
        sink(unsafe_trunc(T,i))
        return
    end
    @veda kernel2(Int, 1.)
end

@testset "splatting" begin
    function kernel(out, a, b)
        out[] = a+b
        return
    end

    out = [0]
    out_dev = VEArray(out)

    @veda kernel(out_dev, 1, 2)
    @test Array(out_dev)[1] == 3

    all_splat = (out_dev, 3, 4)
    @veda kernel(all_splat...)
    @test Array(out_dev)[1] == 7

    partial_splat = (5, 6)
    @veda kernel(out_dev, partial_splat...)
    @test Array(out_dev)[1] == 11
end

@testset "return value by ref" begin
    x = Ref{Int}(123)

    function kernel(r)
        r[] += 10
        return
    end

    @synchronize @veda kernel(x)
    @test x[] == 133
end

@testset "pass and return mutabale struct" begin
    mutable struct xm
        x::Int32
        m::Int64
    end

    s = xm(10, 10)

    function kernel(s)
        s.x += 1
        s.m -= 1
        return
    end

    @synchronize @veda kernel(s)
    @test s.x == 11 && s.m == 9
end

@testset "object invoke" begin
    # this mimics what is generated by closure conversion
    # FAIL on VE, kills VEO context

    @eval struct KernelObject{T} <: Function
        val::T
    end

    function (self::KernelObject)(a)
        a[] = self.val
        return
    end

    function outer(a, val)
       inner = KernelObject(val)
       @veda inner(a)
    end

    a = [1.]
    a_dev = VEArray(a)

    outer(a_dev, 2.)

    @test Array(a_dev) ≈ [2.]
end

@testset "closures" begin
    # FAIL on VE, kills VEO context
    function outer(a_dev, val)
       function inner(a)
            # captures `val`
            a[] = val
            return
       end
       @veda inner(a_dev)
    end

    a = [1.]
    a_dev = VEArray(a)

    outer(a_dev, 2.)

    @test Array(a_dev) ≈ [2.]
end

@testset "conversions" begin
    @eval struct Host   end
    @eval struct Device end

    Adapt.adapt_storage(::VectorEngine.Adaptor, a::Host) = Device()

    Base.convert(::Type{Int}, ::Host)   = 1
    Base.convert(::Type{Int}, ::Device) = 2

    out = [0]

    # convert arguments
    out_dev = VEArray(out)
    let arg = Host()
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg)
            return
        end
        @veda kernel(arg, out_dev)
        @test Array(out_dev) ≈ [2]
    end

    # convert tuples
    out_dev = VEArray(out)
    let arg = (Host(),)
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg[1])
            return
        end
        @veda kernel(arg, out_dev)
        @test Array(out_dev) ≈ [2]
    end

    # convert named tuples
    out_dev = VEArray(out)
    let arg = (a=Host(),)
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg.a)
            return
        end
        @veda kernel(arg, out_dev)
        @test Array(out_dev) ≈ [2]
    end

    # don't convert structs
    out_dev = VEArray(out)
    @eval struct Nested
        a::Host
    end
    let arg = Nested(Host())
        @test Array(out_dev) ≈ [0]
        function kernel(arg, out)
            out[] = convert(Int, arg.a)
            return
        end
        @veda kernel(arg, out_dev)
        @test Array(out_dev) ≈ [1]
    end
end

@testset "argument count" begin
    val = [0]
    val_dev = VEArray(val)
    for i in (1, 10, 20, 34)
        variables = ('a':'z'..., 'A':'Z'...)
        params = [Symbol(variables[j]) for j in 1:i]
        # generate a kernel
        body = quote
            function kernel(arr, $(params...))
                arr[] = $(Expr(:call, :+, params...))
                return
            end
        end
        eval(body)
        args = [j for j in 1:i]
        call = Expr(:call, :kernel, val_dev, args...)
        vedacall = :(@veda $call)
        eval(vedacall)
        @test Array(val_dev)[1] == sum(args)
    end
end

@testset "keyword arguments" begin
    @eval inner_kwargf(foobar;foo=1, bar=2) = nothing

    @veda (()->inner_kwargf(42;foo=1,bar=2))()

    @veda (()->inner_kwargf(42))()

    @veda (()->inner_kwargf(42;foo=1))()

    @veda (()->inner_kwargf(42;bar=2))()

    @veda (()->inner_kwargf(42;bar=2,foo=1))()
end

@testset "captured values" begin
    function f(capture::T) where {T}
        function kernel(ptr)
            ptr[] = capture
            return
        end

        arr = VEArray(zeros(T))
        @veda kernel(arr)

        return Array(arr)[1]
    end

    using Test
    @test f(1) == 1
    @test f(2) == 2
end

end

############################################################################################

#?#@testset "#55: invalid integers created by alloc_opt" begin
#?#    function f(a)
#?#        x = SVector(0f0, 0f0)
#?#        v = MVector{3, Float32}(undef)
#?#        for (i,_) in enumerate(x)
#?#            v[i] = 1.0f0
#?#        end
#?#        a[1] = v[1]
#?#        return nothing
#?#    end
#?#    @veda f(VEArray(zeros(1)))
#?#end

############################################################################################


text:
using FactorizationMachines
using Base.Test

function main()
    train_filename = ASCIIString(ARGS[1])
    test_filename = ASCIIString(ARGS[2])

    info("Running benchmark tests ...")
    info("Train filename: $train_filename")
    info("Test filename : $test_filename")

    (X_train, y_train) = @time read_libsvm(train_filename)
    info("Train dim: $(size(X_train))")
    (X_test, y_test) = @time read_libsvm(test_filename)
    info("Test dim: $(size(X_test))")

    info("Training model...")
    fm = @time train(X_train, y_train,
        method = Methods.sgd(num_epochs = Int(10), alpha = 0.1),
        model_params = Models.gauss(num_factors = 4))

    info("Prediction over test dataset...")
    p = @time predict(fm, X_test)

    rmse = sqrt(sum((p - y_test).^2) / length(p))
    info("RMSE on test data: $rmse")
end

main()


text:
using GridWorlds
using Test
using Random

ENVS = [EmptyGridWorld, FourRooms, GoToDoor]
ACTIONS = [TURN_LEFT, TURN_RIGHT, MOVE_FORWARD]

@testset "GridWorlds.jl" begin
    for env in ENVS
        @testset "$(env)" begin
            w = env()
            @test typeof(w.agent_pos) == CartesianIndex{2}
            @test typeof(w.agent.dir) <: LRUD
            @test size(w.world.world, 1) == length(w.world.objects)
            @test 1 ≤ w.agent_pos[1] ≤ size(w.world.world, 2)
            @test 1 ≤ w.agent_pos[2] ≤ size(w.world.world, 3)

            for _=1:1000
                w = w(rand(ACTIONS))
                @test 1 ≤ w.agent_pos[1] ≤ size(w.world.world, 2)
                @test 1 ≤ w.agent_pos[2] ≤ size(w.world.world, 3)
                @test w.world[WALL, w.agent_pos] == false
                view = get_agent_view(w)
                @test typeof(view) <: BitArray{3}
                @test size(view,1) == length(w.world.objects)
                @test size(view,2) == 7
                @test size(view,3) == 7
            end
        end
    end
end


text:
let
    #= line 1 =#
    var"|Record|b_op_Getter"::Function = function op_Getter(x)
                #= none:1 =#
                #= none:4 =#
                function (a,)
                    #= none:4 =#
                    #= none:7 =#
                    x(a)
                end
            end::Function
    var"|Record|c_tupleFirst"::Any = function (x,)
                #= none:1 =#
                #= none:4 =#
                #= none:4 =# @inbounds x[1]
            end::Any
    var"|Record|d_unit"::Symbol = Symbol("unit@9_|Record|e")::Symbol
    var"|Record|f_println"::Function = println::Function
    (var"|Record|f_println"::Function)((op_Getter(var"|Record|c_tupleFirst"::Any))((1::Int64, "3"::String)::Tuple)::Int64)::Any
    (var"|Record|f_println"::Function)((op_Getter(var"|Record|c_tupleFirst"::Any))(((1::Int64, "3"::String)::Tuple, 2::Int64)::Tuple)::Tuple)::Any
    (var"|Record|f_println"::Function)((op_Getter(var"|Record|c_tupleFirst"::Any))((1::Int64, "3"::String)::Tuple)::Int64)::Any
    (var"|Record|f_println"::Function)((op_Getter(var"|Record|c_tupleFirst"::Any))(((1.0::Float64, 2.0::Float64)::Tuple, 4.0::Float64)::Tuple)::Tuple)::Any
    var"|Record|g"::Function = function (_,)
            #= ../examples/records.mlfsf:1 =#
            var"|Record|f_println"
        end
    var"|Record|h"::Function = function (_,)
            #= ../examples/records.mlfsf:1 =#
            Symbol("unit@9_|Record|e")
        end
    var"|Record|i"::Function = function (_,)
            #= ../examples/records.mlfsf:1 =#
            var"|Record|c_tupleFirst"
        end
    var"|Record|j"::Function = function (_,)
            #= ../examples/records.mlfsf:1 =#
            var"|Record|b_op_Getter"
        end
    var"|Record|"::Any = Symbol("Module Record")
    Record::Any = Symbol("namespace {Module Record, {{unit, |Record|d_unit, Type unit@9_|Record|e}, {println, |Record|f_println, forall a@6600224964343670762.a@6600224964343670762->unit@9_|Record|e}, {tupleFirst, |Record|c_tupleFirst, implicit[forall a@3062623644099412066 b@7339522248752322263.field {{a@3062623644099412066, b@7339522248752322263}, first, a@3062623644099412066}]}, {op_Getter, |Record|b_op_Getter, forall a@12824608043531541150 b@3103303379621457644 c@6082679040560098397.(field {a@12824608043531541150, b@3103303379621457644, c@6082679040560098397})->a@12824608043531541150->c@6082679040560098397}}}")
    0
end


text:
function initialize_search_tasks!(tuning_run::Run)
    instance_id    = 1
    results        = RemoteChannel[]

    tuning_run.cost_values    = zeros(tuning_run.cost_evaluations)
    tuning_run.starting_cost  = tuning_run.measurement_method(tuning_run,
                                                              tuning_run.starting_point)
    initial_result            = Result("Initialize",
                                       tuning_run.starting_point,
                                       tuning_run.starting_point,
                                       tuning_run.starting_cost,
                                       1, 1, 1, false)

    for i = 1:size(tuning_run.methods, 1)
        for j = 1:tuning_run.methods[i, 2]
            push!(results, RemoteChannel(()->ResultChannel(deepcopy(initial_result))))

            channel = results[instance_id]
            @spawn eval(tuning_run.methods[i, 1])(deepcopy(tuning_run), channel)
            instance_id += 1
        end
    end

    results
end


text:
# SteinPosterior
# An abstract subtype of SteinDistribution representing a posterior
# distribution p(x) = pi(x | y_1,...,y_L) induced by a likelihood
# pi(y_1,...,y_L|x) over datapoints y_L and a prior pi(x)
abstract SteinPosterior <: SteinDistribution

# Returns the number of datapoints on which this posterior is conditioned
function numdatapoints(d::SteinPosterior)
    size(d.y, 1)
end

# Returns the mini-batch approximation
# log pi(x) + (L/length(idx)) log pi(x | d.y[idx])
# to the log posterior density
function logdensity(d::SteinPosterior,
                    x::Array{Float64, 1};
                    idx=1:numdatapoints(d))
    batchratio = numdatapoints(d)/length(idx)
    logprior(d,x) + batchratio * loglikelihood(d,x;idx=idx)
end

# For each sample point xs[i,:], returns the mini-batch approximation to the log posterior density
function logdensity(d::SteinPosterior,
                    xs::Array{Float64, 2};
                    idx=1:numdatapoints(d))
    logdensities = zeros(size(xs,1))
    for i in 1:size(x,1)
        logdensities[i] = logdensity(d, vec(xs[i,:]); idx=idx)
    end
    logdensities
end

# A generic to make applying logprior to each row.
function logprior(d::SteinPosterior,
                  xs::Array{Float64, 2})
    logpriors = Array(Float64, size(xs,1))
    for i in size(xs,1)
        logpriors[i] = logprior(d, vec(xs[i,:]))
    end
    logpriors
end

# A generic to make applying loglikelihood to each row.
function loglikelihood(d::SteinPosterior,
                       xs::Array{Float64, 2};
                       idx=1:numdatapoints(d))
    loglikelihoods = Array(Float64, size(xs,1))
    for i in size(xs,1)
        loglikelihoods[i] = loglikelihood(d, vec(xs[i,:]); idx=idx)
    end
    loglikelihoods
end

# For each sample point xs[i,:] returns log likelihood gradient
# grad_{x} log pi(d.y[idx] | xs[i,:])
# over the batch of datapoints d.y[idx]
function gradloglikelihood(d::SteinPosterior,
                           xs::Array{Float64, 2};
                           idx=1:numdatapoints(d))
    q = size(xs, 1)
    if q < 1
        return []
    end
    gradloglikelihoods = zeros(size(xs))
    for i=1:q
        gradloglikelihoods[i,:] = vec(gradloglikelihood(d, vec(xs[i,:]); idx=idx))
    end
    gradloglikelihoods
end

# Returns the mini-batch approximation
# grad_{x} log pi(x) + (L/length(idx)) grad_{x} log pi(x | d.y[idx])
# to the gradient of the log posterior density
function gradlogprior(d::SteinPosterior,
                      x::Array{Float64, 2})
    gradlogpriors = zeros(size(x))
    for i in 1:size(x,1)
        gradlogpriors[i,:] = gradlogprior(d, vec(x[i,:]))
    end
    gradlogpriors
end

# Returns the mini-batch approximation
# grad_{x} log pi(x) + (L/length(idx)) grad_{x} log pi(x | d.y[idx])
# to the gradient of the log posterior density
function gradlogdensity(d::SteinPosterior,
                        x::Array{Float64, 1};
                        idx=1:numdatapoints(d))
    batchratio = numdatapoints(d)/length(idx)
    gradlogprior(d,x) + batchratio .* gradloglikelihood(d, x; idx=idx)
end

# For each sample point xs[i,:] returns the mini-batch approximation
# to the gradient of the log posterior density
function gradlogdensity(d::SteinPosterior,
                        xs::Array{Float64, 2};
                        idx=1:numdatapoints(d))
    q = size(xs, 1)
    if q < 1
        return []
    end
    gradlogdensities = zeros(size(xs))
    for i=1:q
        gradlogdensities[i,:] = vec(gradlogdensity(d, vec(xs[i,:]); idx=idx))
    end
    gradlogdensities
end


text:
using Printf
using ArgParse
using HTTP
using JSON2
using Jtac

# -------- List of supported games ------------------------------------------- #

const GAMES = Dict(
    "tictactoe" => [:TicTacToe]
  , "mnkgame"   => [:MNKGame, Int, Int, Int]
  , "metatac"   => [:MetaTac]
  , "morris"    => [:Morris]
)

# -------- En-/decoding visual / renderable content -------------------------- #

include("rencoder.jl")

# -------- En-/decoding game states ------------------------------------------ #

struct GameState
  typ :: String
  params :: Vector{String}
  game :: Game
end

function GameState(game :: Game)

  typ = string(typeof(game))
  ff = findfirst('{', typ)
  fl = findlast('}', typ)

  if !isnothing(ff) && !isnothing(fl)
    params = split(typ[ff+1:fl-1], ",")
    typ = lowercase(typ[1:ff-1])
  else
    typ = lowercase(typ)
    params = String[]
  end

  GameState(typ, params, game)
end

function JSON2.read(io :: IO, :: Type{GameState})

  fields = JSON2.read(io)
  typ = fields.typ

  gametype = GAMES[typ][1]
  paramtypes = GAMES[typ][2:end]
  params = map(parse, paramtypes, fields.params)

  cexpr = isempty(params) ? gametype : Expr(:curly, gametype, params...)
  game = eval(cexpr)()

  for field in fieldnames(typeof(game))
    value = convert(getfield(game, field) |> typeof, getfield(fields.game, field))
    setfield!(game, field, value)
  end

  GameState(game)

end

# -------- Requests and Responses -------------------------------------------- #

struct TurnRequest
  game :: GameState
  action :: Int
end

TurnRequest(; game, action) = TurnRequest(game, action)
JSON2.@format TurnRequest keywordargs

struct ApplyRequest
  game :: GameState
  model :: String
  power :: Int
  temperature :: Float32
  exploration :: Float32
end

function ApplyRequest(; game, model, power, temperature, exploration)
  ApplyRequest(game, model, power, temperature, exploration)
end

JSON2.@format ApplyRequest keywordargs

struct ApplyResponse
  value :: Float32
  policy :: Vector{Float32}
  features :: Dict{String,Vector{Float32}}
end

# -------- Auxiliary Functions ----------------------------------------------- #

wrap(f, args...) = HTTP.Response(200, JSON2.write(f(args...)))

function load_models(dir, gpu)

  # Get all files/directories in dir
  files = readdir(dir)

  # Extract the files with '.jtm' ending and load them as async (gpu) models
  models = map(filter(x->length(x) > 4 && x[end-3:end] == ".jtm", files)) do fname
    name = fname[1:end - 4] # strip '.jtm'
    path = joinpath(dir, name)
    model = gpu ? load_model(path) |> to_gpu : load_model(path)
    name => Async(model)
  end

  Dict{String,Model}( "rollout" => RolloutModel()
                    , "random" => RandomModel()
                    , models... )

end

# -------- Handling GET Requests --------------------------------------------- #

function GET_index(req::HTTP.Request, index)
  if index == ""
    content ="Welcome the the Jtac AI server. No index file was provided"
  else
    content = read(index)
  end
  HTTP.Response(200, content)
end

function GET_reload(req::HTTP.Request, dir, gpu, models)

  # Remove previous models
  empty!(models)

  # Add new ones
  # TODO: We should be more efficient here and only look for newly added
  # / deleted models
  merge!(models, load_models(dir, gpu))

  "Successfully loaded $(length(models)) models."

end

function GET_models(req::HTTP.Request, models)
  map(collect(models)) do (name, model)
    Dict( "name"     => name, "game"     => string(gametype(model))
        , "params"   => sum([0; length.(Knet.params(model))])
        , "features" => Jtac.feature_name.(Jtac.features(model)))
  end
end

function GET_games(req::HTTP.Request)
  map(keys(GAMES), values(GAMES)) do k, v
    Dict("typ" => k, "params" => v[2:end])
  end
end

function GET_create(req::HTTP.Request)
  typ = HTTP.URIs.splitpath(req.target)[3]
  params = HTTP.URIs.splitpath(req.target)[4:end]

  gametype = GAMES[typ][1]
  paramtypes = GAMES[typ][2:end]
  params = map(parse, paramtypes, params)

  cexpr = isempty(params) ? gametype : Expr(:curly, gametype, params...)
  GameState(eval(cexpr)())
end

function GET_dummyimage(req::HTTP.Request)
  game = TicTacToe()
  for i in rand(1:5) random_turn!(game) end
  Image(game, rand(Float32, 9), 0.5f0) |> Typed
end

function POST_turn(req::HTTP.Request)
  req = JSON2.read(IOBuffer(HTTP.payload(req)), TurnRequest)
  apply_action!(req.game.game, req.action)
  GameState(req.game.game)
end

function POST_apply(req::HTTP.Request, models)

  # Parse target and json and sent via POST
  paths = HTTP.URIs.splitpath(req.target)

  req = JSON2.read(IOBuffer(HTTP.payload(req)), ApplyRequest)

  # Get model, game, power
  model  = models[req.model]
  game   = req.game.game
  power  = req.power

  # A value power <= 1 means that we pick the IntuitionPlayer
  if power <= 1
    player = IntuitionPlayer(model, temperature = req.temperature)
  else
    player = MCTSPlayer( model
                       , power = power
                       , temperature = req.temperature
                       , exploration = req.exploration )
  end

  # TODO: Get (improved) value of state from player!
  v, _, f = apply_features(base_model(model), game) # get value and features
  p = think(player, game)                # get improved policy

  # Bring features in nicer format as dictionary
  feats = Jtac.features(base_model(model))
  fdict = Dict{String,Vector{Float32}}(
            Jtac.feature_name(feat) => Jtac.feature_conv(feat, f[idx])
            for (feat, idx) in zip(feats, Jtac.feature_indices(feats, typeof(game))) )

  # If a visual response is desired, return the corresponding image data
  if length(paths) == 3 && paths[3] == "visual"

    Image(game, p, v) |> Typed

  # Else, send value/policy/feature data directly
  else

    ApplyResponse(v, p, fdict)

  end

end

# -------- Main entry function ----------------------------------------------- #

function aiserver( ip = "127.0.0.1"
                 , port = "4242"
                 , dir = "."
                 , gpu = "false"
                 , index = "" )

  # Load models from model_dir
  models = load_models(dir, gpu)

  # Create a router that maps addresses to function calls
  router = HTTP.Router()

  # Register addresses
  HTTP.@register router "GET"  "/" (r->GET_index(r, index))

  HTTP.@register router "GET"  "/api/reload" (r->wrap(GET_reload, r, dir, gpu, models))
  HTTP.@register router "GET"  "/api/models" (r->wrap(GET_models, r, models))
  HTTP.@register router "GET"  "/api/games"  (r->wrap(GET_games, r))
  HTTP.@register router "GET"  "/api/create" (r->wrap(GET_create, r))

  HTTP.@register router "GET"  "/api/dummyimage" (r->wrap(GET_dummyimage, r))

  HTTP.@register router "POST" "/api/turn"   (r->wrap(POST_turn, r))
  HTTP.@register router "POST" "/api/apply" (r->wrap(POST_apply, r, models))

  # Start the server
  println("Server listens on $ip:$port.")

  try

    HTTP.serve(router, ip, port)

  catch err

    println(err)

  end

end

# -------- Parsing arguments and calling entry point ------------------------- #

s = ArgParseSettings()
@add_arg_table s begin
    "--ip"
    help = "IP address that the server listens to"
    default = "127.0.0.1"
  "--port"
    help = "port that the server listens to"
    arg_type = Int
    default = 4242
  "--dir"
    help = "directory path that is searched for '.jtm' model files"
    default = "."
  "--gpu"
    help = "whether the GPU will be used if CUDA is supported"
    action = :store_true
  "--index"
    help = "statically served index page"
    default = ""
end

args = parse_args(s)

aiserver(args["ip"], args["port"], args["dir"], args["gpu"], args["index"])

max_stars_repo_name: roSievers/jtac


text:
# Helper functions

# TODO: do this properly
# currently it doesn't check for types, except string and
# the code is horrible doing this zip between names (odd positions) and values (even positions)
# since options look like: "param1", value1, "param2", value2, ...
function setparams!(m::SCIPMathProgModel)
    a = m.options
    for (name, value) in zip(a[1][1:2:end],a[1][2:2:end])
        setparameter!(m, name, value)
    end
end

setparameter!(m::SCIPMathProgModel, name::ASCIIString, value::Bool) = _setBoolParam(m, pointer(name), Cint(value))
function setparameter!(m::SCIPMathProgModel, name::ASCIIString, value::Int)
    if _getParamType(m, pointer(name)) == 1 # integer
        return _setIntParam(m, pointer(name), Cint(value))
    else
        return _setLongintParam(m, pointer(name), Clonglong(value))
    end
end
setparameter!(m::SCIPMathProgModel, name::ASCIIString, value::Float64) = _setRealParam(m, pointer(name), Cdouble(value))
setparameter!(m::SCIPMathProgModel, name::ASCIIString, value::Char) = _setCharParam(m, pointer(name), Cchar(value))
setparameter!(m::SCIPMathProgModel, name::ASCIIString, value::ASCIIString) = _setStringParam(m, pointer(name), pointer(value))


text:
using Documenter

makedocs(;
    sitename="DiffusionMoments",
    pages = [
        "Index"=>"index.md"
    ]
)

max_stars_repo_name: SosUts/DiffusionMoments

max_issues_repo_name: SosUts/DiffusionMoments

max_forks_repo_name: SosUts/DiffusionMoments


text:
using LaTeXStrings, PyPlot, MosekTools, MarkovBounds
# use Catalyst to define the reaction network under consideration
# in this case, we investigate simple Michaelis-Menten kinetics
rn = @reaction_network begin
    10.0, S + E --> SE
    0.2, SE --> S + E
    4.0, SE --> P + E
    1.3, P --> S
end
S, E, SE, P = species(rn)
# specify initial molecular counts in system
x0 = Dict(S => 10, E => 10, SE => 0, P => 0)
x_scales = Dict(S => 5, E => 5, SE => 5, P => 5)
# compute bounds on mean molecular counts of the product counts at steady state
# for different truncation orders
orders = 2:6
sol = Dict()
for m in orders
    sol[m] = stationary_mean(rn, [S,P], x0, m, Mosek.Optimizer, auto_scaling = false, scales = x_scales)
end

for s in [S, P]
    fig, ax = subplots()
    ax.plot(orders, [sol[m].bounds[s][1] for m in orders], marker = "*", color="blue", label="lower bound")
    ax.plot(orders, [sol[m].bounds[s][2] for m in orders], marker = "*", color="red", label="upper bound")
    ax.set_xlabel("truncation order")
    ax.set_ylabel(string(L"⟨",s.f.name,L"_{\infty}⟩"))
    ax.legend()
    display(fig)
end

# same for variance
sol = Dict()
for m in orders
    sol[m] = stationary_variance(rn, [S,P], x0, m, Mosek.Optimizer, auto_scaling = false, scales = x_scales)
end

for s in [S, P]
    fig, ax = subplots()
    ax.plot(orders, [sol[m].bounds[s] for m in orders], marker = "*", color="red", label="upper bound")
    ax.set_xlabel("truncation order")
    ax.set_ylabel(string(L"Var(",s.f.name,L"_{\infty})"))
    ax.legend()
    display(fig)
end

# transient moments
nT = 5
traj = []
Tf_range = [0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3,
            0.4, 0.5, 0.625, 0.75, 1.0, 1.25, 1.5, 2.0, 3.0, 4.0, 5.0]
order = 4
for Tf in Tf_range
    trange = range(Tf/nT, stop=Tf, length=nT)
    sol = transient_mean(rn, [S,P], x0, order, trange, Mosek.Optimizer, scales = x_scales)
    push!(traj, sol)
end

for s in [S,P]
    fig, ax = subplots()
    ax.plot(Tf_range, [sol.bounds[s][1] for sol in traj], color="blue", label = "lower bound")
    ax.plot(Tf_range, [sol.bounds[s][2] for sol in traj], color="red", label = "upper bound")
    ax.set_xlabel("time")
    ax.set_ylabel(string(L"⟨", s.f.name,L"⟩"))
    ax.legend()
    display(fig)
end

# variance
nT = 5
traj = []
Tf_range = [0.025, 0.05, 0.075, 0.1, 0.15, 0.2, 0.25, 0.3,
            0.4, 0.5, 0.625, 0.75, 1.0, 1.25, 1.5, 2.0, 3.0, 4.0, 5.0]
order = 4
for Tf in Tf_range
    trange = range(Tf/nT, stop=Tf, length=nT)
    sol = transient_variance(rn, [P], x0, order, trange, Mosek.Optimizer, scales =  x_scales)
    push!(traj, sol)
end

ps = []
for s in [P]
    fig, ax = subplots()
    ax.plot(Tf_range, [sol.bounds[s][1] for sol in traj], color="red", label = "upper bound")
    ax.set_xlabel("time")
    ax.set_ylabel(string(L"Var(", s.f.name,L")"))
    ax.legend()
    display(fig)
end

## volume of confidence ellipsoid
orders = 2:8
sol = Dict()
for m in orders
    sol[m] = stationary_confidence_ellipsoid(rn, x0, [S,P], m, Mosek.Optimizer; scales = x_scales)
end

fig, ax = subplots()
ax.plot(orders, [sol[m].bounds[[S,P]] for m in orders], marker="*", color="red", label = "upper bound")
ax.set_xlabel("order")
ax.set_ylabel("Volume of S-P confidence ellipsoids")
ax.legend()
display(fig)


text:
export ft_setresetpiperetrycount

function ft_setresetpiperetrycount(ft_handle::Culong, count::Integer)
  ft_status = ccall((:FT_SetResetPipeRetryCount, d2xx),
                     Cuint,
                     (Culong, Cuint),
                     ft_handle, count)
  checkstatus(ft_status)
  return nothing
end


text:
module Grids

using ..OceanTurbulenceParameterEstimation

using Oceananigans
using Oceananigans.Grids: Flat, Bounded, Periodic,
                          Face, Center,
                          AbstractGrid,
                          RegularRectilinearGrid

using JLD2

export OneDimensionalEnsembleGrid

include("utils.jl")
include("one_dimensional_ensemble_grid.jl")

end #module


text:
using BenchmarkTools, AMDGPU

DAT = Float32

if DAT==Float64
    sc = 1
    DAT_Int = Int32
elseif DAT==Float32
    sc = 2
    DAT_Int = Int32
end

@inbounds function memcopy_triad_pow_int!(A, B, C, s, pow_int, ni)
    ix = (workgroupIdx().x - 1) * workgroupDim().x + workitemIdx().x
    iy = (workgroupIdx().y - 1) * workgroupDim().y + workitemIdx().y
    for ii=1:ni
        A[ix,iy] = B[ix,iy] + s*C[ix,iy]^pow_int
    end
    return
end

@inbounds function memcopy_triad_pow_float!(A, B, C, s, pow_float, ni)
    ix = (workgroupIdx().x - 1) * workgroupDim().x + workitemIdx().x
    iy = (workgroupIdx().y - 1) * workgroupDim().y + workitemIdx().y
    for ii=1:ni
        A[ix,iy] = B[ix,iy] + s*C[ix,iy]^pow_float
    end
    return
end

function run_bench()
    fact      = 24
    ni        = 100
    nx, ny    = sc*fact*1024, fact*1024
    threads   = (256, 2)
    grid      = (nx, ny)
    A         = ROCArray(zeros(DAT, nx, ny))
    B         = ROCArray( rand(DAT, nx, ny))
    C         = ROCArray( ones(DAT, nx, ny))
    pow_int   = DAT_Int(3)
    pow_float = DAT(3.75)
    s         = rand(DAT)
    lam       = rand(DAT)
    _dx, _dy  = DAT(1.0), DAT(1.0)
    dt        = DAT(1.0/10.0/4.1)
    println("nx, ny = $(nx), $(ny)")
    # run test 2
    println("test $(DAT)^$(DAT_Int)")
    t_it = @belapsed begin wait( @roc groupsize=$threads gridsize=$grid memcopy_triad_pow_int!($A, $B, $C, $s, $pow_int, $ni) ) end
    T_tot = 3*1/1e9*nx*ny*sizeof(DAT)/t_it
    Flops = 3*1/1e9*nx*ny*ni/t_it
    println("T_tot triad2D pow_int   = $(round(T_tot,sigdigits=7)) GB/s")
    println("Flops triad2D pow_int   = $(round(Flops,sigdigits=7)) GFLOP/s")
    # run test 3
    println("test $(DAT)^$(DAT)")
    t_it = @belapsed begin wait( @roc groupsize=$threads gridsize=$grid memcopy_triad_pow_float!($A, $B, $C, $s, $pow_float, $ni) ) end
    T_tot = 3*1/1e9*nx*ny*sizeof(DAT)/t_it
    Flops = 3*1/1e9*nx*ny*ni/t_it
    println("T_tot triad2D pow_float = $(round(T_tot,sigdigits=7)) GB/s")
    println("Flops triad2D pow_float = $(round(Flops,sigdigits=7)) GFLOP/s")
    return
end

run_bench()

max_stars_repo_name: luraess/JuliaGPUPerf

max_issues_repo_name: luraess/JuliaGPUPerf

max_forks_repo_name: luraess/JuliaGPUPerf


text:
using Test, Random, AutoGrad, Statistics
using Knet.Ops21: conv, relu
using Knet.Layers21: BatchNorm
using Knet.KnetArrays: KnetArray
using CUDA: CUDA, CuArray

@testset "ops21/conv" begin

    function convtest(
        ; atype = Array,
        nd = 4,
        activation = nothing,
        normalization = nothing,
        alpha = 1,
        beta = 0,
        channelmajor = false,
        crosscorrelation = false,
        dilation = 1,
        groups = 1,
        padding = 0,
        stride = 1,
        usebias = false,
        usez = false,
        gcheck_test = true,
        approx_test = true,
    )
        cx, cy = 3*groups, 2*groups
        xd = rand(5:10, nd); xd[nd-1] = cx
        wd = ((3*ones(Int,nd-2))..., cx÷groups, cy)
        x = randn(xd...)
        w = randn(wd...)
        b = usebias ? randn(ones(Int,nd-2)...,cy,1) : nothing
        cm(x) = (x===nothing ? nothing : permutedims(x, (nd-1,(1:nd-2)...,nd)))
        if channelmajor; (x,w,b) = cm.((x,w,b)); end
        kw = (; activation, alpha, beta, channelmajor, crosscorrelation, dilation, groups, normalization, padding, stride)
        z = usez ? randn!(conv(w, x; kw...)) : nothing
        ax,aw,ab,az = (i->i===nothing ? nothing : atype(i)).((x,w,b,z))
        px,pw,pb,pz = (i->i===nothing ? nothing : Param(i)).((ax,aw,ab,az))
        r1 = ((isa(aw,Array) || !approx_test) ? true :
              isapprox(conv(w, x; bias=b, z=z, kw...),
                       Array(conv(aw, ax; bias=ab, z=az, kw...))))
        atol,rtol = (eltype(pw) == Float64 ? (0.01, 0.05) : (0.1, 0.5))
        r2 = !gcheck_test ? true : @gcheck conv(pw, px; bias=pb, z=pz, kw...) (;rtol,atol)
        r1 && r2
    end

    @test convtest(; )
    @test convtest(; nd=3)
    @test convtest(; nd=5)
    @test convtest(; activation=relu)
    @test convtest(; normalization=(x->2x.+1))
    @test convtest(; alpha=2)
    @test convtest(; beta=2, usez=true)
    @test convtest(; channelmajor=true)
    @test convtest(; crosscorrelation=true)
    @test convtest(; dilation=2)
    @test convtest(; groups=2)
    @test convtest(; padding=1)
    @test convtest(; stride=2)
    @test convtest(; usebias=true)

    if CUDA.functional()

        @test convtest(; atype=CuArray, )
        @test convtest(; atype=CuArray, nd=3)
        @test convtest(; atype=CuArray, nd=5)
        @test convtest(; atype=CuArray, activation=relu)
        @test convtest(; atype=CuArray, normalization=(x->2x.+1))
        @test convtest(; atype=CuArray, alpha=2)
        @test convtest(; atype=CuArray, beta=2, usez=true)
        # cudnn8 does not support Float64 with channelmajor
        @test convtest(; atype=CuArray{Float32}, channelmajor=true, gcheck_test=false)
        @test convtest(; atype=CuArray, crosscorrelation=true)
        @test convtest(; atype=CuArray, dilation=2)
        @test convtest(; atype=CuArray, groups=2)
        @test convtest(; atype=CuArray, padding=1)
        @test convtest(; atype=CuArray, stride=2)
        @test convtest(; atype=CuArray, usebias=true)

        @test convtest(; atype=KnetArray, )
        @test convtest(; atype=KnetArray, nd=3)
        @test convtest(; atype=KnetArray, nd=5)
        @test convtest(; atype=KnetArray, activation=relu)
        @test convtest(; atype=KnetArray, normalization=(x->2x.+1))
        @test convtest(; atype=KnetArray, alpha=2)
        @test convtest(; atype=KnetArray, beta=2, usez=true)
        @test convtest(; atype=KnetArray{Float32}, channelmajor=true, gcheck_test=false)
        @test convtest(; atype=KnetArray, crosscorrelation=true)
        @test convtest(; atype=KnetArray, dilation=2)
        @test convtest(; atype=KnetArray, groups=2)
        @test convtest(; atype=KnetArray, padding=1)
        @test convtest(; atype=KnetArray, stride=2)
        @test convtest(; atype=KnetArray, usebias=true)

        # Multiple use of the same conv weight
        w = Param(CUDA.randn(Float64, 3,3,3,3))
        x = Param(CUDA.randn(Float64, 8,8,3,1))
        @test @gcheck conv(w,conv(w,x))
    end
end


text:
# This file is a part of SimilaritySearch.jl

"""
    execute_callback(opt::NeighborhoodSize, index)

SearchGraph's callback for adjusting neighborhood strategy
"""
function execute_callback(opt::NeighborhoodSize, index)
    N = index.neighborhood
    N.ksearch = ceil(Int, N.minsize + log(N.logbase, length(index)))
end

"""
    find_neighborhood(index::SearchGraph{T}, item, pools; hints=index.hints)

Searches for `item` neighborhood in the index, i.e., if `item` were in the index whose items should be its neighbors (intenal function).
`res` is always reused since `reduce` creates a new KnnResult from it (a copy if `reduce` in its simpler terms)
"""
function find_neighborhood(index::SearchGraph, item, pools::SearchGraphPools; hints=index.hints)
    n = length(index)
    if n > 0
        res = getknnresult(index.neighborhood.ksearch, pools)
        search(index.search_algo, index, item, res, hints, pools)
        reduce_neighborhood(index.neighborhood.reduce, index, item, res, pools)
    else
        Int32[]
    end
end

"""
    push_neighborhood!(index::SearchGraph, item, neighbors, callbacks; push_item=true)

Inserts the object `item` into the index, i.e., creates an edge from items listed in L and the
vertex created for ìtem` (internal function)
"""
function push_neighborhood!(index::SearchGraph, item, neighbors, callbacks; push_item=true)
    push_item && push!(index.db, item)
    push!(index.links, neighbors)
    push!(index.locks, Threads.SpinLock())
    n = length(index)
    n == 1 && return
    ## vstate = getvisitedvertices(index)
    @inbounds for id in neighbors
        push!(index.links[id], n)  # sat push?
    end

    callbacks !== nothing && execute_callbacks(callbacks, index)

    if index.verbose && length(index) % 100_000 == 0
        println(stderr, "added n=$(length(index)), neighborhood=$(length(neighbors)), $(string(index.search_algo)), $(Dates.now())")
    end
end

const GlobalSatKnnResult = [KnnResult(1)]
function __init__neighborhood()
    for _ in 2:Threads.nthreads()
        push!(GlobalSatKnnResult, KnnResult(1))
    end
end

"""
    SatNeighborhood()

New items are connected with a small set of items computed with a SAT like scheme (**cite**).
It starts with `k` near items that are reduced to a small neighborhood due to the SAT partitioning stage.
"""
struct SatNeighborhood <: NeighborhoodReduction end
Base.copy(::SatNeighborhood) = SatNeighborhood()

"""
    DistalSatNeighborhood()

New items are connected with a small set of items computed with a Distal SAT like scheme (**cite**).
It starts with `k` near items that are reduced to a small neighborhood due to the SAT partitioning stage but in reverse order of distance.
"""
struct DistalSatNeighborhood <: NeighborhoodReduction end
Base.copy(::DistalSatNeighborhood) = DistalSatNeighborhood()

"""
    struct IdentityNeighborhood

It does not modifies the given neighborhood
"""
struct IdentityNeighborhood <: NeighborhoodReduction end
Base.copy(::IdentityNeighborhood) = IdentityNeighborhood()

## functions

function sat_should_push(sat_neighborhood::T, index, item, id, dist, near::KnnResult) where T
    @inbounds obj = index[id]
    dist = dist < 0.0 ? evaluate(index.dist, item, obj) : dist
    push!(near, 0, dist)

    @inbounds for linkID in sat_neighborhood
        d = evaluate(index.dist, index[linkID], obj)
        push!(near, linkID, d)
    end

    argmin(near) == 0
end

function reduce_neighborhood(::IdentityNeighborhood, index::SearchGraph, item, res, pools::SearchGraphPools)
    copy(res.id)
end

"""
    reduce(sat::DistalSatNeighborhood, index::SearchGraph, item, res, pools)

Reduces `res` using the DistSAT strategy.
"""
@inline function reduce_neighborhood(::DistalSatNeighborhood, index::SearchGraph, item, res, pools::SearchGraphPools)
    N = Vector{Int32}(undef, 2)
    resize!(N, 0)
    @inbounds for i in length(res):-1:1  # DistSat => works a little better but also produces a bit larger neighborhoods
        id, dist = getpair(res, i)
        sat_should_push(N, index, item, id, dist, getsatknnresult(pools)) && push!(N, id)
    end

    N
end

@inline function reduce_neighborhood(::SatNeighborhood, index::SearchGraph, item, res, pools::SearchGraphPools)
    N = Vector{Int32}(undef, 2)
    resize!(N, 0)
    @inbounds for (id, dist) in res  # DistSat => works a little better but also produces a bit larger neighborhoods
        sat_should_push(N, index, item, id, dist, getsatknnresult(pools)) && push!(N, id)
    end

    N
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
module Checkerboard

using LinearAlgebra
using SparseArrays

export checkerboard_groups!, checkerboard_groups, checkerboard_order!, checkerboard_order
export checkerboard_mul!, checkerboard_transpose_mul!, checkerboard_inverse_mul!, checkerboard_inverse_transpose_mul!
export checkerboard_matrix

"""
Construct full checkerboard matrix. For code
testing rather than use in the final full simulation.
"""
function checkerboard_matrix(neighbor_table::Matrix{Int},coshtij::Vector{T1}, sinhtij::Vector{T1};transposed::Bool=false,threshold::T2=1e-12) where {T1<:Number,T2<:AbstractFloat}

    # to contain M[row,col]=val info for constructing sparse matrix
    rows = Int[]
    cols = Int[]
    elements = T1[]
    # size of matrix
    N = maximum(neighbor_table)
    # stores columns vector
    colvector = zeros(T1,N)
    # iterating over rows
    for col in 1:N
        # initialize column vector as unit vector
        colvector      .= 0.0
        colvector[col]  = 1.0
        # doing matrix vector multiply
        if transposed==true
            # multiply unit vector by Mᵀ matrix
            checkerboard_transpose_mul!(colvector,neighbor_table,coshtij,sinhtij)
        else
            # multiply unit vector by M matrix
            checkerboard_mul!(colvector,neighbor_table,coshtij,sinhtij)
        end
        # iterate of column vecto
        for row in 1:N
            # if nonzero
            if abs(colvector[row])>threshold
                # save matrix element
                append!(rows,row)
                append!(cols,col)
                append!(elements,colvector[row])
            end
        end
    end
    return rows, cols, elements
end

"""
In-place multiplication of vector with checkerboard matrix.
This method assumes the `neighbor_table` and associated `coshs` and `sinhs` are already ordered correctly
for the checkerboard decomposition.
"""
function checkerboard_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2}, Lτ::Int) where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in 1:length(coshtij)

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 + s*t2
            y[idx_j] = c*t2 + conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Matrix{T2}, sinhtij::Matrix{T2}) where {T1<:Number,T2<:Number}

    # number of bonds
    N = size(coshtij,2)

    # length of imaginary time axis
    Lτ = size(coshtij,1)

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in 1:N

        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            # get matrix elements
            c = coshtij[τ,n]
            s = sinhtij[τ,n]

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            # perform matrix-vector multiply
            y[idx_i] = c*t1 + s*t2
            y[idx_j] = c*t2 + conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2}) where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in 1:length(coshtij)

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        t1 = y[i]
        t2 = y[j]

        y[i] = c*t1 + s*t2
        y[j] = c*t2 + conj(s)*t1
    end

    return nothing
end

"""
In-place multiplication of vector with checkerboard matrix.
This method assumes the `neighbor_table` and associated `coshs` and `sinhs` are already ordered correctly
for the checkerboard decomposition.
"""
function checkerboard_transpose_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2}, Lτ::Int)  where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in length(coshtij):-1:1

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 + s*t2
            y[idx_j] = c*t2 + conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_transpose_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Matrix{T2}, sinhtij::Matrix{T2})  where {T1<:Number,T2<:Number}

    # number of bonds
    N = size(coshtij,2)

    # length of imaginary time axis
    Lτ = size(coshtij,1)

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in N:-1:1

        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            c = coshtij[τ,n]
            s = sinhtij[τ,n]

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 + s*t2
            y[idx_j] = c*t2 + conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_transpose_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2})  where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in length(coshtij):-1:1

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        t1 = y[i]
        t2 = y[j]

        y[i] = c*t1 + s*t2
        y[j] = c*t2 + conj(s)*t1
    end

    return nothing
end

"""
In-place multiplication of vector with the inverse of the checkerboard matrix.
This method assumes the `neighbor_table` and associated `coshs` and `sinhs` are already ordered correctly
for the checkerboard decomposition.
"""
function checkerboard_inverse_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2}, Lτ::Int)  where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in length(coshtij):-1:1

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 - s*t2
            y[idx_j] = c*t2 - conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_inverse_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Matrix{T2}, sinhtij::Matrix{T2})  where {T1<:Number,T2<:Number}

    N  = size(coshtij,2)
    Lτ = size(coshtij,1)

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in N:-1:1

        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            c = coshtij[τ,n]
            s = sinhtij[τ,n]

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 - s*t2
            y[idx_j] = c*t2 - conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_inverse_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2})  where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in length(coshtij):-1:1

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        t1 = y[i]
        t2 = y[j]

        y[i] = c*t1 - s*t2
        y[j] = c*t2 - conj(s)*t1
    end

    return nothing
end

"""
In-place multiplication of vector with the inverse of the checkerboard matrix.
This method assumes the `neighbor_table` and associated `coshs` and `sinhs` are already ordered correctly
for the checkerboard decomposition.
"""
function checkerboard_inverse_transpose_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2}, Lτ::Int)  where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in 1:length(coshtij)

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 - s*t2
            y[idx_j] = c*t2 - conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_inverse_transpose_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Matrix{T2}, sinhtij::Matrix{T2})  where {T1<:Number,T2<:Number}

    N  = size(coshtij,2)
    Lτ = size(coshtij,1)

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in 1:N

        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        # iterating over time slices
        @simd for τ in 1:Lτ

            c = coshtij[τ,n]
            s = sinhtij[τ,n]

            # mapping (τ,site) ==> index
            idx_i = (i-1)*Lτ + τ
            idx_j = (j-1)*Lτ + τ

            t1 = y[idx_i]
            t2 = y[idx_j]

            y[idx_i] = c*t1 - s*t2
            y[idx_j] = c*t2 - conj(s)*t1
        end
    end

    return nothing
end

function checkerboard_inverse_transpose_mul!(y::AbstractVector{T1}, neighbor_table_tij::Matrix{Int}, coshtij::Vector{T2}, sinhtij::Vector{T2})  where {T1<:Number,T2<:Number}

    # iterating over pairs of neighboring sites
    @fastmath @inbounds for n in 1:length(coshtij)

        c = coshtij[n]
        s = sinhtij[n]
        i = neighbor_table_tij[1, n]
        j = neighbor_table_tij[2, n]

        t1 = y[i]
        t2 = y[j]

        y[i] = c*t1 - s*t2
        y[j] = c*t2 - conj(s)*t1
    end

    return nothing
end

"""
    function checkerboard_order(neighbor_table::Matrix{Int})::Vector{Int}

Given a `neighbor_table`, this functions determines the `order` the checkerboard ordering
for the neighboring sites. Assumes the `neighbor_table` has already been sorted.
"""
function checkerboard_order(neighbor_table::Matrix{Int})::Vector{Int}

    nneighbor_table = size(neighbor_table,2)
    groups = zeros(Int,nneighbor_table)
    order = zeros(Int,nneighbor_table)
    checkerboard_groups!(groups,neighbor_table)
    checkerboard_order!(order,groups)
    return order
end

"""
    function checkerboard_order(groups::Vector{Int})::Vector{Int}

Given a checkerboard groupings of neighbor_table `groups`, this functions determines
the `order` the associated neighbor_table should be iterated over for performing Matrix
multiplications.
"""
function checkerboard_order(groups::Vector{Int})::Vector{Int}

    order = zeros(Int,length(groups))
    checkerboard_order!(order,groups)
    return order
end

"""
    function checkerboard_order!(order::Vector{Int},groups::Vector{Int})

Given a checkerboard groupings of neighbor_table `groups`, this functions determines
the `order` the associated neighbor_table should be iterated over for performing Matrix
multiplications.
"""
function checkerboard_order!(order::Vector{Int},groups::Vector{Int})

    sortperm!(order,groups)
    return nothing
end

"""
    checkerboard_groups(neighbor_table::Matrix{Int})::Vector{Int}

Constructs the checkerboard decomposition for a given neighbor table.
Assumes the neighbor table has already been sorted using the sort_neighbor_table!
method from the Lattices module.
Returns the number of groups in the checkerboard decomposition.
"""
function checkerboard_groups(neighbor_table::Matrix{Int})::Vector{Int}

    groups = zeros(Int,size(neighbor_table,2))
    checkerboard_groups!(groups,neighbor_table)
    return groups
end

"""
    checkerboard_groups!(groups::Vector{Int},neighbor_table::Matrix{Int})::Int

Constructs the checkerboard decomposition for a given neighbor table.
Assumes the neighbor table has already been sorted using the sort_neighbor_table! method from the Lattices module.
Returns the number of groups in the checkerboard decomposition.
"""
function checkerboard_groups!(groups::Vector{Int},neighbor_table::Matrix{Int})::Int

    # checking dimensions
    @assert size(neighbor_table,2)==length(groups)
    @assert size(neighbor_table,1)==2
    # getting the total number of neighbor pairs
    nneighbors = size(neighbor_table,2)
    # intially all neighbors are unassigned to a group
    fill!(groups,0)
    # keeps track of which group is being constructed
    group = 0
    # keep track of number of bonds assigned to group
    nassigned = 0
    # while any bond is not assigned to a group
    while nassigned < nneighbors
        # increment to next group
        group += 1
        # iterate over neighbor_table in lattice
        for neighbor in 1:nneighbors
            # if neighbor is not assigned to a group
            if groups[neighbor]==0
                # assign it to current group
                groups[neighbor] = group
                nassigned       += 1
                # iterate over previous neighbors
                for prev_neighbor in 1:neighbor-1
                    # if previous neighbor is a group member
                    if groups[prev_neighbor]==group
                        # if the previous neighbor overlaps with current neighbor
                        if ( neighbor_table[1,neighbor]==neighbor_table[1,prev_neighbor] ||
                             neighbor_table[2,neighbor]==neighbor_table[2,prev_neighbor] ||
                             neighbor_table[1,neighbor]==neighbor_table[2,prev_neighbor] ||
                             neighbor_table[2,neighbor]==neighbor_table[1,prev_neighbor] )
                            # remove current neighbor from group
                            groups[neighbor] = 0
                            nassigned       -= 1
                            break
                        end
                    end
                end
            end
        end
    end
    return group
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
include("slab_gen.jl")
using slab: r2cSlab, c2cSlab

@assert isdefined(:r2cSlab)
@assert isdefined(:c2cSlab)


text:
"""
    serialize_model(model::MM, filename::String)::Serialized{MM} where {MM<:MetabolicModel}

Serialize the `model` to file `filename`, returning a [`Serialized`](@ref)
model that is able to load itself back automatically upon precaching by
[`precache!`](@ref).
"""
function serialize_model(
    model::MM,
    filename::String,
)::Serialized{MM} where {MM<:MetabolicModel}
    open(f -> serialize(f, model), filename, "w")
    Serialized{MM}(nothing, filename)
end

"""
    serialize_model(model::Serialized, filename::String)::Serialized

Specialization of [`serialize_model`](@ref) that prevents nested serialization
of already-serialized models.
"""
function serialize_model(model::Serialized, filename::String)
    precache!(model)
    serialize_model(model.m, filename)
end


text:
"""
    dx!(dx,x,p::Parameters,t)

Function defining change in amount of consumers and resources over time.

Dynamics occur as described by the set of equations:

``\\frac{dC_i}{dt} = C_i \\left(∑\\right)``

# """
function dx!(dx,x,p::Parameters,t)
  #loop over consumers i = 1:N
  for i = 1:p.N
    dx[p.M + i] = -x[p.M + i] * p.Rm[i]
    #loop over metabolites
    for j = 1:p.M
      dx[p.M + i] +=  x[p.M + i] * x[j] * p.u[i,j] * (1 - p.l_sum[j]) # u_ij * (1-l_j)
    end
    if (x[p.M + i] + dx[p.M + i]) < eps()
      dx[p.M + i] = -x[p.M + i]
    end
  end

 #loop over metabolites j = 1:M
  for j = 1:p.M
    dx[j] = p.ρ[j] - x[j] * p.ω[j]#inflow
    #loop over species
    for i = 1:p.N
      dx[j] += -p.u[i,j] * x[p.M+ i] * x[j] #minus uptake
      #loop over other metabolites
      for k = 1:p.M
        dx[j] += p.u[i,k] * x[p.M + i] * x[k] * p.l[k,j] # +leakage
      end
    end
  end

  # dx[dx .< eps()] .= -x[dx .< eps()]

end

function dx_threaded!(dx,x,p::Parameters,t)
  #loop over consumers i = 1:N
  Threads.@threads for i = 1:p.N
    dx[p.M + i] = -x[p.M + i] * p.Rm[i]
    #loop over metabolites
    for j = 1:p.M
      dx[p.M + i] +=  x[p.M + i] * x[j] * p.u[i,j] * (1 - p.l_sum[j]) # u_ij * (1-l_j)
    end
    if (x[p.M + i] + dx[p.M + i]) < eps()
      dx[p.M + i] = -x[p.M + i]
    end
  end

 #loop over metabolites j = 1:M
 Threads.@threads for j = 1:p.M
    dx[j] = p.ρ[j] - x[j] * p.ω[j]#inflow
    #loop over species
    for i = 1:p.N
      dx[j] += -p.u[i,j] * x[p.M+ i] * x[j] #minus uptake
      #loop over other metabolites
      for k = 1:p.M
        dx[j] += p.u[i,k] * x[p.M + i] * x[k] * p.l[k,j] # +leakage
      end
    end
  end

  # dx[dx .< eps()] .= -x[dx .< eps()]

end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

module WorldModel

export
    AbstractWorldModel,

    initialize,
    update,

    updateObjectState,
    getAllObjectStates,
    getAll,
    updateAll,

    AirSpace,
    ASWMState

import CommonInterfaces: initialize, update

using AbstractWorldModelImpl

using AirSpaceImpl

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
type ScalarOrderStatistic{D,S<:MultivariateDistribution} <: AbstractScalarOrderStatistic{D}
    sequence::S
    order::Int64
end

parentsequence(X::ScalarOrderStatistic) = X.sequence
order(X::ScalarOrderStatistic) = X.order

function rand(X::ScalarOrderStatistic)
    samples = sort(rand(parentsequence(X)))
    return samples[X.order]
end

function mean(X::ScalarOrderStatistic)
    f = x -> pdf(X,x)*x
    return quadgk(f,minimum(X),maximum(X))[1]#this seems to be rather ill conditioned
end

#
# IID SEQUENCES
#
function ScalarOrderStatistic{D,S}(sequence::IIDRandomSequence{D,S},order::Int64)
    return ScalarOrderStatistic{D,IIDRandomSequence{D,S}}(sequence,order)
end

typealias IIDScalarOrderStatistic{D<:ValueSupport,T<:UnivariateDistribution} ScalarOrderStatistic{D,IIDRandomSequence{D,T}}
typealias IIDContinuousScalarOrderStatistic{T<:ContinuousUnivariateDistribution} IIDScalarOrderStatistic{Continuous,T}
typealias IIDDiscreteScalarOrderStatistic{T<:DiscreteUnivariateDistribution} IIDScalarOrderStatistic{Discrete,T}

# (having probelsm with Discrete/Continuous DRY due to ambiguity creation)
minimum{T<:ContinuousUnivariateDistribution}(X::IIDContinuousScalarOrderStatistic{T}) = minimum(X.sequence.d)
minimum{T<:DiscreteUnivariateDistribution}(X::IIDDiscreteScalarOrderStatistic{T}) = minimum(X.sequence.d)
maximum{T<:ContinuousUnivariateDistribution}(X::IIDContinuousScalarOrderStatistic{T}) = maximum(X.sequence.d)
maximum{T<:DiscreteUnivariateDistribution}(X::IIDDiscreteScalarOrderStatistic{T}) = maximum(X.sequence.d)

function cdf{T<:ContinuousUnivariateDistribution}(X::IIDContinuousScalarOrderStatistic{T}, x::Real)
  P,n,k = cdf(X.sequence.d,x), length(X.sequence), X.order
  return ccdf(Binomial(n,P),k-1)
end
function cdf{T<:DiscreteUnivariateDistribution}(X::IIDDiscreteScalarOrderStatistic{T}, x::Int64)
  P,n,k = cdf(X.sequence.d,x), length(X.sequence), X.order
  return ccdf(Binomial(n,P),k-1)
end

function pdf{T<:ContinuousUnivariateDistribution}(X::ScalarOrderStatistic{Continuous,IIDRandomSequence{Continuous,T}}, x::Real)
  P,p,n,k = cdf(X.sequence.d,x), pdf(X.sequence.d,x),length(X.sequence),X.order
  coeff = Γ(n+1)/(Γ(k)*Γ(n-k+1))
  return coeff*P^(k-1)*(1.-P)^(n-k)*p
end

#
# INID sequence
#
function ScalarOrderStatistic{D,S}(sequence::INIDRandomSequence{D,S},order::Int64)
    return ScalarOrderStatistic{D,INIDRandomSequence{D,S}}(sequence,order)
end


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

@doc raw"""Represents a Persistent Disk resource in Google Compute Engine.  A GCE PD must exist before mounting to a container. The disk must also be in the same GCE project and zone as the kubelet. A GCE PD can only be mounted as read/write once or read-only many times. GCE PDs support ownership management and SELinux relabeling.

    IoK8sApiCoreV1GCEPersistentDiskVolumeSource(;
        fsType=nothing,
        partition=nothing,
        pdName=nothing,
        readOnly=nothing,
    )

    - fsType::String : Filesystem type of the volume that you want to mount. Tip: Ensure that the filesystem type is supported by the host operating system. Examples: \"ext4\", \"xfs\", \"ntfs\". Implicitly inferred to be \"ext4\" if unspecified. More info: https://kubernetes.io/docs/concepts/storage/volumes#gcepersistentdisk
    - partition::Int32 : The partition in the volume that you want to mount. If omitted, the default is to mount by volume name. Examples: For volume /dev/sda1, you specify the partition as \"1\". Similarly, the volume partition for /dev/sda is \"0\" (or you can leave the property empty). More info: https://kubernetes.io/docs/concepts/storage/volumes#gcepersistentdisk
    - pdName::String : Unique name of the PD resource in GCE. Used to identify the disk in GCE. More info: https://kubernetes.io/docs/concepts/storage/volumes#gcepersistentdisk
    - readOnly::Bool : ReadOnly here will force the ReadOnly setting in VolumeMounts. Defaults to false. More info: https://kubernetes.io/docs/concepts/storage/volumes#gcepersistentdisk
"""
mutable struct IoK8sApiCoreV1GCEPersistentDiskVolumeSource <: SwaggerModel
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    partition::Any # spec type: Union{ Nothing, Int32 } # spec name: partition
    pdName::Any # spec type: Union{ Nothing, String } # spec name: pdName
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly

    function IoK8sApiCoreV1GCEPersistentDiskVolumeSource(;fsType=nothing, partition=nothing, pdName=nothing, readOnly=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1GCEPersistentDiskVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1GCEPersistentDiskVolumeSource, Symbol("partition"), partition)
        setfield!(o, Symbol("partition"), partition)
        validate_property(IoK8sApiCoreV1GCEPersistentDiskVolumeSource, Symbol("pdName"), pdName)
        setfield!(o, Symbol("pdName"), pdName)
        validate_property(IoK8sApiCoreV1GCEPersistentDiskVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        o
    end
end # type IoK8sApiCoreV1GCEPersistentDiskVolumeSource

const _property_map_IoK8sApiCoreV1GCEPersistentDiskVolumeSource = Dict{Symbol,Symbol}(Symbol("fsType")=>Symbol("fsType"), Symbol("partition")=>Symbol("partition"), Symbol("pdName")=>Symbol("pdName"), Symbol("readOnly")=>Symbol("readOnly"))
const _property_types_IoK8sApiCoreV1GCEPersistentDiskVolumeSource = Dict{Symbol,String}(Symbol("fsType")=>"String", Symbol("partition")=>"Int32", Symbol("pdName")=>"String", Symbol("readOnly")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1GCEPersistentDiskVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1GCEPersistentDiskVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1GCEPersistentDiskVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1GCEPersistentDiskVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1GCEPersistentDiskVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1GCEPersistentDiskVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1GCEPersistentDiskVolumeSource)
    (getproperty(o, Symbol("pdName")) === nothing) && (return false)
    true
end

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


text:
using SparseArrays # spzeros
function ybus( buses, lines )

nbus = length(buses)
nline = length(lines)
busset = 1:nbus
lineset = 1:nline

Ybus = spzeros(ComplexF64, nbus, nbus)
yff = spzeros(ComplexF64, nline, 1)
yft = spzeros(ComplexF64, nline, 1)
ytf = spzeros(ComplexF64, nline, 1)
ytt = spzeros(ComplexF64, nline, 1)
tau = spzeros(ComplexF64, nline, 1)
for l in lineset
    tau[l] = ifelse( lines[l].tap == 0, exp( + im * lines[l].shft), lines[l].tap * exp( + im * lines[l].shft) )

    ytt[l] = inv(lines[l].r + im * lines[l].x) + im * lines[l].b / 2
    yff[l] = (inv(lines[l].r + im * lines[l].x) + im * lines[l].b / 2 ) / (tau[l] .* conj(tau[l]))
    yft[l] = - (inv(lines[l].r + im * lines[l].x)) / conj(tau[l])
    ytf[l] = - (inv(lines[l].r + im * lines[l].x)) / tau[l]

    Ybus[lines[l].fbus, lines[l].tbus] = yft[l]
    Ybus[lines[l].tbus, lines[l].fbus] = ytf[l]
    Ybus[lines[l].fbus, lines[l].fbus] += yff[l]
    Ybus[lines[l].tbus, lines[l].tbus] += ytt[l]
end
for b in busset
    Ybus[b,b] += buses[b].Gs + im * buses[b].Bs
end

return Ybus, yff, yft, ytf, ytt
end

max_stars_repo_name: JipKim/ACOPF

max_issues_repo_name: JipKim/ACOPF

max_forks_repo_name: JipKim/ACOPF