text:
import OceanTurb

function free_convection_tke_mass_flux(ds; parameters=OceanTurb.TKEMassFlux.TKEParameters(), Δt=600)

    ρ₀ = 1027.0
    cₚ = 4000.0
    f  = ds.metadata["coriolis_parameter"]
    α  = ds.metadata["thermal_expansion_coefficient"]
    β  = 0.0
    g  = ds.metadata["gravitational_acceleration"]
    constants = OceanTurb.Constants(Float64, ρ₀=ρ₀, cP=cₚ, f=f, α=α, β=β, g=g)

    T = ds["T"]
    wT = ds["wT"]
    N = Nz = size(T, 3)
    zc = znodes(T)
    zf = znodes(wT)
    H = abs(zf[1])
    Nt = size(T, 4)
    times = T.times

    model = OceanTurb.TKEMassFlux.Model(
                      grid = OceanTurb.UniformGrid(N=N, H=H),
                   stepper = :BackwardEuler,
        eddy_diffusivities = OceanTurb.TKEMassFlux.RiDependentDiffusivities(),
                 constants = constants
    )

    # Coarse grain initial condition from LES and set equal
    # to initial condition of parameterization.
    model.solution.T.data[1:N] .= interior(T)[1, 1, :, 1]

    # Set boundary conditions
    FT = ds.metadata["temperature_flux"]
    ∂T∂z = ds.metadata["dθdz_deep"]
    model.bcs.T.top = OceanTurb.FluxBoundaryCondition(FT)
    model.bcs.T.bottom = OceanTurb.GradientBoundaryCondition(∂T∂z)

    solution = zeros(N, Nt)
    flux = zeros(N+1, Nt)

    for n in 1:Nt
        OceanTurb.run_until!(model, Δt, times[n])

        solution[:, n] .= model.solution.T[1:N]

        flux[:, n] .= OceanTurb.diffusive_flux(:T, model)[1:N+1]
        flux[N+1, n] = FT
    end

    return (T=solution, wT=flux)
end


text: rm(test_private_key; force = true, recursive = true)

ext: jl


text:
import .FTFont

mutable struct Font
	font::FTFont.Font
	model::Model
	vertexType::DataType
	textureUniform::Symbol
	charCount::Int

	Font() = new()
end

isvalid(font::Font) = isdefined(font, :font) && isvalid(font.model)

function init(font::Font, ftFont::FTFont.Font, shader::Shader; positionFunc::Function = identity, textureUniform::Symbol = :diffuseTexture, maxCharacters::Int = 2048)
	@assert isvalid(shader)

	fontName = FTFont.fontname(ftFont)
	texture = Texture()
	init(texture, shader.renderer, ftFont.bitmap; id = Symbol("Tex_" * fontName))

	material = Material(shader)
	setuniform(material, textureUniform, texture)
	setstate(material, AlphaBlendSrcAlpha())
	setstate(material, DepthStateDisabled())

	mesh = Mesh()
	init(mesh, shader.renderer, Array{shader.attribType}(undef, maxCharacters * 4), zeros(UInt16, maxCharacters * 6); positionFunc = positionFunc, id = Symbol("Mesh_" * fontName), usage = :dynamic)
	mesh.indexLength = 0

	font.font = ftFont
	font.vertexType = shader.attribType
	font.textureUniform = textureUniform
	font.model = Model(mesh, material)
	setbound(font.model, Shapes.Space{Float32}())
	font.charCount = 0

	finalizer(done, font)
	font
end

function done(font::Font)
	if isvalid(font)
		texture = getuniform(font.model.material, font.textureUniform)
		done(texture)
		done(font.model.mesh)
	end
end

function make_vertex(::Type{V}, pos, col, tex) where V
	t = (pos, col, tex)
	unsafe_load(Ptr{V}(pointer_from_objref(Ref(t))))
end

function drawchar(font::Font, pos::FTFont.Vec2{Float32}, bmp::Array{UInt8, 2}, box::FTFont.Rect{Int}, color::Color)
	@assert font.model.mesh.indexLength < length(font.model.mesh.indices)
	baseIndex = font.charCount * 6
	baseVertex = font.charCount * 4
	indices = font.model.mesh.indices
	vertices = font.model.mesh.vertices
	texWidth, texHeight = size(font.font.bitmap)
	boxF = FTFont.rect(Float32, box.min.x - 1, box.min.y - 1, box.max.x, box.max.y)
	boxFSize = size(boxF)

	vertType = eltype(vertices)
	vecColor = Vec4(color...)
	@assert sizeof(vertType) == sizeof(Tuple{Vec3, Vec4, Vec2})

	vertices[baseVertex+1] = make_vertex(vertType, Vec3(pos.x,			  pos.y,			  0), vecColor, Vec2(boxF.min.x / texWidth, boxF.min.y / texHeight))
	vertices[baseVertex+2] = make_vertex(vertType, Vec3(pos.x,			  pos.y + boxFSize.y, 0), vecColor, Vec2(boxF.min.x / texWidth, boxF.max.y / texHeight))
	vertices[baseVertex+3] = make_vertex(vertType, Vec3(pos.x + boxFSize.x, pos.y,			  0), vecColor, Vec2(boxF.max.x / texWidth, boxF.min.y / texHeight))
	vertices[baseVertex+4] = make_vertex(vertType, Vec3(pos.x + boxFSize.x, pos.y + boxFSize.y, 0), vecColor, Vec2(boxF.max.x / texWidth, boxF.max.y / texHeight))

	#=
	vertices[baseVertex+1] = Base.invokelatest(vertType, Vec3(pos.x,			  pos.y,			  0), vecColor, Vec2(boxF.min.x / texWidth, boxF.min.y / texHeight))
	vertices[baseVertex+2] = Base.invokelatest(vertType, Vec3(pos.x,			  pos.y + boxFSize.y, 0), vecColor, Vec2(boxF.min.x / texWidth, boxF.max.y / texHeight))
	vertices[baseVertex+3] = Base.invokelatest(vertType, Vec3(pos.x + boxFSize.x, pos.y,			  0), vecColor, Vec2(boxF.max.x / texWidth, boxF.min.y / texHeight))
	vertices[baseVertex+4] = Base.invokelatest(vertType, Vec3(pos.x + boxFSize.x, pos.y + boxFSize.y, 0), vecColor, Vec2(boxF.max.x / texWidth, boxF.max.y / texHeight))
	=#

	indices[baseIndex+1] = baseVertex
	indices[baseIndex+2] = baseVertex + 1
	indices[baseIndex+3] = baseVertex + 3
	indices[baseIndex+4] = baseVertex
	indices[baseIndex+5] = baseVertex + 3
	indices[baseIndex+6] = baseVertex + 2

	font.charCount += 1
end

function drawtext(font::Font, cursor::FTFont.TextCursor, s::AbstractString, color)
	@assert isvalid(font)
	drawFunc = (pos::FTFont.Vec2{Float32}, bmp::Array{UInt8, 2}, box::FTFont.Rect{Int})->drawchar(font, pos, bmp, box, color)
	FTFont.drawtext(drawFunc, font.font, cursor, s)
end

function cleartext(font::Font)
	@assert isvalid(font)
	font.model.mesh.indexLength = 0
	font.charCount = 0
end

function updatemodel(font::Font)
	indexLength = font.model.mesh.indexLength
	@assert font.charCount * 6 >= indexLength
	if font.charCount * 6 > indexLength
		updatebuffers(font.model.mesh, div(indexLength*4, 6)+1:font.charCount*4, indexLength+1:font.charCount*6)
		font.model.mesh.indexLength = font.charCount*6
	end
end

function add(renderer::Renderer, font::Font)
	updatemodel(font)
	add(renderer, font.model)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1ad3919a-1f79-11ec-1a44-d55ae420c8a6
begin
	using Pkg
	Pkg.activate("../")
	using nevanlinna_ac
	using Plots
	using DelimitedFiles
	using Printf
	"packages"
end

# ╔═╡ 8f75c874-0834-48a2-93a7-5521177826e4
η = 0.001

# ╔═╡ 8153e87c-fa69-44e0-8895-11ef2ffddf7c
omega = [i for i in range(0,5,length=1000)]

# ╔═╡ ad7da214-6055-4d57-b010-7fb4dbe96ea6
md"""
### J = 0.0; Γ = 1.0
"""

# ╔═╡ 8027d1ea-1188-4fd3-b1b9-99a9640c59f4
begin
	a0 =readdlm("../data/a_j_0.0.txt")
	g0 =make_input("../data/g_j_0.0.txt")
	"a0, g0"
end

# ╔═╡ ab3829e7-4ffa-422b-a62c-6f99e08bcc1d
isNevanlinnasolvable(g0)

# ╔═╡ fdad034c-6a53-43d3-b888-e38e2291b328
begin
	A0 = [spectrum_density(w,η,g0) for w in omega]
	g0m =make_input("../data/gm_j_0.0.txt")
	g0n =make_input("../data/gn_j_0.0.txt")
	A0mn = [spectrum_density(w,η,g0m) - spectrum_density(w,η,g0n) for w in omega]

	plot(a0[:,1], a0[:,2], lw=2, label="theory")
	plot!(omega, A0,line=(:dash,2), label="interpolate")
	plot!(omega, A0mn,line=(:dash,2), label="interpolate-2")
end

# ╔═╡ 0f5abb05-38e7-4009-93b2-f2ff7a4669f8
begin
	plot([g0.GF[i].ωn for i = 1:g0.number],[g0.GF[i].val |> imag for i = 1:g0.number],
		line=(2), label="imag")
	plot!([g0.GF[i].ωn for i = 1:g0.number],[g0.GF[i].val |> real for i = 1:g0.number],
		line=(2), label="real")
	plot!(xlabel="ωn", ylabel="GF value", legend=:bottomright)
end

# ╔═╡ 16aca10e-f308-488b-9b40-9dff09885d55
md"""
### J = Γ = 1.0
"""

# ╔═╡ bc61d277-7b0f-459f-bce9-8f0573a09f5b
begin
	a1 =readdlm("../data/a_j_1.0.txt")
	g1 =make_input("../data/g_j_1.0.txt")
	"a1 = A(ω), g1 = GF"
end

# ╔═╡ da4f1bcb-2ece-49ad-b9de-adf9e5d01e67
begin
	plot([g1.GF[i].ωn for i = 1:g1.number],[g1.GF[i].val |> imag for i = 1:g1.number],
		line=(2), label="imag")
	plot!([g1.GF[i].ωn for i = 1:g1.number],[g1.GF[i].val |> real for i = 1:g1.number],
		line=(2), label="real")
	plot!(xlabel="ωn", ylabel="GF value", legend=:bottomright)
end

# ╔═╡ 7b349ab2-f729-4a5c-99ae-6384ed3b3867
isNevanlinnasolvable(g1)

# ╔═╡ 1224a52a-55d0-4382-b678-544aa5ffe73d
A1 = [spectrum_density(w,η,g1) for w in omega]

# ╔═╡ 94d0a88f-f491-40c5-87f1-f4138f284fe3
begin
	g1m =make_input("../data/gm_j_1.0.txt")
	#g1n =make_input("../data/gn_j_1.0.txt")
	"∑_mn Cmn(iωn)/(Em - En)"
end

# ╔═╡ 43de4aed-4d19-4a8c-85f9-2ca5cd2cd2a8
begin
	plot([g1m.GF[i].ωn for i = 1:g1m.number],[g1m.GF[i].val |> imag for i = 1:g1m.number],
		line=(2), label="imag")
	plot!([g1m.GF[i].ωn for i = 1:g1m.number],[g1m.GF[i].val |> real for i = 1:g1m.number],
		line=(2), label="real")
	plot!(xlabel="ωn", ylabel="GF value", legend=:bottomright)
end

# ╔═╡ 96788e8e-0851-4a74-9214-9e6f7ae01852
A1m = [spectrum_density(w,η,g1m)*w for w in omega]

# ╔═╡ 54064a72-1a7b-4c7d-99f6-5ce292df7970
plot(omega, A1m)

# ╔═╡ 9415fb29-4294-47ec-94ed-28d778c504aa
begin
	n = 5
	b = 30 - n + 1
	g1ms = MasubaraGF(n, g1m.GF[b:30])
	A1ms = [spectrum_density(w,η,g1ms)*w for w in omega]
end

# ╔═╡ 1ab3f1f3-2c74-4717-91a7-f77863945d69
plot(omega, A1ms)

# ╔═╡ bd45648e-e070-40d2-b4e9-6b16f2c965af
isNevanlinnasolvable(g1m)

# ╔═╡ e150dbde-effb-417f-91ab-0ecc18d09945
function spectrum(ω::Real, η::Real, G::MasubaraGF)
    z = ω + 1.0im * η
    return -2im*evaluation(z, G) |> imag
end

# ╔═╡ 765f4fb6-43cb-44dc-b37a-e4eaf7e18197
begin
	ig1 =make_input("../data/ig_j_1.0.txt")
	"../data/ig_j_1.0.txt"
end

# ╔═╡ 75ed89a9-01f7-4587-bd21-c0ced4f8b819
ig1s = MasubaraGF(4, ig1.GF[27:30])

# ╔═╡ 13d87ae5-923c-476b-84b2-56bb44e46713
begin
	plot([ig1.GF[i].ωn for i = 1:ig1.number],[ig1.GF[i].val |> imag for i=1:ig1.number],
		line=(2), label="imag")
	plot!([ig1.GF[i].ωn for i = 1:ig1.number],
		  [ig1.GF[i].val |> real for i=1:ig1.number],
		  line=(2), label="real")
	plot!(xlabel="ωn", ylabel="iGF value", legend=:bottomright)
end

# ╔═╡ a6049bf4-715f-4507-a43e-4566d1fa6dc9
isNevanlinnasolvable(ig1)

# ╔═╡ 59a7814e-d8ea-436d-ae13-5139c1dea323
isNevanlinnasolvable(ig1s)

# ╔═╡ 23de7386-30c1-406f-89bc-d9cb3ae47868
begin
	g1n =make_input("../data/gn_j_1.0.txt")
	"iωn * G(iωn)"
end

# ╔═╡ 85e0fbba-80e2-4f14-836a-3579e23becaa
isNevanlinnasolvable(g1n)

# ╔═╡ 8ec7f3b6-b45e-448e-b7c6-fe12c9309a2b
begin
	plot([g1n.GF[i].ωn for i = 1:g1n.number],[g1n.GF[i].val |> imag for i = 1:g1n.number],
		line=(2), label="imag")
	plot!([g1n.GF[i].ωn for i = 1:g1n.number],[g1n.GF[i].val |> real for i = 1:g1n.number],
		line=(2), label="real")
	plot!(xlabel="ωn", ylabel="GF value", legend=:bottomright)
end

# ╔═╡ d710e801-cb4b-4ff5-a125-13578c1fafb5
A1ns = [spectrum_density(w,η,g1n) for w in omega]

# ╔═╡ 758d84d0-f844-4c8a-9c49-f6ae20eec7b6
begin

	#plot!(omega, A1,line=(:dash,2), label="interpolate")
	#plot!(omega, A1m,line=(2), label="C/dE")
	plot(omega, A1ms,line=(3), label="interpolate")
	plot!(a1[:,1], a1[:,2], line=(:black,:dash,3), label="theory")
	#plot!(omega, A1n,line=(:dash,2), label="C*ω")
	plot!(title="g=1.0, β=20")
	plot!(ylim=(-0.1,10))
	plot!(xlabel="ω", ylabel="A(ω)",
	xtickfont=font(12),
    ytickfont=font(12),
    guidefont=font(12),
    legendfont=font(12))
end

# ╔═╡ 02fb2a24-5aa8-4e6e-9bff-68cefaa21f3b
md"""
$A(x) = \frac{x}{√(2π)} e^{-x^2 /2 }$
"""

# ╔═╡ 11f6a78a-4bf1-4a12-8176-43bcd9f1516a
begin
	A = readdlm("../data/ag1_bose.txt")
	"../data/ag1_bose.txt"
end

# ╔═╡ 791a5b9b-b2b3-4dcb-96c6-9d4ff5a0b083
begin
	g = make_input("../data/ig1_bose.txt")
	"../data/ig1_bose.txt"
end

# ╔═╡ 54a57c63-ff28-44b6-946c-6baeb5b5c56e
isNevanlinnasolvable(g1)

# ╔═╡ c243728f-64d2-4994-b686-44ef14b5dbaa
a = [spectrum_density(w,η,g)*w for w in omega]

# ╔═╡ 1600847a-de36-4c4a-bf74-d868b3bedaad
begin
	plot(A[:,1], A[:,2], lw=2, label="theory")
	plot!(omega, a,line=(:dash,2), label="interpolate")
	plot!(xlim=(0,5))
end

# ╔═╡ 286f864c-ab38-4ae0-8be6-66dea42c0056
begin
	plot([g.GF[i].ωn for i = 1:g.number],[g.GF[i].val |> imag for i = 1:g.number],
		line=(2), label="imag")
	plot!([g.GF[i].ωn for i = 1:g.number],[g.GF[i].val |> real for i = 1:g.number],
		line=(2), label="real")
	plot!(xlabel="ωn", ylabel="GF value", legend=:bottomright)
end

# ╔═╡ 84a0e459-36c9-462e-bd3c-14f24cb486c0
md"""
$A(x) = \frac{x^3}{√(2π)} e^{-x^2 /2 }$
"""

# ╔═╡ 3ca7de5c-9a31-41db-bad2-90b46c7111c5
begin
	A3 = readdlm("../data/ag3_bose.txt")
	"../data/ag3_bose.txt"
end

# ╔═╡ f21b14b4-fab8-40db-b76c-015f9f3c9e8d
begin
	g3 = make_input("../data/ig3_bose.txt")
	"../data/ig3_bose.txt"
end

# ╔═╡ f8e48a07-5d9f-486c-8599-08ea50747508
isNevanlinnasolvable(g3)

# ╔═╡ 9a3aeef6-d92f-49ee-ba9b-90e9f79d648b
a3 = [spectrum_density(w,η,g3)*w for w in omega]

# ╔═╡ 423c967d-3cfc-4650-9a1a-6d5dffedec0a
begin
	plot(A3[:,1], A3[:,2], lw=2, label="theory")
	plot!(omega, a3,line=(:dash,2), label="interpolate")
	plot!(xlim=(0,5))
end

# ╔═╡ d46b79d4-bfc7-4f85-a801-409e414e8439
begin
	plot([g3.GF[i].ωn for i = 1:g3.number],[g3.GF[i].val |> imag for i = 1:g3.number],
		line=(2), label="imag")
	plot!([g3.GF[i].ωn for i = 1:g3.number],[g3.GF[i].val |> real for i = 1:g3.number],
		line=(2), label="real")
	plot!(xlabel="ωn", ylabel="GF value", legend=:bottomright)
end

# ╔═╡ d44b60bb-6cff-4e81-9682-7bcf5643717e
pwd()

# ╔═╡ Cell order:
# ╟─8f75c874-0834-48a2-93a7-5521177826e4
# ╠═8153e87c-fa69-44e0-8895-11ef2ffddf7c
# ╟─ad7da214-6055-4d57-b010-7fb4dbe96ea6
# ╟─8027d1ea-1188-4fd3-b1b9-99a9640c59f4
# ╟─ab3829e7-4ffa-422b-a62c-6f99e08bcc1d
# ╠═fdad034c-6a53-43d3-b888-e38e2291b328
# ╟─0f5abb05-38e7-4009-93b2-f2ff7a4669f8
# ╟─16aca10e-f308-488b-9b40-9dff09885d55
# ╠═bc61d277-7b0f-459f-bce9-8f0573a09f5b
# ╟─da4f1bcb-2ece-49ad-b9de-adf9e5d01e67
# ╟─7b349ab2-f729-4a5c-99ae-6384ed3b3867
# ╠═1224a52a-55d0-4382-b678-544aa5ffe73d
# ╠═94d0a88f-f491-40c5-87f1-f4138f284fe3
# ╟─43de4aed-4d19-4a8c-85f9-2ca5cd2cd2a8
# ╠═96788e8e-0851-4a74-9214-9e6f7ae01852
# ╠═54064a72-1a7b-4c7d-99f6-5ce292df7970
# ╠═9415fb29-4294-47ec-94ed-28d778c504aa
# ╟─1ab3f1f3-2c74-4717-91a7-f77863945d69
# ╟─bd45648e-e070-40d2-b4e9-6b16f2c965af
# ╠═e150dbde-effb-417f-91ab-0ecc18d09945
# ╠═765f4fb6-43cb-44dc-b37a-e4eaf7e18197
# ╠═75ed89a9-01f7-4587-bd21-c0ced4f8b819
# ╟─13d87ae5-923c-476b-84b2-56bb44e46713
# ╠═85e0fbba-80e2-4f14-836a-3579e23becaa
# ╟─a6049bf4-715f-4507-a43e-4566d1fa6dc9
# ╟─59a7814e-d8ea-436d-ae13-5139c1dea323
# ╠═23de7386-30c1-406f-89bc-d9cb3ae47868
# ╟─8ec7f3b6-b45e-448e-b7c6-fe12c9309a2b
# ╠═d710e801-cb4b-4ff5-a125-13578c1fafb5
# ╠═758d84d0-f844-4c8a-9c49-f6ae20eec7b6
# ╟─02fb2a24-5aa8-4e6e-9bff-68cefaa21f3b
# ╟─11f6a78a-4bf1-4a12-8176-43bcd9f1516a
# ╟─791a5b9b-b2b3-4dcb-96c6-9d4ff5a0b083
# ╟─54a57c63-ff28-44b6-946c-6baeb5b5c56e
# ╠═c243728f-64d2-4994-b686-44ef14b5dbaa
# ╠═1600847a-de36-4c4a-bf74-d868b3bedaad
# ╠═286f864c-ab38-4ae0-8be6-66dea42c0056
# ╟─84a0e459-36c9-462e-bd3c-14f24cb486c0
# ╟─3ca7de5c-9a31-41db-bad2-90b46c7111c5
# ╟─f21b14b4-fab8-40db-b76c-015f9f3c9e8d
# ╟─f8e48a07-5d9f-486c-8599-08ea50747508
# ╟─9a3aeef6-d92f-49ee-ba9b-90e9f79d648b
# ╠═423c967d-3cfc-4650-9a1a-6d5dffedec0a
# ╟─d46b79d4-bfc7-4f85-a801-409e414e8439
# ╠═1ad3919a-1f79-11ec-1a44-d55ae420c8a6
# ╠═d44b60bb-6cff-4e81-9682-7bcf5643717e


text:
using NeuralThermalODE

# This is either the raw HTC or a spline fit.
htc_file = "checkpoint_htc_360.csv"
htc_recovered = read_htc(htc_file)

# Experimental LN2 HTC
exp_file = "ln2.csv"
exp_recovered = read_htc(exp_file)

const testfile = "Fig_67.csv_Result.csv"
const test_data = read_data(testfile)

const N = 10
const t = 0.00635
const dx = t / N

const sinkT = test_data[end, 2]
const ambT = test_data[1, 2]
const u20 = [ambT for i in 1:N]

tspan = (Float64(test_data[1, 1]), Float64(test_data[end, 1]))

# Solve ODE.
const p = Vector(htc_recovered[2])

prob = ODEProblem(heat_transfer, u20, tspan, p, saveat = 1.0)
sol = solve(prob, Tsit5(), saveat = 1.0)
#plot(sol.t[:], sol[1, :], label = "Exposed")
#plot!(sol.t[:], sol[end, :], label = "Unexposed")

exposedT = sol[1, :]
vapourT = (exposedT .+ sinkT) / 2
ΔT = abs.(sinkT .- vapourT)
plot(xlims = (1, 120))
#plot!(xscale = :log, yscale = :log)
plot!(xlabel = "\\Delta T", ylabel = "HTC W/mK")
plot!(ΔT, p, label = "Recovered")
plot!(exp_recovered[1], exp_recovered[2], label = "Small scale experiment")

max_stars_repo_name: NikStoyanov/NeuralThermalODE

max_issues_repo_name: NikStoyanov/NeuralThermalODE

max_forks_repo_name: NikStoyanov/NeuralThermalODE


text:
include("maze.jl")

pkgs = ["ModernGL", "GLAbstraction", "GeometryTypes", "GLWindow", "Images", "FileIO", "Reactive", "ImageMagick"]
for pkg in pkgs; Pkg.installed(pkg) == nothing && Pkg.add(pkg); end

using ModernGL, GeometryTypes, GLAbstraction, GLWindow, Images, FileIO, Reactive, ArgParse

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin
		"--epochs"
			help = "number of epochs"
			default = 20
			arg_type=Int
		"--eps"
			help = "epsilon parameter for the random behaviour"
			default = 0.05
			arg_type = Float64
		"--gamma"
			help = "the discount factor"
			default = 0.9
			arg_type = Float64
		"--alpha"
			help = "learning step"
			default = 0.8
			arg_type = Float64
		"--dims"
			help = "dimensions of the maze"
			nargs='+'
			default = [3,3]
			arg_type = Int

	end
	return parse_args(s)
end

function rectangle(start, width, height, pos, cols, elms)
	indx = length(pos)

	push!(pos, start)
	push!(pos, (start[1] + width, start[2]))
	push!(pos, (start[1], start[2] - height))
	push!(pos, (start[1] + width, start[2] - height))

	push!(elms, (indx, indx+1, indx+2))
	push!(elms, (indx+1, indx+2, indx+3))

	for i=1:4; push!(cols, (0.0, 0.0, 0.0)); end
end

function agent_model(start, width, height, pos, cols, elms)
	indx = length(pos)
	push!(pos, (start[1] + width*0.5, start[2] - height*0.25))
	push!(pos, (start[1] + width*0.75, start[2] - height*0.75))
	push!(pos, (start[1] + width*0.25, start[2] - height*0.75))

	push!(cols, (1.0, 0.0, 0.0))
	push!(cols, (0.0, 0.0, 1.0))
	push!(cols, (0.0, 0.0, 1.0))

	push!(elms, (indx, indx+1, indx+2))
end

function paint(start, width, height, pos, cols, elms, clr)
	indx = length(pos)
	push!(pos, start)
	push!(pos, (start[1] + width, start[2]))
	push!(pos, (start[1], start[2] - height))
	push!(pos, (start[1] + width, start[2] - height))

	push!(elms, (indx, indx+1, indx+2))
	push!(elms, (indx+1, indx+2, indx+3))

	for i=1:4; push!(cols, clr); end
end

function key_callback(window, key, scancode, action, mode)
	global user
	if key == GLFW.KEY_ESCAPE && action == GLFW.PRESS
		GLFW.SetWindowShouldClose(window, true)
	else
		if key == GLFW.KEY_L && action == GLFW.PRESS
			left()
		elseif key == GLFW.KEY_R && action == GLFW.PRESS
			right()
		elseif key == GLFW.KEY_M && action == GLFW.PRESS
			move()
		elseif key == GLFW.KEY_S && action == GLFW.PRESS
			reset()
		elseif key == GLFW.KEY_U && action == GLFW.PRESS
			user = !user
		elseif key == GLFW.KEY_1 && action == GLFW.PRESS
			global speed
			speed = speed * 0.5
		elseif key == GLFW.KEY_2 && action == GLFW.PRESS
			global speed
			speed = speed * 2
		elseif key == GLFW.KEY_T && action == GLFW.PRESS
			global train
			train = !train
		end
	end
end

function reset()
	global direction
	global ro
	global a_ro
	global s_ro
	global f_ro
	global signal
	global window
	global a_model
	global env
	hm, wm,_ = size(env.maze)
	w = 2.0 / wm
	h = 2.0 / hm

	a_ro[:view] = translationmatrix(Vec{3, Float32}((-(wm / 2 - 0.5)*w, (hm/2 + 0.5 - env.start[1])*h , 0.0f0))) * scalematrix(Vec{3, Float32}((1.0/wm, 1.0/hm, 1.0f0)))

	push!(signal, rotationmatrix_z(deg2rad((direction-1)*90)))
	direction = 1

	Reactive.run_till_now()
	glClear(GL_COLOR_BUFFER_BIT)
	render(ro)
	render(s_ro)
	render(f_ro)
	render(a_ro)
	GLFW.SwapBuffers(window)
	sleep(0.01)
end

function rot(pos=1.0)
	global direction
	global ro
	global a_ro
	global s_ro
	global f_ro
	global signal
	global window
	global speed

	if pos == 1.0
		direction = direction == 1 ? 4 : direction - 1
	else
		direction = direction == 4 ? 1 : direction + 1
	end

	for i=1:15
		push!(signal, rotationmatrix_z(deg2rad(6*pos)))
		Reactive.run_till_now()
		glClear(GL_COLOR_BUFFER_BIT)
		render(ro)
		render(s_ro)
		render(f_ro)
		render(a_ro)
		GLFW.SwapBuffers(window)
		sleep(speed)
	end

end

left() = rot(1.0)
right() = rot(-1.0)

function move()
	global env
	global direction
	global signal
	global ro
	global a_ro
	global s_ro
	global f_ro
	global window
	global speed

	hm, wm,_ = size(env.maze)
	w = 2.0 / wm
	h = 2.0 / hm
	wdx = direction == 2 ? 1.0 : direction == 4 ? -1.0 : 0.0
	hdx = direction == 1 ? 1.0 : direction == 3 ? -1.0 : 0.0
	for i=1:15
		a_ro[:view] = translationmatrix(Vec{3, Float32}((w/15.0 * wdx, h/15.0 * hdx,0.0f0))) * a_ro[:view]
		glClear(GL_COLOR_BUFFER_BIT)
		render(ro)
		render(s_ro)
		render(f_ro)
		render(a_ro)
		GLFW.SwapBuffers(window)
		sleep(speed)
	end
end

function create_maze(wm, hm)
	imgstart = load("start.png")
	imgend = load("checkered.png")

	w = 2.0 / wm
	h = 2.0 / hm

	positions = Point{2, Float32}[]
	clrs = Vec3f0[]
	elements = Face{3, UInt32, -1}[]

	a_positions = Point{2, Float32}[]
	a_clrs = Vec3f0[]
	a_elements = Face{3, UInt32, -1}[]

	agent_model((-1, 1), 2.0, 2.0, a_positions, a_clrs, a_elements)

	for i=1:wm
		for j=1:hm
			start = (-1 + (i-1)*w, 1 - (j-1)*h)
			if env.maze[j, i, 1] == 0
				rectangle(start, w, h / 10.0, positions, clrs, elements)
			end
			if env.maze[j, i, 2] == 0 && i == wm
				s = (start[1] + w * 0.9, start[2])
				rectangle(s, w / 10.0, h, positions, clrs, elements)
			end
			if env.maze[j, i, 3] == 0 && j == hm
				s = (start[1], start[2] - h * 0.9)
				rectangle(s, w, h / 10.0, positions, clrs, elements)
			end
			if env.maze[j, i, 4] == 0
				rectangle(start, w / 10.0, h, positions, clrs, elements)
			end
		end
	end

	s_vertex_positions = Point{2,Float32}[(-1.0, 1.0), (1.0, 1.0), (1.0, -1.0), (-1.0, -1.0)]

	# The colors assigned to each vertex
	vertex_colors = Vec3f0[(1, 1, 1),      # top-left
		(1, 1, 1),                     # top-right
		(1, 1, 1),                     # bottom-right
		(1, 1, 1)]                     # bottom-left

	# The texture coordinates of each vertex
	vertex_texcoords = Vec2f0[(0, 0),
	                          (1, 0),
				  (1, 1),
				  (0, 1)]

	vertex_elements = Face{3,UInt32,-1}[(0,1,2),          # the first triangle
	                             (2,3,0)]          # the second triangle

	vertex_source= vert"""
	# version 150
	in vec2 position;
	in vec3 color;

	uniform mat4 model;
	uniform mat4 view;
	uniform mat4 proj;

	out vec3 Color;
	void main()
	{
	Color = color;
	gl_Position = proj * view * model * vec4(position,0.0, 1.0);
	}
	"""

	tex_vertex_source= vert"""
	# version 150
	in vec2 position;
	in vec3 color;
	in vec2 texcoord;

	uniform mat4 model;
	uniform mat4 view;
	uniform mat4 proj;

	out vec3 Color;
	out vec2 Texcoord;

	void main()
	{
	Color = color;
	Texcoord = texcoord;
	gl_Position = proj * view * model * vec4(position,0.0, 1.0);
	}
	"""

	fragment_source = frag"""
	# version 150
	in vec3 Color;
	out vec4 outColor;
	void main()
	{
	outColor = vec4(Color, 1.0);
	}
	"""

	tex_fragment_source = frag"""
	# version 150
	in vec3 Color;
	in vec2 Texcoord;

	out vec4 outColor;

	uniform sampler2D tex;
	void main()
	{
	outColor = texture(tex, Texcoord) * vec4(Color, 1.0);
	}
	"""

	global signal = Signal(rotate(0f0, Vec((0,0,1f0))))
	model = rotate(0f0, Vec((0,0,1f0)))
	view = rotate(0f0, Vec((0,0,1f0)))
	global a_model = foldp(*, rotate(0f0, Vec((0,0,1f0))), signal)
	a_view = translationmatrix(Vec((-(wm / 2 - 0.5)*w, (hm/2 + 0.5 - env.start[1])*h , 0.0f0))) * scalematrix(Vec((1.0/wm, 1.0/hm, 1.0f0)))

	proj = orthographicprojection(Float32, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0)

	bufferdict = Dict(:position=>GLBuffer(positions),
	:color=>GLBuffer(clrs),
	:indexes=>indexbuffer(elements),
	:model=>model,
	:view=>view,
	:proj=>proj,
	)

	a_bufferdict = Dict(:position=>GLBuffer(a_positions),
	:color=>GLBuffer(a_clrs),
	:indexes=>indexbuffer(a_elements),
	:model=>a_model,
	:view=>a_view,
	:proj=>proj
	)

	s_view = translationmatrix(Vec((-(wm / 2 - 0.5)*w, (hm/2 + 0.5 - env.start[1])*h , 0.0f0))) * scalematrix(Vec((1.0/(wm+3), 1.0/(hm+3), 1.0f0)))
	s_bufferdict = Dict(:position=>GLBuffer(s_vertex_positions),
	:color=>GLBuffer(vertex_colors),
	:indexes=>indexbuffer(vertex_elements),
	:model=>model,
	:view=>s_view,
	:proj=>proj,
	:tex=>Texture(data(imgstart)),
	:texcoord=>GLBuffer(vertex_texcoords),
	)

	f_view = translationmatrix(Vec((-(wm / 2 + 0.5 - env.goal[2])*w, (hm/2 + 0.5 - env.goal[1])*h , 0.0f0))) * scalematrix(Vec((1.0/(wm+3), 1.0/(hm+3), 1.0f0)))
	f_bufferdict = Dict(:position=>GLBuffer(s_vertex_positions),
	:color=>GLBuffer(vertex_colors),
	:indexes=>indexbuffer(vertex_elements),
	:model=>model,
	:view=>f_view,
	:proj=>proj,
	:tex=>Texture(data(imgend)),
	:texcoord=>GLBuffer(vertex_texcoords),
	)

	global ro = std_renderobject(bufferdict, LazyShader(vertex_source, fragment_source))
	global a_ro = std_renderobject(a_bufferdict, LazyShader(vertex_source, fragment_source))
	global s_ro = std_renderobject(s_bufferdict, LazyShader(tex_vertex_source, tex_fragment_source))
	global f_ro = std_renderobject(f_bufferdict, LazyShader(tex_vertex_source, tex_fragment_source))
end

function apply(action)
	global direction
	if action == MazeAction(UP)
		if direction == 1
			move()
		elseif direction == 2
			left()
			move()
		elseif direction == 3
			right()
			right()
			move()
		else
			right()
			move()
		end
	elseif action == MazeAction(RIGHT)
		if direction == 1
			right()
			move()
		elseif direction == 2
			move()
		elseif direction == 3
			left()
			move()
		else
			left()
			left()
			move()
		end
	elseif action == MazeAction(DOWN)
		if direction == 1
			right()
			right()
			move()
		elseif direction == 2
			right()
			move()
		elseif direction == 3
			move()
		else
			left()
			move()
		end

	else
		if direction == 1
			left()
			move()
		elseif direction == 2
			right()
			right()
			move()
		elseif direction == 3
			right()
			move()
		else
			move()
		end
	end
end

function main()
	args = parse_commandline()
	hm, wm = args["dims"]
	global env = MazeEnv((hm,wm))
	agent = QLearner(env;ε=args["eps"], α=args["alpha"], Ɣ=args["gamma"])

	global direction = 1
	global train = true
	global user = false

	global window = create_glcontext("Maze Solving", resolution=(800, 600))
	global speed = 0.01

	vao = glGenVertexArrays()
	glBindVertexArray(vao)

	#Creating the environment
	println("Creating the environment")
	create_maze(wm, hm)

	global ro
	global a_ro
	GLFW.SetKeyCallback(window, key_callback)

	glClearColor(1.0,1.0,1.0,1.0)
	state = getInitialState(env)
	totalRewards = 0.0
	numOfStates = 1
	epoch = 1
	numOfEpochs = args["epochs"]

	while !GLFW.WindowShouldClose(window)
		glClear(GL_COLOR_BUFFER_BIT)
		if !user
			if isTerminal(state, env)
				println("Epoch: $epoch, totalReward: $totalRewards")
				state = getInitialState(env)
				totalRewards = 0.0
				numOfStates = 1
				epoch += 1
				if epoch > numOfEpochs
					GLFW.SetWindowShouldClose(window, true)
				end
				reset()
			end
			action = play(agent, state, env; learn=train)
			apply(action)
			state, reward = transfer(env, state, action)
			observe(agent, state, reward, env; learn=train)
			totalRewards += reward
			numOfStates += 1
		else
			render(ro)
			render(s_ro)
			render(f_ro)
			render(a_ro)
		end
		GLFW.SwapBuffers(window)
		GLFW.PollEvents()
	end
end

main()

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
abstract type AbstractSubSequence end

struct Subsequence{VT <: AbstractVector} <: AbstractSubSequence
    seq::VT
    onset::Int
    type::String
end

struct Motif <: AbstractSubSequence
    seqs::Vector{Subsequence}
end

const SubSeqType = Union{AbstractSubSequence, AbstractVector{<:AbstractSubSequence}}

onsets(m::Motif) = getfield.(m.seqs, :onset)
seqs(m::Motif) = getfield.(m.seqs, :seq)
seqlength(m::Motif) = length(m.seqs[1].seq)
subseqtype(::Union{Motif, Vector{<:Motif}}) = "Motif"

onsets(m::Subsequence) = m.onset
seqs(m::Subsequence) = [m.seq]
seqlength(m::Subsequence) = length(m.seq)
subseqtype(m::Subsequence) = m.type

onsets(m::Array{<:Subsequence}) = onsets.(m)
seqs(m::Array{<:Subsequence}) = getfield.(m, :seq)
seqlength(m::Array{<:Subsequence}) = length(m[1].seq)
subseqtype(m::Array{<:Subsequence}) = subseqtype(m[1])

function motifs(p::Profile, k, found_motifs = Motif[]; r=2, th = p.m, dist = ZEuclidean())
    length(found_motifs) == k && return found_motifs
    m = p.m
    P = copy(p.P)
    remove_found_motifs!(P, found_motifs, th)
    d, i = findmin(P)
    distance_profile!(P, dist, getwindow(p.T, m, i), p.T)
    remove_found_motifs!(P, found_motifs, th)
    perm = sortperm(P)
    onsets = [perm[1]]
    j = 2
    while P[perm[j]] < (d + 1e-5)*r
        all(abs.(perm[j] .- onsets) .> th) && push!(onsets, perm[j])
        j += 1
    end
    sort!(onsets)
    push!(found_motifs, Motif(Subsequence(p, onsets, "Motif")))
    motifs(p::Profile, k, found_motifs; r=r, th=th, dist=dist)
end

function Subsequence(p::Profile, onsets, type::String = "Subsequence")
    m    = p.m
    inds = 0:m-1
    n    = length(onsets)
    vecs = map(1:n) do i
        p.T[inds .+ onsets[i]] # Don't use getwindow
    end
    Subsequence.(vecs, onsets, type)
end

function remove_found_motifs!(P, found_motifs, th)
    for mo in found_motifs
        for o in onsets(mo)
            inds = max(1, o-th):min(length(P), o+th)
            P[inds] .= typemax(eltype(P))
        end
    end
end

function anomalies(p::Profile, k)
    onsets = partialsortperm(p.P, 1:k, rev=true)
    Subsequence(p, copy(onsets), "Anomaly")
end


text:
function COSMO_parse_Pi(param::SingleParam{String})
    Vec = Vector{Vector{Float64}}(undef,0)
    n = length(param.components)
    default = Float64[]
    for i in 1:n
        if !param.ismissingvalues[i]
            push!(Vec,parse.(Float64,split(param.values[i])))
        else
            push!(Vec,default)
        end
    end
    SingleParam(param.name,param.components,Vec,param.ismissingvalues,param.sourcecsvs,param.sources)
end

cosmo_tol(_new,_old) = mapreduce((x,y) -> abs(x/y -1.0),+,_new,_old)


text:
using Distributed, SharedArrays, Test
@everywhere using BrainWave

sampRate = 8192; nTaps=512
nChans=8
rec, evtTab = simulateRecording(nChans=nChans, dur=500, sampRate=sampRate)
rec2 = deepcopy(rec)
chans = collect(1:nChans)

filterContinuous!(rec, sampRate, "highpass", nTaps, [30], channels=chans, transitionWidth=0.2)

sharedRec = convert(SharedArray, rec2);
filterContinuous!(sharedRec, sampRate, "highpass", nTaps, [30], channels=chans, transitionWidth=0.2);
rec2 = convert(Array, sharedRec)

@test unique(rec2 .== rec)[1] .== true


text:
################################################################################
#
#          NfOrd/Elem.jl : Elements of orders of number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export ==, +, -, *, ^, add!, conjugates_arb, conjugates_arb_log, discriminant,
       divexact, elem_in_nf, coordinates, isone, iszero, minkowski_map, mod,
       mul!, norm, one, parent, powermod, rand, rand!, representation_matrix,
       show, trace, t2, zero

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(x::NfAbsOrdElem{S, T}, dict::IdDict) where {S, T}
  z = parent(x)()
  z.elem_in_nf = Base.deepcopy_internal(x.elem_in_nf, dict)
  if x.has_coord
    z.has_coord = true
    z.coordinates = Base.deepcopy_internal(x.coordinates, dict)::Vector{fmpz}
  end
  return z
end

################################################################################
#
#  Conversion from matrix
#
################################################################################

function elem_from_mat_row(O::NfAbsOrd, M::fmpz_mat, i::Int, d::fmpz = fmpz(1))
  return O(fmpz[M[i, j] for j=1:degree(O)])
end

################################################################################
#
#  Parent object overloading
#
################################################################################

@doc Markdown.doc"""
      (O::NfOrd)(a::nf_elem, check::Bool = true) -> NfAbsOrdElem

Given an element $a$ of the ambient number field of $\mathcal O$, this
function coerces the element into $\mathcal O$. It will be checked that $a$
is contained in $\mathcal O$ if and only if `check` is `true`.
"""
(O::NfAbsOrd{S, T})(a::T, check::Bool = true) where {S, T} = begin
  if nf(O) !== parent(a)
    throw(error("Underlying number fields not equal"))
  end
  if check
    (x, y) = _check_elem_in_order(a,O)
    !x && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(a), y)
  else
    return NfAbsOrdElem(O, deepcopy(a))
  end
end

@doc Markdown.doc"""
      (O::NfOrd)(a::NfAbsOrdElem, check::Bool = true) -> NfAbsOrdElem

Given an element $a$ of some order in the ambient number field of
$\mathcal O$, this function coerces the element into $\mathcal O$. It
will be checked that $a$ is contained in $\mathcal O$ if and only if
`check` is `true`.
"""
(O::NfAbsOrd{S, T})(a::NfAbsOrdElem{S, T}, check::Bool = true) where {S, T} = begin
  b = nf(parent(a))(a)
  if check
    (x, y) = _check_elem_in_order(b,O)
    !x && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(b), y)
  else
    return NfAbsOrdElem(O, deepcopy(b))
  end
end

(O::NfAbsOrd{S, T})(a::T, arr::Vector{fmpz}, check::Bool = false) where {S, T} = begin
  if check
    (x, y) = _check_elem_in_order(a,O)
    (!x || arr != y ) && error("Number field element not in the order")
    return NfAbsOrdElem(O, deepcopy(a), y)
  else
    return NfAbsOrdElem(O, deepcopy(a), deepcopy(arr))
  end
end

@doc Markdown.doc"""
      (O::NfOrd)(a::Union{fmpz, Integer}) -> NfAbsOrdElem

Given an element $a$ of type `fmpz` or `Integer`, this
function coerces the element into $\mathcal O$. It will be checked that $a$
is contained in $\mathcal O$ if and only if `check` is `true`.
"""
(O::NfAbsOrd)(a::Union{fmpz, Integer}) = begin
  return NfAbsOrdElem(O, nf(O)(a))
end

@doc Markdown.doc"""
      (O::NfOrd)(arr::Array{fmpz, 1})

Returns the element of $\mathcal O$ with coefficient vector `arr`.
"""
(O::NfAbsOrd)(arr::Array{fmpz, 1}) = begin
  return NfAbsOrdElem(O, deepcopy(arr))
end

@doc Markdown.doc"""
      (O::NfOrd)(arr::Array{Integer, 1})

Returns the element of $\mathcal O$ with coefficient vector `arr`.
"""
(O::NfAbsOrd)(arr::Array{S, 1}) where {S <: Integer} = begin
  return NfAbsOrdElem(O, deepcopy(arr))
end

@doc Markdown.doc"""
      (O::NfOrd)() -> NfAbsOrdElem

This function constructs a new element of $\mathcal O$ which is set to $0$.
"""
(O::NfAbsOrd)() = NfAbsOrdElem(O)

################################################################################
#
#  Parent
#
################################################################################

@doc Markdown.doc"""
    parent(a::NfAbsOrdElem) -> NfOrd

Returns the order of which $a$ is an element.
"""
parent(a::NfAbsOrdElem) = a.parent

################################################################################
#
#  Parent check
#
################################################################################

function check_parent(x::NfAbsOrdElem{S, T}, y::NfAbsOrdElem{S, T}) where {S, T}
  return parent(x) === parent(y)
end

################################################################################
#
#  Element in number field
#
################################################################################

@doc Markdown.doc"""
    elem_in_nf(a::NfAbsOrdElem) -> nf_elem

Returns the element $a$ considered as an element of the ambient number field.
"""
function elem_in_nf(a::NfAbsOrdElem; copy::Bool = true)
  if isdefined(a, :elem_in_nf)
    if copy
      return deepcopy(a.elem_in_nf)
    else
      return a.elem_in_nf
    end
  end
  error("Not a valid order element")
end

_elem_in_algebra(a::NfAbsOrdElem; copy::Bool = true) = elem_in_nf(a, copy = copy)

################################################################################
#
#  "Assure" functions for fields
#
################################################################################

function assure_has_coord(a::NfAbsOrdElem)
  if a.has_coord
    return nothing
  else
    (x, y) = _check_elem_in_order(a.elem_in_nf, parent(a))
    !x && error("Not a valid order element")
    a.coordinates = y
    a.has_coord = true
    return nothing
  end
end

################################################################################
#
#  Coordinates
#
################################################################################

@doc Markdown.doc"""
    coordinates(a::NfAbsOrdElem) -> Array{fmpz, 1}

Returns the coefficient vector of $a$.
"""
function coordinates(a::NfAbsOrdElem; copy::Bool = true)
  assure_has_coord(a)
  @hassert :NfOrd 2 a == dot(a.coordinates, basis(parent(a), copy = false))
  if copy
    return deepcopy(a.coordinates)
  else
    return a.coordinates
  end
end

################################################################################
#
#  Discriminant
#
################################################################################

@doc Markdown.doc"""
    discriminant(B::Array{NfAbsOrdElem, 1}) -> fmpz

Returns the discriminant of the family $B$.
"""
function discriminant(B::Array{NfAbsOrdElem{S, T}, 1}) where {S, T}
  length(B) == 0 && error("Number of elements must be non-zero")
  length(B) != degree(parent(B[1])) &&
        error("Number of elements must be $(degree(parent(B[1])))")
  O = parent(B[1])
  A = zero_matrix(FlintZZ, degree(O), degree(O))
  for i in 1:degree(O)
    el = tr(B[i]^2)
    A[i, i] = el
    for j in 1:degree(O)
      el = tr(B[i] * B[j])
      A[i, j] = el
      A[j, i] = el
    end
  end
  return det(A)
end

################################################################################
#
#  Hashing
#
################################################################################

Base.hash(x::NfAbsOrdElem, h::UInt) = Base.hash(x.elem_in_nf, h)

################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(x::NfAbsOrdElem, y::NfAbsOrdElem) -> Bool

Returns whether $x$ and $y$ are equal.
"""
 ==(x::NfAbsOrdElem, y::NfAbsOrdElem) = parent(x) === parent(y) &&
                                            x.elem_in_nf == y.elem_in_nf

################################################################################
#
#  Characteristic and minimal polynomial
#
################################################################################
@doc Markdown.doc"""
    charpoly(a::NfAbsOrdElem) -> fmpz_poly

    charpoly(a::NfAbsOrdElem, FlintZZ) -> fmpz_poly
The characteristic polynomial of $a$.
"""
function charpoly(a::NfAbsOrdElem, Zx::FmpzPolyRing = FmpzPolyRing(FlintZZ, :x, false))
  return Zx(charpoly(elem_in_nf(a)))
end

@doc Markdown.doc"""
    minpoly(a::NfAbsOrdElem) -> fmpz_poly
The minimal polynomial of $a$.
"""
function minpoly(a::NfAbsOrdElem, Zx::FmpzPolyRing = FmpzPolyRing(FlintZZ, :x, false))
  return Zx(minpoly(elem_in_nf(a)))
end

################################################################################
#
#  Special elements
#
################################################################################

@doc Markdown.doc"""
    zero(O::NfOrd) -> NfAbsOrdElem

Returns the zero element of $\mathcal O$.
"""
zero(O::NfAbsOrd) = O(fmpz(0))

@doc Markdown.doc"""
    one(O::NfOrd) -> NfAbsOrdElem

Returns the one element of $\mathcal O$.
"""
one(O::NfAbsOrd) = O(fmpz(1))

@doc Markdown.doc"""
    zero(a::NfAbsOrdElem) -> NfAbsOrdElem

Returns the zero element of the parent of $a$.
"""
zero(a::NfAbsOrdElem) = parent(a)(0)

@doc Markdown.doc"""
    one(O::NfOrd) -> NfAbsOrdElem

Returns the one element of the parent of $a$.
"""
one(a::NfAbsOrdElem) = one(parent(a))

################################################################################
#
#  isone/iszero
#
################################################################################

@doc Markdown.doc"""
    isone(a::NfOrd) -> Bool

Tests if $a$ is one.
"""
isone(a::NfAbsOrdElem) = isone(a.elem_in_nf)

@doc Markdown.doc"""
    iszero(a::NfOrd) -> Bool

Tests if $a$ is zero.
"""
iszero(a::NfAbsOrdElem) = iszero(a.elem_in_nf)

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::NfAbsOrdElem)
  print(io, a.elem_in_nf)
end

################################################################################
#
#  Unary operations
#
################################################################################

@doc Markdown.doc"""
    -(x::NfAbsOrdElem) -> NfAbsOrdElem

Returns the additive inverse of $x$.
"""
function -(x::NfAbsOrdElem)
  z = parent(x)()
  z.elem_in_nf = - x.elem_in_nf
  if x.has_coord
    z.coordinates = map(y -> -y, x.coordinates)
    z.has_coord = true
  end
  return z
end

###############################################################################
#
#  Binary operations
#
###############################################################################

@doc Markdown.doc"""
    *(x::NfAbsOrdElem, y::NfAbsOrdElem) -> NfAbsOrdElem

Returns $x \cdot y$.
"""
function *(x::NfAbsOrdElem{S, T}, y::NfAbsOrdElem{S, T}) where {S, T}
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf*y.elem_in_nf
  return z
end

@doc Markdown.doc"""
    +(x::NfAbsOrdElem, y::NfAbsOrdElem) -> NfAbsOrdElem

Returns $x + y$.
"""
function +(x::NfAbsOrdElem, y::NfAbsOrdElem)
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf + y.elem_in_nf
  if x.has_coord && y.has_coord
    z.coordinates =
      [ x.coordinates[i] + y.coordinates[i] for i = 1:degree(parent(x))]
    z.has_coord = true
  end
  return z
end

@doc Markdown.doc"""
    -(x::NfAbsOrdElem, y::NfAbsOrdElem) -> NfAbsOrdElem

Returns $x - y$.
"""
function -(x::NfAbsOrdElem, y::NfAbsOrdElem)
  !check_parent(x, y) && error("Wrong parents")
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf - y.elem_in_nf
  if x.has_coord && y.has_coord
    z.coordinates =
      [ x.coordinates[i] - y.coordinates[i] for i = 1:degree(parent(x))]
    z.has_coord = true
  end
  return z
end

@doc Markdown.doc"""
    divexact(x::NfAbsOrdElem, y::NfAbsOrdElem, check::Bool) -> NfAbsOrdElem

Returns $x/y$. It is assumed that $x/y$ is an element of the same order
as $x$.
"""
function divexact(x::NfAbsOrdElem, y::NfAbsOrdElem, check::Bool = true)
  !check_parent(x, y) && error("Wrong parents")
  a = divexact(x.elem_in_nf, y.elem_in_nf)
  if check
    if !in(a, parent(x))
      error("Quotient not an element of the order")
    end
  end
  z = parent(x)()
  z.elem_in_nf = a
  return z
end

################################################################################
#
#  Ad hoc operations
#
################################################################################

function *(x::NfAbsOrdElem, y::Integer)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf * y
  if x.has_coord
    z.coordinates = map(z -> y*z, x.coordinates)
    z.has_coord = true
  end
  return z
end

*(x::Integer, y::NfAbsOrdElem) = y * x

function *(x::NfAbsOrdElem, y::fmpz)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf * y
  if x.has_coord
    z.coordinates = map(z -> y*z, x.coordinates)
    z.has_coord = true
  end
  return z
end

*(x::fmpz, y::NfAbsOrdElem) = y * x

for T in [Integer, fmpz]
  @eval begin
    function +(x::NfAbsOrdElem, y::$T)
      z = parent(x)()
      z.elem_in_nf = x.elem_in_nf + y
      return z
    end

    +(x::$T, y::NfAbsOrdElem) = y + x

    function -(x::NfAbsOrdElem, y::$T)
      z = parent(x)()
      z.elem_in_nf = x.elem_in_nf - y
      return z
    end

    function -(x::$T, y::NfAbsOrdElem)
      z = parent(y)()
      z.elem_in_nf = x - y.elem_in_nf
      return z
    end
  end
end

################################################################################
#
#  Adhoc exact division
#
################################################################################

for T in [Integer, fmpz]
  @eval begin
    function divexact(x::NfAbsOrdElem, y::$T, check::Bool = true)
      a = divexact(x.elem_in_nf, y)
      if check
        if !in(a, parent(x))
          error("Quotient not an element of the order")
        end
      end
      z = parent(x)()
      z.elem_in_nf = a
      return z
    end
  end
end

################################################################################
#
#  Adhoc division
#
################################################################################

function //(x::nf_elem, y::NfOrdElem)
  check_parent(x, y.elem_in_nf)
  return x//y.elem_in_nf
end

function //(y::NfOrdElem, x::nf_elem)
  check_parent(x, y.elem_in_nf)
  return y.elem_in_nf//x
end

################################################################################
#
#  Exponentiation
#
################################################################################

@doc Markdown.doc"""
    ^(x::NfAbsOrdElem, y::Int)

Returns $x^y$.
"""
function ^(x::NfAbsOrdElem, y::Int)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf^y
  return z
end

################################################################################
#
#  Modular reduction
#
################################################################################

@doc Markdown.doc"""
    mod(a::NfAbsOrdElem, m::Union{fmpz, Int}) -> NfAbsOrdElem

Reduces the coefficient vector of $a$ modulo $m$ and returns the corresponding
element. The coefficient vector of the result will have entries $x$ with
$0 \leq x \leq m$.
"""
function mod(a::NfAbsOrdElem, m::Union{fmpz, Int})
  assure_has_coord(a)
  d = degree(parent(a))
  ar = Vector{fmpz}(undef, d)
  for i in 1:d
    ar[i] = mod(a.coordinates[i], m)
  end
  return NfAbsOrdElem(parent(a), ar) # avoid making a copy of ar
end

################################################################################
#
#  Modular exponentiation
#
################################################################################

@doc Markdown.doc"""
    powermod(a::NfAbsOrdElem, i::fmpz, m::Union{fmpz, Int}) -> NfAbsOrdElem

Returns an element $a^i$ modulo $m$.
"""
function powermod(a::NfAbsOrdElem, i::fmpz, p::fmpz)

  #if contains_equation_order(parent(a))#This doesn't work!
  if issimple(nf(parent(a))) && isdefining_polynomial_nice(nf(parent(a)))
    return powermod_fast(a, i, p)
  else
    return powermod_gen(a, i, p)
  end
end

function powermod_gen(a::NfAbsOrdElem, i::fmpz, p::fmpz)
  if i == 0
    return one(parent(a))
  end
  if i == 1
    b = mod(a,p)
    return b
  end
  if mod(i,2) == 0
    j = div(i, 2)
    b = powermod(a, j, p)
    b = b^2
    b = mod(b,p)
    return b
  end
  b = mod(a * powermod(a, i - 1, p), p)
  return b
end

function powermod_fast(a::NfAbsOrdElem, i::fmpz, p::fmpz)
  if i == 0
    return one(parent(a))
  end
  if i == 1
    bb = mod(a, p)
    return bb
  end

  b = a.elem_in_nf
  d = denominator(b)
  f, e = ppio(d, p) # f = p^k, and e is a unit mod p
  b *= e
  e = invmod(e, p)
  e = powmod(e, i, p)

  y = one(parent(b))
  while i > 1
    if iszero(b)
      return zero(parent(a))
    end
    if iseven(i)
      b = mod(b*b, p)
    else
      y = mod(b*y, p)
      b = mod(b*b, p)
    end
    i = div(i, 2)
  end
  b = mod(b*y, p)

  return mod(parent(a)(b)*e, p)
end

@doc Markdown.doc"""
    powermod(a::NfAbsOrdElem, i::fmpz, m::Union{fmpz, Int}) -> NfAbsOrdElem

Returns an element $a^i$ modulo $m$.
"""
function powermod(a::NfOrdElem, i::fmpz, I::NfOrdIdl)
  if i == 0
    return one(parent(a))
  end
  b = mod(a, I)
  if i == 1
    return b
  end
  if iszero(b)
    return b
  end
  y = one(parent(b))
  while i > 1
    if iseven(i)
      b = mod(b*b, I)
    else
      y = mod(b*y, I)
      b = mod(b*b, I)
    end
    i = div(i, 2)
  end
  b = mod(b*y, I)

  return mod(b, I)
end

denominator(a::NfAbsNSElem) = denominator(a.data)

function mod(a::NfAbsNSElem, p::fmpz)
  b = copy(a)
  for i=1:length(b.data)
    el = coeff(b.data, i)
    dnew, cp = ppio(denominator(el), p)
    el *= cp
    n = mod(numerator(el), dnew * p)
    setcoeff!(b.data, i, fmpq(n, dnew))
  end
  return b
end

@doc Markdown.doc"""
    powermod(a::NfAbsOrdElem, i::Integer, m::Integer) -> NfAbsOrdElem

Returns the element $a^i$ modulo $m$.
"""
powermod(a::NfAbsOrdElem, i::Integer, m::Integer) = powermod(a, fmpz(i), fmpz(m))

@doc Markdown.doc"""
    powermod(a::NfAbsOrdElem, i::fmpz, m::Integer) -> NfAbsOrdElem

Returns the element $a^i$ modulo $m$.
"""
powermod(a::NfAbsOrdElem, i::fmpz, m::Integer)  = powermod(a, i, fmpz(m))

@doc Markdown.doc"""
    powermod(a::NfAbsOrdElem, i::Integer, m::fmpz) -> NfAbsOrdElem

Returns the element $a^i$ modulo $m$.
"""
powermod(a::NfAbsOrdElem, i::Integer, m::fmpz)  = powermod(a, fmpz(i), m)

################################################################################
#
#  Representation matrices
#
################################################################################

@doc Markdown.doc"""
    representation_matrix(a::NfAbsOrdElem) -> fmpz_mat

Returns the representation matrix of the element $a$.
"""
function representation_matrix(a::NfAbsOrdElem)
  O = parent(a)
  assure_has_basis_matrix(O)
  assure_has_basis_mat_inv(O)
  A = representation_matrix(a, nf(O))
  mul!(A, O.basis_matrix, A)
  mul!(A, A, O.basis_mat_inv)
  !(A.den == 1) && error("Element not in order")
  return A.num
end

@doc Markdown.doc"""
    representation_matrix(a::NfAbsOrdElem, K::AnticNumberField) -> FakeFmpqMat

Returns the representation matrix of the element $a$ considered as an element
of the ambient number field $K$. It is assumed that $K$ is the ambient number
field of the order of $a$.
"""
function representation_matrix(a::NfAbsOrdElem{S, T}, K::S) where {S, T}
  nf(parent(a)) != K && error("Element not in this field")
  A, d = Nemo.representation_matrix_q(a.elem_in_nf)
  A.base_ring = FlintZZ
  z = FakeFmpqMat(A, d)
  return z
end

@doc Markdown.doc"""
***
    representation_matrix_mod(a::NfAbsOrdElem, d::fmpz) -> fmpz_mat

Returns the representation matrix of the element $a$ with entries reduced mod d.
"""
function representation_matrix_mod(a::NfAbsOrdElem, d::fmpz)
  O = parent(a)
  assure_has_basis_matrix(O)
  assure_has_basis_mat_inv(O)
  A = representation_matrix(a, nf(O))
  d2 = O.basis_matrix.den * O.basis_mat_inv.den*A.den
  d2c, d2nc = ppio(d2, d)
  d1 = d * d2c
  A1 = A.num
  mod!(A.num, d1)
  M1 = mod(O.basis_matrix.num, d1)
  mul!(A1, M1, A1)
  M2 = mod(O.basis_mat_inv.num, d1)
  mul!(A1, A1, M2)
  mod!(A1, d1)
  divexact!(A1, A1, d2c)
  inver = invmod(d2nc, d1)
  mul!(A1, A1, inver)
  mod!(A1, d)
  return A1
end

################################################################################
#
#  Trace
#
################################################################################

@doc Markdown.doc"""
    tr(a::NfAbsOrdElem) -> fmpz

Returns the trace of $a$.
"""
function tr(a::NfAbsOrdElem)
  return FlintZZ(tr(a.elem_in_nf))
end

################################################################################
#
#  Norm
#
################################################################################

@doc Markdown.doc"""
    norm(a::NfAbsOrdElem) -> fmpz

Returns the norm of $a$.
"""
function norm(a::NfAbsOrdElem)
  return FlintZZ(norm(a.elem_in_nf))
end

################################################################################
#
#  Random element generation
#
################################################################################

# TODO: Make this faster, don't allocate the ar array ...
function rand!(z::NfAbsOrdElem{S, T}, B::Vector{NfAbsOrdElem{S, T}}, R) where {S, T}
  O = parent(z)
  y = O()
  for i in 1:degree(O)
    mul!(y, rand(R), B[i])
    add!(z, z, y)
  end
  return z
end

function rand!(z::NfAbsOrdElem, O::NfAbsOrd, R::UnitRange{T}) where T <: Integer
  y = O()
  ar = rand(R, degree(O))
  B = basis(O, copy = false)
  for i in 1:degree(O)
    mul!(y, rand(R), B[i])
    add!(z, z, y)
  end
  return z
end

@doc Markdown.doc"""
    rand(O::NfOrd, R::UnitRange{Integer}) -> NfAbsOrdElem

Computes a coefficient vector with entries uniformly distributed in `R` and returns
the corresponding element of the order.
"""
function rand(O::NfOrd, R::UnitRange{T}) where T <: Integer
  z = O()
  rand!(z, O, R)
  return z
end

function rand!(z::NfAbsOrdElem, O::NfOrd, n::Union{Integer, fmpz})
  return rand!(z, O, -n:n)
end

@doc Markdown.doc"""
    rand(O::NfOrd, n::Union{Integer, fmpz}) -> NfAbsOrdElem

Computes a coefficient vector with entries uniformly distributed in
$\{-n,\dotsc,-1,0,1,\dotsc,n\}$ and returns the corresponding element of the
order $\mathcal O$.
"""
function rand(O::NfOrd, n::Integer)
  return rand(O, -n:n)
end

function rand(O::NfOrd, n::fmpz)
  z = O()
  rand!(z, O, BigInt(n))
  return z
end

function rand!(z::NfAbsOrdElem, O::NfOrd, n::fmpz)
  return rand!(z, O, BigInt(n))
end

################################################################################
#
#  Unsafe operations
#
################################################################################

@inline function add!(z::NfAbsOrdElem, x::NfAbsOrdElem, y::NfAbsOrdElem)
  add!(z.elem_in_nf, x.elem_in_nf, y.elem_in_nf)
  z.has_coord = false
  return z
end

@inline function sub!(z::NfAbsOrdElem, x::NfAbsOrdElem, y::NfAbsOrdElem)
  sub!(z.elem_in_nf, x.elem_in_nf, y.elem_in_nf)
  z.has_coord = false
  return z
end

@inline function mul!(z::NfAbsOrdElem, x::NfAbsOrdElem, y::NfAbsOrdElem)
  mul!(z.elem_in_nf, x.elem_in_nf, y.elem_in_nf)
  z.has_coord = false
  return z
end

function addeq!(z::NfAbsOrdElem, x::NfAbsOrdElem)
  addeq!(z.elem_in_nf, x.elem_in_nf)
  if x.has_coord && z.has_coord
    for i in 1:degree(parent(z))
      add!(z.coordinates[i], z.coordinates[i], x.coordinates[i])
    end
  end
  return z
end

################################################################################
#
#  Unsafe ad hoc operations
#
################################################################################

# ad hoc
for T in [Integer, fmpz]
  @eval begin
    @inline function mul!(z::NfAbsOrdElem, x::NfAbsOrdElem, y::$T)
      z.elem_in_nf = mul!(z.elem_in_nf, x.elem_in_nf, y)
      z.has_coord = false
      return z
    end

    mul!(z::NfAbsOrdElem, x::$T, y::NfAbsOrdElem) = mul!(z, y, x)
  end
end

for T in [Integer, fmpz]
  @eval begin
    @inline function add!(z::NfAbsOrdElem, x::NfAbsOrdElem, y::$T)
      z.elem_in_nf = add!(z.elem_in_nf, x.elem_in_nf, y)
      z.has_coord = false
      return z
    end

    add!(z::NfAbsOrdElem, x::$T, y::NfAbsOrdElem) = add!(z, y, x)
  end
end

for T in [Integer, fmpz]
  @eval begin
    @inline function sub!(z::NfAbsOrdElem, x::NfAbsOrdElem, y::$T)
      z.elem_in_nf = sub!(z.elem_in_nf, x.elem_in_nf, y)
      z.has_coord = false
      return z
    end

    sub!(z::NfAbsOrdElem, x::$T, y::NfAbsOrdElem) = add!(z, y, x)
  end
end

################################################################################
#
#  Base cases for dot product of vectors
#
################################################################################

dot(x::NfAbsOrdElem, y::Integer) = x * y

dot(x::Integer, y::NfAbsOrdElem) = y * x

dot(x::NfAbsOrdElem, y::fmpz) = x * y

dot(x::fmpz, y::NfAbsOrdElem) = y * x

################################################################################
#
#  Conversion
#
################################################################################

(K::AnticNumberField)(x::NfAbsOrdElem{AnticNumberField, nf_elem}) = elem_in_nf(x)

(K::NfAbsNS)(x::NfAbsOrdElem{NfAbsNS, NfAbsNSElem}) = elem_in_nf(x)

################################################################################
#
#  Minkowski embedding
#
################################################################################

@doc Markdown.doc"""
    minkowski_map(a::NfAbsOrdElem, abs_tol::Int) -> Array{arb, 1}

Returns the image of $a$ under the Minkowski embedding.
Every entry of the array returned is of type `arb` with radius less then
`2^-abs_tol`.
"""
function minkowski_map(a::NfAbsOrdElem, abs_tol::Int = 32)
  # Use a.elem_in_nf instead of elem_in_nf(a) to avoid copying the data.
  # The function minkowski_map does not alter the input!
  return minkowski_map(a.elem_in_nf, abs_tol)
end

################################################################################
#
#  Conjugates
#
################################################################################

@doc Markdown.doc"""
    conjugates_arb(x::NfAbsOrdElem, abs_tol::Int) -> Array{acb, 1}

Compute the the conjugates of `x` as elements of type `acb`.
Recall that we order the complex conjugates
$\sigma_{r+1}(x),...,\sigma_{r+2s}(x)$ such that
$\sigma_{i}(x) = \overline{\sigma_{i + s}(x)}$ for $r + 2 \leq i \leq r + s$.
>
Every entry `y` of the array returned satisfies `radius(real(y)) < 2^-abs_tol`,
`radius(imag(y)) < 2^-abs_tol` respectively.
"""
function conjugates_arb(x::NfAbsOrdElem, abs_tol::Int = 32)
  # Use a.elem_in_nf instead of elem_in_nf(a) to avoid copying the data.
  # The function minkowski_map does not alter the input!
  return conjugates_arb(x.elem_in_nf, abs_tol)
end

@doc Markdown.doc"""
    conjugates_arb_log(x::NfAbsOrdElem, abs_tol::Int) -> Array{arb, 1}

Returns the elements
$(\log(\lvert \sigma_1(x) \rvert),\dotsc,\log(\lvert\sigma_r(x) \rvert),
\dotsc,2\log(\lvert \sigma_{r+1}(x) \rvert),\dotsc,
2\log(\lvert \sigma_{r+s}(x)\rvert))$ as elements of type `arb` radius
less then `2^-abs_tol`.
"""
function conjugates_arb_log(x::NfAbsOrdElem, abs_tol::Int = 32)
  return conjugates_arb_log(x.elem_in_nf, abs_tol)
end

################################################################################
#
#  T2
#
################################################################################

@doc Markdown.doc"""
    t2(x::NfAbsOrdElem, abs_tol::Int = 32) -> arb

Return the $T_2$-norm of $x$. The radius of the result will be less than
`2^-abs_tol`.
"""
function t2(x::NfAbsOrdElem, abs_tol::Int = 32)
  return t2(x.elem_in_nf, abs_tol)
end

################################################################################
#
#  Promotion
#
################################################################################

Nemo.promote_rule(::Type{NfAbsOrdElem{S, T}}, ::Type{U}) where {S, T, U <: Integer} = NfAbsOrdElem{S, T}

Nemo.promote_rule(::Type{NfAbsOrdElem{S, T}}, ::Type{fmpz}) where {S, T} = NfAbsOrdElem{S, T}

Nemo.promote_rule(::Type{NfAbsOrdElem{S, T}}, ::Type{T}) where {S, T} = T

Nemo.promote_rule(::Type{T}, ::Type{NfAbsOrdElem{S, T}}) where {S, T} = T

################################################################################
#
#  Factorization
#
################################################################################

@doc Markdown.doc"""
    factor(a::NfOrdElem) -> Fac{NfOrdElem}

Computes a factorization of `a` into irreducible elements. The return value
is a factorization `fac`, which satisfies `a = unit(fac) * prod(p^e for (p, e)
in fac)`.

The function requires that `a` is non-zero and that all prime ideals containing
`a` are principal, which is for example satisfied if class group of the order
of `a` is trivial.
"""
function factor(a::NfOrdElem)
  iszero(a) && throw(error("Element must be non-zero"))
  OK = parent(a)
  I = a * OK
  D = Dict{NfOrdElem, Int}()
  u = a
  for (p, e) in factor(I)
    b, c = isprincipal(p)
    !b && throw(error("Prime ideal dividing the element not principal"))
    D[c] = e
    u = divexact(u, c^e)
  end
  return Nemo.Fac(u, D)
end


text:
"""
    optimize_seed_set(h::Hypergraph, S::Set{Int}, metaV, metaE; printme::Bool=true)

Remove all unnecessary nodes from the evaluated seed set.
"""
function optimize_seed_set(h::Hypergraph, S::Set{Int}, metaV, metaE; printme::Bool=true)
    seed = deepcopy(S)

    degrees = Dict{Int, Int}([v => length(gethyperedges(h, v)) for v in S])
    s_degrees = sort(collect(degrees), by = x -> x[2])

    ids = [pair[1] for pair in s_degrees]

    for v in ids
        actE = zeros(Bool, nhe(h))
        actV = zeros(Bool, nhv(h))

        # remove v
        delete!(seed, v)

        # activate all nodes in seed\v
        for s in seed
           actV[s] = true
        end

        #simulate diffusion using seed\v
        simres = simulate!(h, actV, actE, metaV, metaE; printme = false)

        # cannot activate all nodes
        if simres.actvs != nhv(h)
            # re-add v to seed
            push!(seed, v)
        end
    end

    printme && println("S $(length(S)) -- new seed set $(length(seed))")

    seed
end


text:
@doc raw"""
    submanifold_component(M::AbstractManifold, p, i::Integer)
    submanifold_component(M::AbstractManifold, p, ::Val(i)) where {i}
    submanifold_component(p, i::Integer)
    submanifold_component(p, ::Val(i)) where {i}

Project the product array `p` on `M` to its `i`th component. A new array is returned.
"""
submanifold_component(::Any...)
@inline function submanifold_component(M::AbstractManifold, p, i::Integer)
    return submanifold_component(M, p, Val(i))
end
@inline submanifold_component(M::AbstractManifold, p, i::Val) = submanifold_component(p, i)
@inline submanifold_component(p, ::Val{I}) where {I} = p.parts[I]
@inline submanifold_component(p::ArrayPartition, ::Val{I}) where {I} = p.x[I]
@inline submanifold_component(p, i::Integer) = submanifold_component(p, Val(i))

@doc raw"""
    submanifold_components(M::AbstractManifold, p)
    submanifold_components(p)

Get the projected components of `p` on the submanifolds of `M`. The components are returned in a Tuple.
"""
submanifold_components(::Any...)
@inline submanifold_components(::AbstractManifold, p) = submanifold_components(p)
@inline submanifold_components(p) = p.parts
@inline submanifold_components(p::ArrayPartition) = p.x

function _show_component(io::IO, v; pre="", head="")
    sx = sprint(show, "text/plain", v, context=io, sizehint=0)
    sx = replace(sx, '\n' => "\n$(pre)")
    return print(io, head, pre, sx)
end

function _show_component_range(io::IO, vs, range; pre="", sym="Component ")
    for i in range
        _show_component(io, vs[i]; pre=pre, head="\n$(sym)$(i) =\n")
    end
    return nothing
end

function _show_product_repr(io::IO, x; name="Product representation", nmax=4)
    n = length(x.parts)
    print(io, "$(name) with $(n) submanifold component$(n == 1 ? "" : "s"):")
    half_nmax = div(nmax, 2)
    pre = "  "
    sym = " Component "
    if n ≤ nmax
        _show_component_range(io, x.parts, 1:n; pre=pre, sym=sym)
    else
        _show_component_range(io, x.parts, 1:half_nmax; pre=pre, sym=sym)
        print(io, "\n ⋮")
        _show_component_range(io, x.parts, (n - half_nmax + 1):n; pre=pre, sym=sym)
    end
    return nothing
end

"""
    ProductRepr(parts)

A more general but slower representation of points and tangent vectors on
a product manifold.

# Example:

A product point on a product manifold `Sphere(2) × Euclidean(2)` might be
created as

    ProductRepr([1.0, 0.0, 0.0], [2.0, 3.0])

where `[1.0, 0.0, 0.0]` is the part corresponding to the sphere factor
and `[2.0, 3.0]` is the part corresponding to the euclidean manifold.
"""
struct ProductRepr{TM<:Tuple}
    parts::TM
end

ProductRepr(points...) = ProductRepr{typeof(points)}(points)

Base.:(==)(x::ProductRepr, y::ProductRepr) = x.parts == y.parts

@inline function number_eltype(x::ProductRepr)
    @inline eti_to_one(eti) = one(number_eltype(eti))
    return typeof(sum(map(eti_to_one, x.parts)))
end

allocate(x::ProductRepr) = ProductRepr(map(allocate, submanifold_components(x))...)
function allocate(x::ProductRepr, ::Type{T}) where {T}
    return ProductRepr(map(t -> allocate(t, T), submanifold_components(x))...)
end
allocate(p::ProductRepr, ::Type{T}, s::Size{S}) where {S,T} = Vector{T}(undef, S)
allocate(p::ProductRepr, ::Type{T}, s::Integer) where {S,T} = Vector{T}(undef, s)

Base.copy(x::ProductRepr) = ProductRepr(map(copy, x.parts))

function Base.copyto!(x::ProductRepr, y::ProductRepr)
    map(copyto!, submanifold_components(x), submanifold_components(y))
    return x
end

function Base.:+(v1::ProductRepr, v2::ProductRepr)
    return ProductRepr(map(+, submanifold_components(v1), submanifold_components(v2))...)
end

function Base.:-(v1::ProductRepr, v2::ProductRepr)
    return ProductRepr(map(-, submanifold_components(v1), submanifold_components(v2))...)
end
Base.:-(v::ProductRepr) = ProductRepr(map(-, submanifold_components(v)))

Base.:*(a::Number, v::ProductRepr) = ProductRepr(map(t -> a * t, submanifold_components(v)))
Base.:*(v::ProductRepr, a::Number) = ProductRepr(map(t -> t * a, submanifold_components(v)))

Base.:/(v::ProductRepr, a::Number) = ProductRepr(map(t -> t / a, submanifold_components(v)))

function Base.convert(::Type{ProductRepr{TPR}}, x::ProductRepr) where {TPR<:Tuple}
    if @isdefined TPR
        return ProductRepr(convert(TPR, submanifold_components(x)))
    else
        return x
    end
end

function Base.show(io::IO, ::MIME"text/plain", x::ProductRepr)
    return _show_product_repr(io, x; name="ProductRepr")
end

ManifoldsBase._get_vector_cache_broadcast(::ProductRepr) = Val(false)

# Tuple-like broadcasting of ProductRepr

function Broadcast.BroadcastStyle(::Type{<:ProductRepr})
    return Broadcast.Style{ProductRepr}()
end
function Broadcast.BroadcastStyle(
    ::Broadcast.AbstractArrayStyle{0},
    b::Broadcast.Style{ProductRepr},
)
    return b
end

Broadcast.instantiate(bc::Broadcast.Broadcasted{Broadcast.Style{ProductRepr},Nothing}) = bc
function Broadcast.instantiate(bc::Broadcast.Broadcasted{Broadcast.Style{ProductRepr}})
    Broadcast.check_broadcast_axes(bc.axes, bc.args...)
    return bc
end

Broadcast.broadcastable(v::ProductRepr) = v

@inline function Base.copy(bc::Broadcast.Broadcasted{Broadcast.Style{ProductRepr}})
    dim = axes(bc)
    length(dim) == 1 || throw(DimensionMismatch("ProductRepr only supports one dimension"))
    N = length(dim[1])
    return ProductRepr(ntuple(k -> @inbounds(Broadcast._broadcast_getindex(bc, k)), Val(N)))
end

Base.@propagate_inbounds Broadcast._broadcast_getindex(v::ProductRepr, I) = v.parts[I[1]]

Base.axes(v::ProductRepr) = axes(v.parts)

@inline function Base.copyto!(
    dest::ProductRepr,
    bc::Broadcast.Broadcasted{Broadcast.Style{ProductRepr}},
)
    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
    # Performance optimization: broadcast!(identity, dest, A) is equivalent to copyto!(dest, A) if indices match
    if bc.f === identity && bc.args isa Tuple{ProductRepr} # only a single input argument to broadcast!
        A = bc.args[1]
        if axes(dest) == axes(A)
            return copyto!(dest, A)
        end
    end
    bc′ = Broadcast.preprocess(dest, bc)
    # Performance may vary depending on whether `@inbounds` is placed outside the
    # for loop or not. (cf. https://github.com/JuliaLang/julia/issues/38086)
    @inbounds @simd for I in eachindex(bc′)
        copyto!(dest.parts[I], bc′[I])
    end
    return dest
end

## ArrayPartition

ManifoldsBase._get_vector_cache_broadcast(::ArrayPartition) = Val(false)


text:
"""
`$(current_module())` is a module containing all of the code regarding conduction of experiments.
"""
module Experiments

using MicroLogging
using LaTeXStrings
using Plots
import PyPlot # workaround due to https://github.com/JuliaPlots/Plots.jl/issues/1047
using DataFrames
using DataArrays
using CSV

importall RobRecSolver

export
    # experiments
    runExperiments,
    exportKnapsackResults,
    exportAssignmentResults,
    getProperties,

    # data generators
    generateData

    files = [
            "data_generators",
            "experiment",
            "export" # Add more files here
        ]

    for file in files
        include("$(file).jl")
    end

end