text:
struct Point
    x::Float64
    y::Float64
end
Point(vec) = Point(vec[1], vec[2])

import Base.show ; Base.show(io::IO, p::Point) = print(io, "P($(round(p.x, 2)), $(round(p.y, 2)))")
import Base.first ; first(p::Point) = p.x
import Base.last ; last(p::Point) = p.y

struct Segment{T}
    p1::Point
    p2::Point
    c::Point
    a::Float64
	b::Float64
	Segment{T}(p1::P, p2::P, c::P, a, b) where {T<:Sense, P<:Point} = new(p1, p2, c, a, b)
	Segment{T}(p1, p2, c) where T = Segment{T}(Point(p1), Point(p2), Point(c))
	function Segment{T}(p1::Point, p2::Point, c::Point) where T
		a = (p2.y - p1.y) / (p2.x - p1.x)
		b = p1.y - a*p1.x
		Segment{T}(p1, p2, c, a, b)
	end
end
Base.show(io::IO, s::Segment) = print(io, "Sgm[$(s.p1), $(s.p2)]")

nadir(v::Vector{Segment{Min}}) = (maximum(x-> x.c.x, v), maximum(x-> x.c.y, v))
nadir(v::Vector{Segment{Max}}) = (minimum(x-> x.c.x, v), minimum(x-> x.c.y, v))

function filterSegment(s::Segment{Min}, u0::Point)
	u = Point(u0.x - 0.5, u0.y - 0.5)
	S = Segment{Min}[]
	p1, p2, c = s.p1, s.p2, s.c
	if u.y >= s.a * u.x + s.b || u.y >= p1.y || u.x >= p2.x #u does not dominate any point on segment s
		if u.x <= p1.x #u is above and to the left of p1
			push!(S, Segment{Min}(p1, p2, Point(c.x, min(c.y, u.y))))
		elseif u.y <= p2.y #u is below and to the right of p2
			push!(S, Segment{Min}(p1, p2, Point(min(c.x, u.x), c.y)))
		else
			push!(S, s)
		end
	else #at least one part of s is dominated by u
		if u.x > p1.x #u does not dominate p1
			push!(S, Segment{Min}(p1, Point(u.x,s.a * u.x + s.b), Point(u.x,c.y)))
		end
		if u.y > p2.y #u does not dominate p2
			push!(S, Segment{Min}(Point((u.y - s.b) / s.a, u.y), p2,Point(c.x, u.y)))
		end
	end
	return S
end

function filterSegment(s::Segment{Max}, u0::Point)
	u = Point(u0.x + 0.5, u0.y + 0.5)
	S = Segment{Max}[]
	p1, p2, c = s.p1, s.p2, s.c
	if u.y <= s.a * u.x + s.b || u.y <= p2.y || u.x <= p1.x #u does not dominate any point on segment s
		if u.y >= p1.y #u is above and to the left of p1
			push!(S, Segment{Max}(p1, p2, Point(max(c.x, u.x), c.y)))
		elseif u.x >= p2.x #u is below and to the right of p2
			push!(S, Segment{Max}(p1, p2, Point(c.x, max(u.y, c.y))))
		else
			push!(S, s)
		end
	else #at least one part of s is dominated by u
		if u.y < p1.y #u does not dominate p1
			push!(S, Segment{Max}(p1, Point((u.y - s.b) / s.a, u.y), Point(c.x,u.y)))
		end
		if u.x < p2.x #u does not dominate p2
			push!(S, Segment{Max}(Point(u.x, s.a * u.x + s.b), p2, Point(u.x, c.y)))
		end
	end
	return S
end

function coversInteger(s::Segment{Min})
	p1, p2, a, b, c = s.p1, s.p2, s.a, s.b, s.c
	return (p1 == p2 || floor(c.y) >= a * floor(c.x) + b) && floor(c.y) >= ceil(p2.y) && floor(c.x) >= ceil(p1.x)
end

function coversInteger(s::Segment{Max})
	p1, p2, a, b, c = s.p1, s.p2, s.a, s.b, s.c
	return (p1 == p2 || ceil(c.y) <= a * ceil(c.x) + b) && ceil(c.y) <= floor(p1.y) && ceil(c.x) <= ceil(p2.x)
end

function filterDualBound(dual::Vector{Segment{T}}, primal, showplot=false) where T<:Sense
	S = dual
	for u in primal
		S_prime = Segment{T}[]
		for s in S
			append!(S_prime, filterSegment(s, u))
		end
		S = S_prime
	end
	res = groupby_continuous(S)

	#Plotting
	if showplot
		clf()
		plotdualbound(S)
		plot(first.(primal), last.(primal), "ko", markersize=4)
		nadirsline_x=[primal[1].x] ; nadirsline_y=[primal[1].y]
		for i = 2:length(primal)
			if T==Min
				append!(nadirsline_x, [primal[i].x, primal[i].x]) ; append!(nadirsline_y, [primal[i-1].y, primal[i].y])
			else
				append!(nadirsline_x, [primal[i-1].x, primal[i].x]) ; append!(nadirsline_y, [primal[i].y, primal[i].y])
			end
		end
		plot(nadirsline_x, nadirsline_y, "k-", linewidth=1) ; plot(nadirsline_x, nadirsline_y, "ks", markersize=2)
		sleep(0.2)
	end

	for r in res
		filter!(coversInteger, r)
	end
	filter!(!isempty, res)
	res
end

function toSegmentList(lb, nadir, sense)
	res = Segment{sense}[]
	if length(lb) == 1
		push!(res, Segment{sense}(lb[1], lb[1], nadir))
	else
		for i = 1:length(lb)-1
			c = if sense==Min
				Point(i==length(lb)-1 ? nadir[1] : lb[i+1][1], nadir[2])
			else
				Point(nadir[1], i==length(lb)-1 ? nadir[2] : lb[i+1][2])
			end
			push!(res, Segment{sense}(Point(lb[i]...), Point(lb[i+1]...), c))
		end
	end
	res
end

toPointList(ub) = [Point(u...) for u in ub]

function groupby_continuous(v)
	isempty(v) && return Vector{Segment}[]
	res = [[v[1]]]
	for i = 2:length(v)
		if v[i].p1 == res[end][end].p2
			push!(res[end], v[i])
		else
			push!(res, [v[i]])
		end
	end
	res
end

function apply_cuts!(n::NodeParragh, cstrData, lift_covers, nb_try = 15)
	for i = 1:nb_try
		res = @suppress solve(n.m, method=:dicho, relax=true)
		if res != :Optimal
			return #peut arriver si le triangle ne contient aucune autre solution que ses deux points extremes, qui sont élminés par les deux contraintes sur les objectifs(+1)
		end
		n.x = [[round(getvalue(JuMP.Variable(n.m, i), j), 8) for i = 1:n.m.numCols] for j = 1:length(getY_N(n.m))]
		success = find_cover_cuts(n, cstrData, lift_covers)
		!success && return
	end
	return
end

function solve_parragh(model ;  showplot = false, docovercuts = true, global_branch = false, use_nsga = true, global_nsga = true, lift_covers = true, time_limit=300, fill_triangles=true, preprocess=true)

	tic()
	vm = copy(model)
	vd = getvOptData(vm)
	@assert length(vd.objs) == 2
	@assert all(isonlybinary.(vm, vd.objs))
	@assert vd.objSenses[1] == vd.objSenses[2]

	sense = vd.objSenses[1] == :Max ? Max : Min
	JuMP.setobjectivesense(vm, vd.objSenses[1])
	cstrData = ConstraintData(vm)

	z1, z2 = vd.objs
	status = global_branch ? solve(vm, method=:lexico, verbose=false) : solve(vm, method=:dicho, round_results=true)
	status == :Optimal || return status

	YN_convex = map(x->round.(Int, x), getY_N(vm))
	XE_convex = [[round(getvalue(JuMP.Variable(vm, i), j)) for i = 1:vm.numCols] for j = 1:length(YN_convex)]

	if length(YN_convex) == 1 || YN_convex[1] == YN_convex[end]
		return @NT(YN = YN_convex, XE = XE_convex, nodes=0)
	end

	if !issorted(YN_convex, by = first)
		s = sortperm(YN_convex, by = first)
        YN_convex = YN_convex[s]
        XE_convex = XE_convex[s]
	end

	if use_nsga
		modelnsga = copy(vm)
		modelnsga.ext[:vOpt].objs = [z1, z2]

		if global_nsga
			ns = nsga_binary(50, 1000, modelnsga, seed=XE_convex, pmut=0.3, showprogress=false)
		else
			ns = union((nsga_binary(20, 500, modelnsga, seed=[XE_convex[i], XE_convex[i+1]], pmut=0.3, showprogress=false) for i = 1:length(XE_convex)-1)...)
			#if we do one nsga per triangle, we can have dominated solutions here.
			NSGAII.fast_non_dominated_sort!(ns, sense==Max ? NSGAII.Max() : NSGAII.Min())
			filter!(x->x.rank == 1, ns)
		end

		ns = unique(x->x.y, ns)
		sort!(ns, by=x->x.y[1])
		LN_NSGA = NonDomPoints(sense, map(x->x.pheno, ns), map(x->Tuple(x.y), ns))

		if fill_triangles && use_nsga
			diffs = [YN_convex[i] .- YN_convex[i+1] for i = 1:length(YN_convex)-1]
			areas = abs.(map(prod, diffs))
			while maximum(areas) > 2*mean(areas)
				i = indmax(areas)
				list_to_add = filter!(x-> YN_convex[i][1] < x.y[1] < YN_convex[i+1][1], ns)
				isempty(list_to_add) && break
				sol_to_add = list_to_add[length(list_to_add)÷2 + 1]
				insert!(YN_convex, i+1, sol_to_add.y)
				insert!(XE_convex, i+1, sol_to_add.pheno)
				diffs = [YN_convex[i] .- YN_convex[i+1] for i = 1:length(YN_convex)-1]
				areas = abs.(map(prod, diffs))
			end
		end
	end

	t = toq()

	nbNodesTotal = 0
	LNGlobal = use_nsga ? LN_NSGA :  NonDomPoints(sense, XE_convex, Tuple.(YN_convex))

	for i = 1:length(YN_convex)-1

		t > time_limit && break
		tic()

		#println("exploring triangle ◬ ($(YN_convex[i]) - $(YN_convex[i+1])) area : $(areas[i])")

		m = copy(vm)
		LN = NonDomPoints(sense, [XE_convex[i], XE_convex[i+1]], [Tuple(YN_convex[i]), Tuple(YN_convex[i+1])])

		if use_nsga
			for j = 1:length(LN_NSGA.yn)
				if LN.yn[1][1] < LN_NSGA.yn[j][1] < LN.yn[end][1]
					push!(LN, LN_NSGA.xe[j], LN_NSGA.yn[j])
				end
			end
		end

		Ƶ1, Ƶ2 = copy(z1, m), copy(z2, m)

		if sense == Max
			bound1 = LN.yn[1][1] + 1
			bound2 = LN.yn[end][2] + 1
			cstr1 = @constraint(m, Ƶ1.aff >= bound1)
			cstr2 = @constraint(m, Ƶ2.aff >= bound2)
			# @constraint(m, Ƶ.aff <= LN.λ[1]*LN.yn[1][1] + LN.λ[2]*LN.yn[1][2])
		else
			bound1 = LN.yn[end][1] - 1
			bound2 = LN.yn[1][2] - 1
			cstr1 = @constraint(m, Ƶ1.aff <= bound1)
			cstr2 = @constraint(m, Ƶ2.aff <= bound2)
		end

		#Stack of nodes to evaluate
		S = [NodeParragh(m, bound1, bound2, cstr1, cstr2)]

		if docovercuts
			apply_cuts!(S[1], cstrData, lift_covers)
		end

		if preprocess
			Ƶ = LN.λ[1]*Ƶ1 + LN.λ[2]*Ƶ2
			JuMP.setobjective(m, vd.objSenses[1], Ƶ.aff)
			preprocess!(S[1], LN, Ƶ1, Ƶ2)
		end

		t += toq()
		#Solve while there are nodes to process
		cpt = 0
		while !isempty(S) && t < time_limit
			tic()
			sort!(S, by = x->x.z_avg, rev=(sense==Min))
			process_node_parragh(pop!(S), S, sense, LN, Ƶ1, Ƶ2, LNGlobal, cstrData, showplot)
			nbNodesTotal += 1
			cpt += 1
			t += toq()
		end

		for i = 1:length(LN.xe)
			push!(LNGlobal, LN.xe[i], LN.yn[i])
		end

	end

	return @NT(YN = LNGlobal.yn	, XE = LNGlobal.xe, nodes = nbNodesTotal, timeout = t > time_limit)
end

function process_node_parragh(n::NodeParragh, S, sense, LN, obj1, obj2, LNGlobal, cstrData, showplot)

	fix_variables!(n)
	setRHS!(n)
	res = @suppress solve(n.m, method=:dicho, relax=true)
	unfix_variables!(n)
	res != :Optimal && return

	YN_relax = getY_N(n.m)
	n.z_avg = mean(x -> LN.λ[1]*x[1] + LN.λ[2]*x[2], YN_relax)
	n.x = [[round(getvalue(JuMP.Variable(n.m, i), j), 8) for i = 1:n.m.numCols] for j = 1:length(YN_relax)]

	new_int_found = false
	inds = find(x -> all(isinteger, x), n.x)
	for i in inds
		z1, z2 = YN_relax[i]
		if !isweaklydominated(z1, z2, LN)
			push!(LN, n.x[i], (round(z1), round(z2)))
			new_int_found = true
		end
	end

	dual = toSegmentList(YN_relax, local_nadir(LN), sense);
	primal = toPointList(LN.yn)

	if showplot
		clf() ;
		for s in dual
			plot([s.p1.x, s.p2.x], [s.p1.y, s.p2.y], "g-", linewidth=1)
			plot([s.p1.x, s.p2.x], [s.p1.y, s.p2.y], "gs", markersize=4)
		end
		plot(first.(primal), last.(primal), "ko", markersize=4)
		nadirsline_x=[primal[1].x]
		nadirsline_y=[primal[1].y]
		for i = 2:length(primal)
			if sense==Min
				append!(nadirsline_x, [primal[i].x, primal[i].x])
				append!(nadirsline_y, [primal[i-1].y, primal[i].y])
			else
				append!(nadirsline_x, [primal[i-1].x, primal[i].x])
				append!(nadirsline_y, [primal[i].y, primal[i].y])
			end
		end
		plot(nadirsline_x, nadirsline_y, "k-", linewidth=1)
		plot(nadirsline_x, nadirsline_y, "ks", markersize=2)
		sleep(0.01)
	end

	filteredDual = MOBC.filterDualBound(dual, primal, showplot)

	if isempty(filteredDual)
		# println("dual bound set empty, discarded")
		#Dual bound set is dominated by the primal set, discard that node
		return
	elseif length(filteredDual) == 1
		if new_int_found
			# println("int found -> pareto-branching")
			push!(S, parragh_branch(n, first(filteredDual), obj1, obj2))
		else
			#branch classique
			# println("branch classique")
			n1, n2 = basicbranch_parragh(n, filteredDual[1])
			push!(S, n1)
			push!(S, n2)
		end
	else
		# println("pareto branching")
		for segmentList in filteredDual
			push!(S, parragh_branch(n, segmentList, obj1, obj2))
		end
	end
	return
end

function parragh_branch(n, segmentList::Vector{Segment{T}}, obj1, obj2) where T
	node = shallowcopy(n)
	boundz1, boundz2 = nadir(segmentList)

	if T == Max
		node.bound1 = max(node.bound1, boundz1)
		node.bound2 = max(node.bound2, boundz2)
	else
		node.bound1 = min(node.bound1, boundz1)
		node.bound2 = min(node.bound2, boundz2)
	end
	node
end

function basicbranch_parragh(n, dualbound)
	frac = zeros(length(n.x[1]))
	for i = 1:length(n.x)
		for j = 1:length(n.x[i])
			if !isinteger(n.x[i][j])
				frac[j] += min(n.x[i][j], 1. - n.x[i][j])
			end
		end
	end

	ind = if maximum(frac) > 1E-6
		indmax(frac)
	else
		findfirst(x->!(x in n.f0) && !(x in n.f1), 1:n.m.numCols)
	end

	n1, n2 = shallowcopy(n), shallowcopy(n)

	push!(n1.f0, ind)
	push!(n2.f1, ind)

	# println("branching on variable $ind")
	n1, n2
end

function bench(m)
	println("global branch")
	println("\tParragh")
	@time a = solve_parragh(m, Inf, use_nsga=true, global_branch=true)
	println("\tStidsen")
	@time b = solve_stidsen(m, Inf, use_nsga=true, global_branch=true)
	println("\nlocal branch")
	println("\tParragh")
	@time c = solve_parragh(m, Inf, use_nsga=true, global_branch=false)
	println("\tStidsen")
	@time d = solve_stidsen(m, Inf, use_nsga=true, global_branch=false)
	@assert a.YN == b.YN == c.YN == d.YN
	nothing
end


text:
using .SPARK
using .VPRK

#*****************************************************************************#
# General initialization functions for all integrators                        #
#*****************************************************************************#

"""
```julia
Integrator(equation, tableau, Δt)
```

Universal constructor for Runge-Kutta (standard, partitioned, variational),
SPARK and splitting integrators that automatically selects the appropriate
integrator based on the equation and tableau types.

"""
function Integrator end

# Print error for integrators not implemented, yet.
function Integrator(equation::Equation, tableau::Union{AbstractTableau,Tableau}, Δt)
    error("No integrator found for equation ", equation, " and tableau ", tableau)
end

#*****************************************************************************#
# Initialization functions for deterministic integrators                      #
#*****************************************************************************#

# Create integrator for Runge-Kutta tableau.
function Integrator(equation::ODE, tableau::Tableau, Δt)
    if isexplicit(tableau)
        # Create integrator for explicit Runge-Kutta tableau
        IntegratorERK(equation, tableau, Δt)
    elseif isdiagnonallyimplicit(tableau)
        # Create integrator for diagonally implicit Runge-Kutta tableau
        IntegratorDIRK(equation, tableau, Δt)
    elseif isfullyimplicit(tableau)
        # Create integrator for fully implicit Runge-Kutta tableau
        IntegratorFIRK(equation, tableau, Δt)
    end
end

# Create integrator for explicit partitioned Runge-Kutta tableau.
function Integrator(equation::Union{PODE,HODE}, tableau::PartitionedTableau, Δt)
    if isexplicit(tableau)
        IntegratorEPRK(equation, tableau, Δt)
    else
        IntegratorIPRK(equation, tableau, Δt)
    end
end

function Integrator(equation::Union{PODE,HODE}, tableau::Tableau, Δt)
    Integrator(equation, PartitionedTableau(tableau), Δt)
end

# Create integrator for implicit partitioned Runge-Kutta tableau.
function Integrator(equation::Union{IODE,LODE}, tableau::PartitionedTableau, Δt)
    IntegratorPRKimplicit(equation, tableau, Δt)
end

# Create integrator for variational partitioned Runge-Kutta tableau.
function Integrator(equation::Union{IODE,LODE}, tableau::TableauVPRK, Δt)
    IntegratorVPRK(equation, tableau, Δt)
end

# Create integrator for formal Lagrangian Runge-Kutta tableau.
function Integrator(equation::LODE, tableau::Tableau, Δt)
    IntegratorFLRK(equation, tableau, Δt)
end

# Create integrator for Projected Gauss-Legendre Runge-Kutta tableau.
function Integrator(equation::Union{IODE,LODE}, tableau::CoefficientsPGLRK, Δt)
    IntegratorPGLRK(equation, tableau, Δt)
end

# Create integrator for variational partitioned additive Runge-Kutta tableau.
function Integrator(equation::Union{IDAE,LDAE}, tableau::TableauVPARK, Δt)
    IntegratorVPARK(equation, tableau, Δt)
end

# Create integrator for special partitioned additive Runge-Kutta tableau.
function Integrator(equation::Union{IDAE,LDAE}, tableau::TableauSPARK, Δt)
    IntegratorSPARK(equation, tableau, Δt)
end

# Create integrator for variational special partitioned additive Runge-Kutta tableau.
function Integrator(equation::Union{IDAE,LDAE}, tableau::TableauVSPARK, Δt)
    IntegratorVSPARK(equation, tableau, Δt)
end

# Create integrator for variational special partitioned additive Runge-Kutta tableau with projection on primary constraint.
function Integrator(equation::Union{IDAE,LDAE}, tableau::TableauVSPARKprimary, Δt)
    IntegratorVSPARKprimary(equation, tableau, Δt)
end

# Create integrator for variational special partitioned additive Runge-Kutta tableau with projection on secondary constraint.
function Integrator(equation::LDAE, tableau::TableauVSPARKsecondary, Δt)
    IntegratorVSPARKsecondary(equation, tableau, Δt)
end

# Create integrator for Hamiltonian partitioned additive Runge-Kutta tableau.
function Integrator(equation::Union{PDAE,HDAE}, tableau::TableauHPARK, Δt)
    IntegratorHPARK(equation, tableau, Δt)
end

# Create integrator for Hamiltonian special partitioned additive Runge-Kutta tableau.
function Integrator(equation::Union{PDAE,HDAE}, tableau::TableauHSPARK, Δt)
    IntegratorHSPARK(equation, tableau, Δt)
end

# Create integrator for Hamiltonian special partitioned additive Runge-Kutta tableau with projection on primary constraint.
function Integrator(equation::Union{PDAE,HDAE}, tableau::TableauHSPARKprimary, Δt)
    IntegratorHSPARKprimary(equation, tableau, Δt)
end

# Create integrator for splitting tableau.
function Integrator(equation::SODE, tableau::AbstractTableauSplitting, Δt)
    IntegratorSplitting(equation, tableau, Δt)
end

#*****************************************************************************#
# Constructor wrappers for deterministic integrators                          #
#*****************************************************************************#

"""
```julia
IntegratorConstructor(DT, D)
IntegratorConstructor(DT, D, tableau)
```

Returns a closure for various integrator constructors that is needed for
composition methods based on the tableau type. If not tableau is provided,
a closure for the constructor of an exact solution is returned.

- `DT`: data type of the solution / state vector
- `D`: dimension of the state vector

"""
function IntegratorConstructor end

# Create integrator constructor for exact solution.
function IntegratorConstructor(DT, D)
    (v::Function, Δt::Number; kwargs...) -> IntegratorExactODE{DT,D}(v, Δt; kwargs...)
end

# Create integrator constructor for Runge-Kutta tableau.
function IntegratorConstructor(DT, D, tableau::Tableau)
    if isexplicit(tableau)
        # Create integrator constructor for explicit Runge-Kutta tableau
        (v::Function, Δt::Number; kwargs...) -> IntegratorERK{DT,D}(v, tableau, Δt; kwargs...)
    elseif isdiagnonallyimplicit(tableau)
        # Create integrator constructor for diagonally implicit Runge-Kutta tableau
        (v::Function, Δt::Number; kwargs...) -> IntegratorDIRK{DT,D}(v, tableau, Δt; kwargs...)
    elseif isfullyimplicit(tableau)
        # Create integrator constructor for fully implicit Runge-Kutta tableau
        (v::Function, Δt::Number; kwargs...) -> IntegratorFIRK{DT,D}(v, tableau, Δt; kwargs...)
    end
end

#*****************************************************************************#
# General integration functions for all integrators                           #
#*****************************************************************************#

"""
```julia
integrate(equation, integrator, ntime; kwargs...)
integrate(equation, tableau, Δt, ntime; kwargs...)
```

Integrate an `equation` with `integrator` for `ntime` time steps
and return the solution.
If a `tableau` and a time step `Δt` is passed instead of an
`integrator`, the appropriate integrator is created automatically.

Some convenience methods exist for the integration of ODEs,
```julia
integrate(v, x₀, tableau, Δt, ntime; t₀=0., kwargs...)
```
where `v` is the function for the vector field, `x₀` the initial condition
and `t₀` the initial time, and for PODEs
```julia
integrate(v, f, q₀, p₀, tableau, Δt, ntime; t₀=0., kwargs...)
```
with vector fields `v` and `f` and initial conditions `q₀` and `p₀`.
"""
function integrate end

# Apply integrator for ntime time steps and return solution.
function integrate(equation::Equation, integrator::Integrator, ntime::Int; kwargs...)
    solution = Solution(equation, timestep(integrator), ntime; kwargs...)
    integrate!(integrator, solution)
    return solution
end

# Integrate given equation with given tableau for ntime time steps and return solution.
function integrate(equation::Equation, tableau::Union{AbstractTableau,Tableau,PartitionedTableau}, Δt, ntime; kwargs...)
    return integrate(equation, Integrator(equation, tableau, Δt), ntime; kwargs...)
end

# Integrate ODE specified by vector field and initial condition with given tableau for ntime time steps and return solution.
function integrate(f::Function, x₀::Vector, tableau::Union{AbstractTableau,Tableau}, Δt, ntime; t₀=0., kwargs...)
    return integrate(ODE(f, t₀, x₀), tableau, Δt, ntime; kwargs...)
end

# Integrate PODE specified by two vector fields and initial conditions with given tableau for ntime time steps and return solution.
function integrate(v::Function, f::Function, q₀::Vector, p₀::Vector, tableau::Union{AbstractTableau,Tableau,PartitionedTableau}, Δt, ntime; t₀=0., kwargs...)
    return integrate(PODE(v, f, t₀, q₀, p₀), tableau, Δt, ntime; kwargs...)
end

#*****************************************************************************#
# Integration functions for deterministic integrators                         #
#*****************************************************************************#

"""
```julia
integrate!(integrator, solution)
```

Perform the actual integration for a given `integrator` and `solution`.
Solve all time steps for all initial conditions.

```julia
integrate!(integrator, solution, m₁, m₂)
```

Solve all time steps for initial conditions m with m₁ ≤ m ≤ m₂.

```julia
integrate!(integrator, solution, m₁, m₂, n₁, n₂)
```

Solve time steps n with n₁ ≤ n ≤ n₂ for initial conditions m with m₁ ≤ m ≤ m₂.

```julia
integrate!(integrator, solution, atomic_solution, m, n)
```

Solve one time step n for one initial condition m.

"""
function integrate! end

# Parts of one integration step that are common to deterministic and stochastic equations.
function integrate_common!(int::Integrator{DT,TT}, sol::Solution{AT,TT}, asol::AtomicSolution{DT,TT}, m::Int, n::Int) where {DT, TT, AT <: AbstractArray{DT}}
    # integrate one initial condition for one time step
    integrate_step!(int, asol)

    # take care of periodic solutions
    cut_periodic_solution!(asol, periodicity(sol))

    # copy solution from cache to solution
    set_solution!(sol, asol, n, m)
end

# Integrate one time step n for one initial condition m.
function integrate!(int::Integrator{DT,TT}, sol::Solution{AT,TT}, asol::AtomicSolution{DT,TT}, m::Int, n::Int) where {DT, TT, AT <: AbstractArray{DT}}
    integrate_common!(int, sol, asol, m, n)
end

# Integrate `equation` for all initial conditions and all time steps.
function integrate!(int::Integrator, sol::Solution)
    integrate!(int, sol, 1, nsamples(sol))
end

# Integrate equation for initial conditions m with m₁ ≤ m ≤ m₂.
function integrate!(int::Integrator, sol::Solution, m₁, m₂)
    # integrate samples m with m₁ ≤ m ≤ m₂ for all time steps
    integrate!(int, sol, m₁, m₂, 1, ntime(sol))
end

# Integrate equation for initial conditions m with m₁ ≤ m ≤ m₂ for time steps n with n₁ ≤ n ≤ n₂.
function integrate!(int::Integrator{DT,TT}, sol::Solution{AT,TT}, m₁::Int, m₂::Int, n₁::Int, n₂::Int) where {DT, TT, AT <: AbstractArray{DT}}
    @assert m₁ ≥ 1
    @assert m₂ ≥ m₁
    @assert m₂ ≤ nsamples(sol)

    @assert n₁ ≥ 1
    @assert n₂ ≥ n₁
    @assert n₂ ≤ ntime(sol)

    asol = AtomicSolution(sol, int)

    # loop over initial conditions showing progress bar
    for m in m₁:m₂
        # get cache from solution
        get_initial_conditions!(sol, asol, m, n₁)
        initialize!(int, asol)

        # loop over time steps
        for n in n₁:n₂
            integrate!(int, sol, asol, m, n)

            # try
            #     integrate!(int, sol, asol, m, n)
            # catch ex
            #     tstr = " in time step " * string(n)
            #
            #     if m₁ ≠ m₂
            #         tstr *= " for initial condition " * string(m)
            #     end
            #
            #     tstr *= "."
            #
            #     if isa(ex, DomainError)
            #         @warn("Domain error" * tstr)
            #     elseif isa(ex, ErrorException)
            #         @warn("Simulation exited early" * tstr)
            #         @warn(ex.msg)
            #     else
            #         @warn(string(typeof(ex)) * tstr)
            #         throw(ex)
            #     end
            # end
        end
    end
end

ext: jl


text:
# A palindromic number reads the same both ways.
# The largest palindrome made from the product
# of two 2-digit numbers is 9009 = 91 × 99.
# Find the largest palindrome made from the product of two 3-digit numbers.

function is_palindrome(n)
    s = string(n)
    return s == reverse(s)
end

@assert is_palindrome(2)
@assert is_palindrome(11)
@assert is_palindrome(9009)
@assert is_palindrome(90109)
@assert is_palindrome(90209)

function find_next_palindrome_decreasing(n)
    while !is_palindrome(n)
        n -= 1
    end
    return n
end

@assert find_next_palindrome_decreasing(10000) == 9999
@assert find_next_palindrome_decreasing(9998) == 9889

function is_divisible_by_digits(n, digits)
    # Checks if a number is divisible by two x digits number
    max = (10 ^ digits) - 1
    min = minimum([max, trunc(Int, n / (10 ^ digits) + 1)])

    for i in range(max, -1, max-min)
        if n % i == 0 && (n / i) <= max
            return true
        end
    end
    return false
end

@assert is_divisible_by_digits(9999, 2) == false
@assert is_divisible_by_digits(9779, 2) == false
@assert is_divisible_by_digits(9009, 2) == true

function compute_palindrome_by_factors(n, digit)
    while !is_palindrome(n) || !is_divisible_by_digits(n, digit)
        n = find_next_palindrome_decreasing(n-1)
    end
    return n
end
@assert compute_palindrome_by_factors(99*99,2) == 9009
println(compute_palindrome_by_factors(999*999, 3))

max_stars_repo_name: ddugue/project-euler

max_issues_repo_name: ddugue/project-euler

max_forks_repo_name: ddugue/project-euler


text:
__precompile__()

module DiffEqIO

using Reexport, DiffEqBase
@reexport using IterableTables

# package code goes here

end # module


text:
module ODEPrecompileTest
    using ModelingToolkit

    function system(; kwargs...)
        # Define some variables
        @parameters t σ ρ β
        @variables x(t) y(t) z(t)
        D = Differential(t)

        # Define a differential equation
        eqs = [D(x) ~ σ*(y-x),
            D(y) ~ x*(ρ-z)-y,
            D(z) ~ x*y - β*z]

        de = ODESystem(eqs)
        return ODEFunction(de, [x,y,z], [σ,ρ,β]; kwargs...)
    end

    # Build an ODEFunction as part of the module's precompilation. These cases
    # will not work, because the generated RGFs are put into the ModelingToolkit cache.
    const f_bad = system()
    const f_noeval_bad = system(; eval_expression=false)

    # Setting eval_expression=false and eval_module=[this module] will ensure
    # the RGFs are put into our own cache, initialised below.
    using RuntimeGeneratedFunctions
    RuntimeGeneratedFunctions.init(@__MODULE__)
    const f_noeval_good = system(; eval_expression=false, eval_module=@__MODULE__)
end


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CovidSurvey
using CSV
using DataFrames
using Dates
# using RCall
using Turing
using StatsPlots
using Plots.PlotMeasures
using Random
## ===========================================================================
# load data and model
# ============================================================================
fdir = projectdir("reports/2021-04-28/")
fname = "cases=2020-06-01_chains=10_hospit=true_model=contacts2_observ=2020-12-16_preds=CF,CC,CR,CS_seed=30000_sero=true_steps=1000_warmup=1000"
# fname = "cases=2020-06-01_chains=10_hospit=true_model=contacts2_observ=2021-02-07_preds=CF,CC,CR,CS_seed=20000_sero=true_steps=2000_warmup=1000"
## read chain
chain = read(fdir*"chains_"*fname*".jls", Chains)
chain = chain[1:10:end,:,:]
# chain = chain[1:10:end,:,[1,2,3,4,6,7,8,9,10]]
diagnostics = gelmandiag(chain)
## read params
_, ps, _ = parse_savename(fname, parsetypes = (Int, Float64, Bool))
!("preds" in keys(ps)) && (ps["preds"] = nothing)
ps = NamedTuple{Tuple(Symbol.(keys(ps)))}(values(ps))
## read data
data = National.load_data(ps.observ, ps.preds |> CovidSurvey.parse_predictors, ps.cases, ps.hospit, ps.sero; update=false)
turing_data = data.turing_data;
## ===========================================================================
# plot results
# ============================================================================
# make predictions
m_pred = National.model_contacts_v2(turing_data..., true; link=KLogistic(3.), invlink=KLogit(3.) )
gq = Turing.generated_quantities(m_pred, chain)
generated_posterior = vectup2tupvec( reshape(gq, length(gq)) );
observed_daily_cases, expected_daily_deaths, expected_daily_hospit, Rt, expected_seropositive, iar = generated_posterior;
#---------------------------------------------------------------------------
# plot posterior predictions
plotlyjs()
p = National.plot_results(data, generated_posterior...)
savefig(p, "/home/and/tmp/figures/POSTERIOR_"*fname*".html")
# savefig(p, projectdir("figures/tmp", fname*".png") )
run(`firefox $("/home/and/tmp/figures/POSTERIOR_"*fname*".html")`, wait=false)
#---------------------------------------------------------------------------
# plot chain statistics
plotlyjs()
n = filter( x->!occursin(r"(?<!α)\[", x), String.(names(chain)))
p = plot(chain[n]);
savefig(p, "/home/and/tmp/figures/CHAINS_"*fname*".html" )
run(`firefox  $("/home/and/tmp/figures/CHAINS_"*fname*".html")`, wait=false)
#---------------------------------------------------------------------------
# plot more chain statistics (meanplot)
using MCMCChains
plotlyjs()
p = meanplot(chain[n]);
savefig(p, projectdir("/home/and/tmp/figures/MEANPLOT_", fname*".html") )
run(`firefox $(projectdir("/home/and/tmp/figures/MEANPLOT_", fname*".html"))`, wait=false)
## ===========================================================================
# evaluate effect size
# ============================================================================
pgfplotsx()
using LaTeXStrings
## posterior predictions
# scalefontsizes(2.)
p = let
    ocases  = hcat(observed_daily_cases...)#hcat([observed_daily_cases[t] for t = 1:data.num_total]...)
    ehospit  = hcat(expected_daily_hospit...)

    tickdates = data.dates[1]:Month(1):data.dates[end]
    p3 = plot(; legend = :outertopright, ylims=(0,7000), bottom_margin=3mm)
    National.plot_cases_posterior!(p3, data, ocases)
    p4 = plot(; legend = :outertopright, ylims=(0,500), top_margin=3mm)
    National.plot_hospit_posterior!(p4, data, ehospit)
    plot(p3, p4, layout=(2,1), sharex=true, link=:x, size=(1000, 700)) #, size=(1500, 1000)
end
# savefig(p, projectdir("figures/", "posterior_check.png") )
#---------------------------------------------------------------------------
# plot posterior effect size distribution vs. prior
pgfplotsx()
p = let
    xs = chain["α[1]"]|>vec
    p1 = histogram(xs, normalize=true, title="contacts to colleagues", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.03), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = chain["α[2]"]|>vec
    p2 = histogram(xs, normalize=true, title="contacts to strangers", label=nothing, fc=:gray, lw=.1, xlabel="effect size", bottom_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.03), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = chain["α[3]"]|>vec
    p3 = histogram(xs, normalize=true, title="contacts to family", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.03), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    xs = chain["α[4]"]|>vec
    p4 = histogram(xs, normalize=true, title="contacts to friends", label=nothing, fc=:gray, lw=.1, xlabel="effect size", top_margin=3mm);
    plot!(Beta(1,100), xlims=(0,0.03), lw=2, label="prior", legend=:topright);
    vline!([mean(xs)], lw=2, label="mean effect = $(round(mean(xs), sigdigits=2))", c=:black, ls=:dash)
    plot(p1,p2,p3,p4, layout=(2,2))
end
# savefig(p, projectdir("figures/", "total_effect_analysis.png") )
#-----------------------------------------------------------------------------
# evaluate Rt
rt = let
    ib = findfirst(==(Date("2020-05-15")), data.dates)
    ie = findfirst(==(Date(ps.observ)), data.dates)
    Rt_array = hcat(Rt...)[ib:ie,:]

    qs = [quantile(v, [0.025, 0.25, 0.5, 0.75, 0.975]) for v in eachrow(Rt_array)]
    llq, lq, mq, uq, uuq = (eachrow(hcat(qs...))..., )

    date = data.dates[ib:ie]
    DataFrame((;date, llq, lq, mq, uq, uuq))
end
plot(rt.date, rt.mq, ribbon=(rt.mq - rt.llq, rt.uuq - rt.mq), linewidth=0, label="95% quantiles")
plot!(rt.date, rt.mq, ribbon=(rt.mq - rt.lq, rt.uq - rt.mq), linewidth=2, label="50% quantiles", hover=mq)
##
save(projectdir(fdir, "Rt_"*fname*".csv"), rt)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
function generator_kernel_two_level_proxal_first(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u::CuDeviceArray{Float64,1}, x::CuDeviceArray{Float64,1},
    z::CuDeviceArray{Float64,1}, l::CuDeviceArray{Float64,1},
    rho::CuDeviceArray{Float64,1},
    pgmin::CuDeviceArray{Float64,1}, pgmax::CuDeviceArray{Float64,1},
    qgmin::CuDeviceArray{Float64,1}, qgmax::CuDeviceArray{Float64,1},
    c2::CuDeviceArray{Float64,1}, c1::CuDeviceArray{Float64,1},
    tau_prox::Float64, rho_prox::Float64, pg_ref_prox::CuDeviceArray{Float64,1},
    l_next_prox::CuDeviceArray{Float64,1}, pg_next_prox::CuDeviceArray{Float64,1}
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if I <= ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        orig_b = c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx])
        prox_b = -tau_prox*pg_ref_prox[I] + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(orig_b + prox_b) / (2*c2[I]*(baseMVA^2) + rho[pg_idx] + tau_prox + rho_prox))))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))
    end
    return
end

function generator_kernel_two_level_proxal_between(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u::CuDeviceArray{Float64,1}, x::CuDeviceArray{Float64,1},
    z::CuDeviceArray{Float64,1}, l::CuDeviceArray{Float64,1},
    rho::CuDeviceArray{Float64,1},
    pgmin::CuDeviceArray{Float64,1}, pgmax::CuDeviceArray{Float64,1},
    qgmin::CuDeviceArray{Float64,1}, qgmax::CuDeviceArray{Float64,1},
    c2::CuDeviceArray{Float64,1}, c1::CuDeviceArray{Float64,1},
    tau_prox::Float64, rho_prox::Float64,
    pg_ref_prox::CuDeviceArray{Float64,1},
    l_next_prox::CuDeviceArray{Float64,1}, pg_next_prox::CuDeviceArray{Float64,1},
    l_prev_prox::CuDeviceArray{Float64,1}, pg_prev_prox::CuDeviceArray{Float64,1},
    smax_prox::CuDeviceArray{Float64,1}, smin_prox::CuDeviceArray{Float64,1},
    s::CuDeviceArray{Float64,1}
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if I <= ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        # Fill out the solution of qg first.
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))

        orig_b = c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx])
        prox_b = -tau_prox*pg_ref_prox[I] + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        pg_val1 = max(pgmin[I],
                     min(pgmax[I],
                         -(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + rho_prox)))
        s_val1 = -(l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val1)) / rho_prox

        s_val2 = smax_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val2)) + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        pg_val2 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + 2*rho_prox))))
        lg_val2 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val2 + s_val2)

        s_val3 = smin_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val3)) + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        pg_val3 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + 2*rho_prox))))
        lg_val3 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val3 + s_val3)

        if smin_prox[I] <= s_val1 && s_val1 <= smax_prox[I]
            u[pg_idx] = pg_val1
            s[I] = s_val1
        elseif lg_val2 <= 0
            u[pg_idx] = pg_val2
            s[I] = s_val2
        elseif lg_val3 >= 0
            u[pg_idx] = pg_val3
            s[I] = s_val3
        else
            # ERROR!
        end
    end
    return
end

function generator_kernel_two_level_proxal_last(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u::CuDeviceArray{Float64,1}, x::CuDeviceArray{Float64,1},
    z::CuDeviceArray{Float64,1}, l::CuDeviceArray{Float64,1},
    rho::CuDeviceArray{Float64,1},
    pgmin::CuDeviceArray{Float64,1}, pgmax::CuDeviceArray{Float64,1},
    qgmin::CuDeviceArray{Float64,1}, qgmax::CuDeviceArray{Float64,1},
    c2::CuDeviceArray{Float64,1}, c1::CuDeviceArray{Float64,1},
    tau_prox::Float64, rho_prox::Float64,
    pg_ref_prox::CuDeviceArray{Float64,1},
    l_prev_prox::CuDeviceArray{Float64,1}, pg_prev_prox::CuDeviceArray{Float64,1},
    smax_prox::CuDeviceArray{Float64,1}, smin_prox::CuDeviceArray{Float64,1},
    s::CuDeviceArray{Float64,1}
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if I <= ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        # Fill out the solution of qg first.
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))

        orig_b = c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx])
        prox_b = -tau_prox*pg_ref_prox[I]
        pg_val1 = max(pgmin[I],
                     min(pgmax[I],
                         -(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox)))
        s_val1 = -(l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val1)) / rho_prox

        s_val2 = smax_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val2))
        pg_val2 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + rho_prox))))
        lg_val2 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val2 + s_val2)

        s_val3 = smin_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val3))
        pg_val3 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + rho_prox))))
        lg_val3 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val3 + s_val3)

        if smin_prox[I] <= s_val1 && s_val1 <= smax_prox[I]
            u[pg_idx] = pg_val1
            s[I] = s_val1
        elseif lg_val2 <= 0
            u[pg_idx] = pg_val2
            s[I] = s_val2
        elseif lg_val3 >= 0
            u[pg_idx] = pg_val3
            s[I] = s_val3
        else
            # ERROR!
        end
    end
    return
end

#=
function generator_kernel_two_level(
    gen_mod::ProxALGeneratorModel{CuArray{Float64,1}},
    baseMVA::Float64, u::CuArray{Float64,1}, xbar::CuArray{Float64,1},
    zu::CuArray{Float64,1}, lu::CuArray{Float64,1}, rho_u::CuArray{Float64,1}
)
    nblk = div(gen_mod.ngen, 32, RoundUp)
    if gen_mod.t_curr == 1
        tgpu = CUDA.@timed @cuda threads=32 blocks=nblk generator_kernel_two_level_proxal_first(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                    u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                    gen_mod.tau, gen_mod.rho, gen_mod.pg_ref, gen_mod.l_next, gen_mod.pg_next)
    elseif gen_mod.t_curr < gen_mod.T
        tgpu = CUDA.@timed @cuda threads=32 blocks=nblk generator_kernel_two_level_proxal_between(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                    u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                    gen_mod.tau, gen_mod.rho, gen_mod.pg_ref, gen_mod.l_next, gen_mod.pg_next,
                    gen_mod.l_prev, gen_mod.pg_prev, gen_mod.smax, gen_mod.smin, gen_mod.s_curr)
    else
        tgpu = CUDA.@timed @cuda threads=32 blocks=nblk generator_kernel_two_level_proxal_last(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                    u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                    gen_mod.tau, gen_mod.rho, gen_mod.pg_ref, gen_mod.l_prev, gen_mod.pg_prev,
                    gen_mod.smax, gen_mod.smin, gen_mod.s_curr)
    end

    return tgpu

end
=#

function generator_kernel_two_level_proxal_first(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, x, z, l, rho,
    pgmin, pgmax, qgmin, qgmax, c2, c1,
    tau_prox::Float64, rho_prox::Float64, pg_ref_prox,
    l_next_prox, pg_next_prox
)
    for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        orig_b = c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx])
        prox_b = -tau_prox*pg_ref_prox[I] + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        u[pg_idx] = max(pgmin[I],
                        min(pgmax[I],
                            (-(orig_b + prox_b) / (2*c2[I]*(baseMVA^2) + rho[pg_idx] + tau_prox + rho_prox))))
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))
    end
end

function generator_kernel_two_level_proxal_between(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, x, z, l, rho,
    pgmin, pgmax, qgmin, qgmax, c2, c1,
    tau_prox::Float64, rho_prox::Float64,
    pg_ref_prox, l_next_prox, pg_next_prox, l_prev_prox, pg_prev_prox,
    smax_prox, smin_prox, s
)
    for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        # Fill out the solution of qg first.
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))

        orig_b = c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx])
        prox_b = -tau_prox*pg_ref_prox[I] + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        pg_val1 = max(pgmin[I],
                     min(pgmax[I],
                         -(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + rho_prox)))
        s_val1 = -(l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val1)) / rho_prox

        s_val2 = smax_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val2)) + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        pg_val2 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + 2*rho_prox))))
        lg_val2 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val2 + s_val2)

        s_val3 = smin_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val3)) + l_next_prox[I] + rho_prox*(-pg_next_prox[I])
        pg_val3 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + 2*rho_prox))))
        lg_val3 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val3 + s_val3)

        if smin_prox[I] <= s_val1 && s_val1 <= smax_prox[I]
            u[pg_idx] = pg_val1
            s[I] = s_val1
        elseif lg_val2 <= 0
            u[pg_idx] = pg_val2
            s[I] = s_val2
        elseif lg_val3 >= 0
            u[pg_idx] = pg_val3
            s[I] = s_val3
        else
            @printf("[%d] ERROR! Couldn't find a solution for generator at 1<t<T.\n", I)
            @printf("  pg_val1 = %.6e    s_val1 = %.6e\n", pg_val1, s_val1)
            @printf("  pg_val2 = %.6e    s_val2 = %.6e\n", pg_val2, s_val2)
            @printf("  pg_val3 = %.6e    s_val3 = %.6e\n", pg_val3, s_val3)
            # ERROR!
        end
    end
end

function generator_kernel_two_level_proxal_last(
    baseMVA::Float64, ngen::Int, gen_start::Int,
    u, x, z, l, rho,
    pgmin::Array{Float64,1}, pgmax::Array{Float64,1},
    qgmin::Array{Float64,1}, qgmax::Array{Float64,1},
    c2::Array{Float64,1}, c1::Array{Float64,1},
    tau_prox::Float64, rho_prox::Float64,
    pg_ref_prox, l_prev_prox, pg_prev_prox,
    smax_prox, smin_prox, s
)
    for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        # Fill out the solution of qg first.
        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-x[qg_idx] + z[qg_idx]))) / rho[qg_idx]))

        orig_b = c1[I]*baseMVA + l[pg_idx] + rho[pg_idx]*(-x[pg_idx] + z[pg_idx])
        prox_b = -tau_prox*pg_ref_prox[I]
        pg_val1 = max(pgmin[I],
                     min(pgmax[I],
                         -(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox)))
        s_val1 = -(l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val1)) / rho_prox

        s_val2 = smax_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val2))
        pg_val2 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + rho_prox))))
        lg_val2 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val2 + s_val2)

        s_val3 = smin_prox[I]
        prox_b = -tau_prox*pg_ref_prox[I] - (l_prev_prox[I] + rho_prox*(pg_prev_prox[I] + s_val3))
        pg_val3 = max(pgmin[I],
                      min(pgmax[I],
                          (-(orig_b + prox_b) / (2*c2[I]*(baseMVA)^2 + rho[pg_idx] + tau_prox + rho_prox))))
        lg_val3 = l_prev_prox[I] + rho_prox*(pg_prev_prox[I] - pg_val3 + s_val3)

        if smin_prox[I] <= s_val1 && s_val1 <= smax_prox[I]
            u[pg_idx] = pg_val1
            s[I] = s_val1
        elseif lg_val2 <= 0
            u[pg_idx] = pg_val2
            s[I] = s_val2
        elseif lg_val3 >= 0
            u[pg_idx] = pg_val3
            s[I] = s_val3
        else
            @printf("[%d] ERROR! Couldn't find a solution for generator at t=T.\n", I)
            @printf("  pg_val1 = %.6e    s_val1 = %.6e\n", pg_val1, s_val1)
            @printf("  pg_val2 = %.6e    s_val2 = %.6e\n", pg_val2, s_val2)
            @printf("  pg_val3 = %.6e    s_val3 = %.6e\n", pg_val3, s_val3)
            # ERROR!
        end

#        @printf("[%d]: PG[t-1] = %.6e PG = %.6e QG = %.6e s = %.e\n", I, pg_prev_prox[I], u[pg_idx], u[qg_idx], s[I])
    end
end

function generator_kernel_two_level_proxal(ngen::Int, gen_start::Int,
    u::CuDeviceArray{Float64,1}, xbar::CuDeviceArray{Float64,1}, z::CuDeviceArray{Float64,1},
    l::CuDeviceArray{Float64,1}, rho::CuDeviceArray{Float64,1},
    pgmin::CuDeviceArray{Float64,1}, pgmax::CuDeviceArray{Float64,1},
    qgmin::CuDeviceArray{Float64,1}, qgmax::CuDeviceArray{Float64,1},
    smin::CuDeviceArray{Float64,1}, smax::CuDeviceArray{Float64,1}, s::CuDeviceArray{Float64,1},
    _A::CuDeviceArray{Float64,1}, _c::CuDeviceArray{Float64,1})

    tx = threadIdx().x
    I = blockIdx().x

    if I <= ngen
        n = 2
        x = @cuDynamicSharedMem(Float64, n)
        xl = @cuDynamicSharedMem(Float64, n, n*sizeof(Float64))
        xu = @cuDynamicSharedMem(Float64, n, (2*n)*sizeof(Float64))

        A = @cuDynamicSharedMem(Float64, (n,n), (13*n+3)*sizeof(Float64)+(3*n+3)*sizeof(Int))
        c = @cuDynamicSharedMem(Float64, n, (13*n+3+3*n^2)*sizeof(Float64)+(3*n+3)*sizeof(Int))

        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-xbar[qg_idx] + z[qg_idx]))) / rho[qg_idx]))

        A_start = 4*(I-1)
        c_start = 2*(I-1)
        if tx <= n
            @inbounds begin
                for j=1:n
                    A[tx,j] = _A[n*(j-1)+tx + A_start]
                end
                c[tx] = _c[tx + c_start]

                if tx == 1
                    A[1,1] += rho[pg_idx]
                    c[1] += l[pg_idx] + rho[pg_idx]*(-xbar[pg_idx] + z[pg_idx])
                end
            end
        end
        CUDA.sync_threads()

        @inbounds begin
            xl[1] = pgmin[I]
            xu[1] = pgmax[I]
            xl[2] = smin[I]
            xu[2] = smax[I]
            x[1] = min(xu[1], max(xl[1], u[pg_idx]))
            x[2] = min(xu[2], max(xl[2], s[I]))
            CUDA.sync_threads()

            status, minor_iter = tron_qp_kernel(n, 500, 200, 1e-6, 1.0, x, xl, xu, A, c)

            u[pg_idx] = x[1]
            s[I] = x[2]
        end
    end

    return
end

function generator_kernel_two_level_proxal(ngen::Int, gen_start::Int,
    u, xbar, z, l, rho,
    pgmin::Array{Float64,1}, pgmax::Array{Float64,1}, qgmin::Array{Float64,1}, qgmax::Array{Float64,1},
    smin::Array{Float64,1}, smax::Array{Float64,1}, s::Array{Float64,1},
    _A::Array{Float64,1}, _c::Array{Float64,1})

    n = 2
    x = zeros(n)
    xl = zeros(n)
    xu = zeros(n)
    A = zeros(n,n)
    c = zeros(n)

    @inbounds for I=1:ngen
        pg_idx = gen_start + 2*(I-1)
        qg_idx = gen_start + 2*(I-1) + 1

        u[qg_idx] = max(qgmin[I],
                        min(qgmax[I],
                            (-(l[qg_idx] + rho[qg_idx]*(-xbar[qg_idx] + z[qg_idx]))) / rho[qg_idx]))

        A_start = 4*(I-1)
        c_start = 2*(I-1)
        for i=1:n
            for j=1:n
                A[i,j] = _A[n*(j-1)+i + A_start]
            end
            c[i] = _c[i + c_start]
        end

        A[1,1] += rho[pg_idx]
        c[1] += l[pg_idx] + rho[pg_idx]*(-xbar[pg_idx] + z[pg_idx])

        xl[1] = pgmin[I]
        xu[1] = pgmax[I]
        xl[2] = smin[I]
        xu[2] = smax[I]
        x[1] = min(xu[1], max(xl[1], u[pg_idx]))
        x[2] = min(xu[2], max(xl[2], s[I]))

        function eval_f_cb(x)
            f = 0.0
            for j=1:n
                for i=1:n
                    f += x[i]*A[i,j]*x[j]
                end
            end
            f *= 0.5
            for i=1:n
                f += c[i]*x[i]
            end
            return f
        end

        function eval_g_cb(x, g)
            for i=1:n
                g[i] = 0.0
                for j=1:n
                    g[i] += A[i,j]*x[j]
                end
                g[i] += c[i]
            end
            return
        end

        function eval_h_cb(x, mode, rows, cols, _scale, lambda, values)
            if mode == :Structure
                rows[1] = 1; cols[1] = 1
                rows[2] = 2; cols[2] = 1
                rows[3] = 2; cols[3] = 2
            else
                values[1] = A[1,1]
                values[2] = A[2,1]
                values[3] = A[2,2]
            end
            return
        end

        nele_hess = 3
        tron = ExaTron.createProblem(n, xl, xu, nele_hess, eval_f_cb, eval_g_cb, eval_h_cb;
                                     :tol => 1e-6, :matrix_type => :Dense, :max_minor => 200,
                                     :frtol => 1e-12)

        tron.x .= x
        status = ExaTron.solveProblem(tron)

        u[pg_idx] = tron.x[1]
        s[I] = tron.x[2]
    end

    return
end

function generator_kernel_two_level(
    gen_mod::ProxALGeneratorModel{CuArray{Float64,1}},
    baseMVA::Float64, u::CuArray{Float64,1}, xbar::CuArray{Float64,1},
    zu::CuArray{Float64,1}, lu::CuArray{Float64,1}, rho_u::CuArray{Float64,1})

    n = 2
    shmem_size = sizeof(Float64)*(14*n+3+3*n^2) + sizeof(Int)*(3*n+3)

    tgpu1 = build_qp_problem!(baseMVA, gen_mod)
    tgpu2 = CUDA.@timed @cuda threads=32 blocks=gen_mod.ngen shmem=shmem_size generator_kernel_two_level_proxal(
                gen_mod.ngen, gen_mod.gen_start,
                u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax,
                gen_mod.smin, gen_mod.smax, gen_mod.s_curr, gen_mod.Q, gen_mod.c)
    return tgpu2
end

function generator_kernel_two_level(
    gen_mod::ProxALGeneratorModel{Array{Float64,1}},
    baseMVA::Float64, u, xbar, zu, lu, rho_u)
    build_qp_problem!(baseMVA, gen_mod)
    tcpu = @timed generator_kernel_two_level_proxal(gen_mod.ngen, gen_mod.gen_start,
                u, xbar, zu, lu, rho_u,
                gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax,
                gen_mod.smin, gen_mod.smax, gen_mod.s_curr, gen_mod.Q, gen_mod.c)
    return tcpu
end

function build_qp_problem!(baseMVA, gen_mod::ProxALGeneratorModel{TD}) where TD
    for g in 1:gen_mod.ngen
        shift_Q = 4*(g-1)
        shift_c = 2*(g-1)

        # Q[1, 1]
        gen_mod.Q[shift_Q + 1] = gen_mod.Q_ref[shift_Q + 1]
        # Q[2, 1]
        gen_mod.Q[shift_Q + 2] = gen_mod.Q_ref[shift_Q + 2]
        # Q[1, 2]
        gen_mod.Q[shift_Q + 3] = gen_mod.Q_ref[shift_Q + 3]
        # Q[2, 2]
        gen_mod.Q[shift_Q + 4] = gen_mod.Q_ref[shift_Q + 4]

        # c[1]
        gen_mod.c[shift_c + 1] = gen_mod.c_ref[shift_c + 1]
        # c[2]
        gen_mod.c[shift_c + 2] = gen_mod.c_ref[shift_c + 2]

        # τ/2 (pg - pgₚ)^2
        gen_mod.Q[shift_Q + 1] += gen_mod.tau + 2.0 * gen_mod.c2[g]*baseMVA^2
        gen_mod.c[shift_c + 1] += - gen_mod.tau * gen_mod.pg_ref[g] + gen_mod.c1[g]*baseMVA

        # λ₊ (pg - pg₊) + ρ/2 * (pg - pg₊)^2
        if (gen_mod.t_curr < gen_mod.T)
            # Pg
            gen_mod.Q[shift_Q + 1] += gen_mod.rho
            gen_mod.c[shift_c + 1] += gen_mod.l_next[g]
            gen_mod.c[shift_c + 1] -= gen_mod.rho * gen_mod.pg_next[g]
        end

        # λ₋ (pg₋ + s - pg) + ρ/2 * (pg₋ + s - pg)^2
        if (gen_mod.t_curr > 1)
            # Pg
            gen_mod.Q[shift_Q + 1] += gen_mod.rho
            gen_mod.c[shift_c + 1] -= gen_mod.l_prev[g]
            gen_mod.c[shift_c + 1] -= gen_mod.rho * gen_mod.pg_prev[g]

            # Slack
            gen_mod.Q[shift_Q + 4] += gen_mod.rho
            gen_mod.c[shift_c + 2] += gen_mod.l_prev[g]
            gen_mod.c[shift_c + 2] += gen_mod.rho * gen_mod.pg_prev[g]

            # Cross-term (Slack * Pg)
            gen_mod.Q[shift_Q + 2] -= gen_mod.rho
            gen_mod.Q[shift_Q + 3] -= gen_mod.rho
        end
    end
end

function _build_qp_kernel!(
    Q::CuDeviceArray{Float64, 1}, c::CuDeviceArray{Float64, 1}, t_curr, T,
    Q_ref::CuDeviceArray{Float64, 1}, c_ref::CuDeviceArray{Float64, 1},
    baseMVA::Float64, c2::CuDeviceArray{Float64,1}, c1::CuDeviceArray{Float64,1},
    tau_prox::Float64, rho_prox::Float64,
    pg_ref_prox::CuDeviceArray{Float64,1},
    l_next_prox::CuDeviceArray{Float64,1}, pg_next_prox::CuDeviceArray{Float64,1},
    l_prev_prox::CuDeviceArray{Float64,1}, pg_prev_prox::CuDeviceArray{Float64,1},
    ngen::Int,
)
    I = threadIdx().x + (blockDim().x * (blockIdx().x - 1))

    if (I <= ngen)
        shift_Q = 4*(I-1)
        shift_c = 2*(I-1)

        # Q[1, 1]
        Q[shift_Q + 1] = Q_ref[shift_Q + 1]
        # Q[2, 1]
        Q[shift_Q + 2] = Q_ref[shift_Q + 2]
        # Q[1, 2]
        Q[shift_Q + 3] = Q_ref[shift_Q + 3]
        # Q[2, 2]
        Q[shift_Q + 4] = Q_ref[shift_Q + 4]

        # c[1]
        c[shift_c + 1] = c_ref[shift_c + 1]
        # c[2]
        c[shift_c + 2] = c_ref[shift_c + 2]

        # τ/2 (pg - pgₚ)^2
        Q[shift_Q + 1] += tau_prox + 2.0 * c2[I]*baseMVA^2
        c[shift_c + 1] += - tau_prox * pg_ref_prox[I] + c1[I]*baseMVA

        # λ₊ (pg - pg₊) + ρ/2 * (pg - pg₊)^2
        if (t_curr < T)
            # Pg
            Q[shift_Q + 1] += rho_prox
            c[shift_c + 1] += l_next_prox[I]
            c[shift_c + 1] -= rho_prox * pg_next_prox[I]
        end

        # λ₋ (pg₋ + s - pg) + ρ/2 * (pg₋ + s - pg)^2
        if (t_curr > 1)
            # Pg
            Q[shift_Q + 1] += rho_prox
            c[shift_c + 1] -= l_prev_prox[I]
            c[shift_c + 1] -= rho_prox * pg_prev_prox[I]

            # Slack
            Q[shift_Q + 4] += rho_prox
            c[shift_c + 2] += l_prev_prox[I]
            c[shift_c + 2] += rho_prox * pg_prev_prox[I]

            # Cross-term (Slack * Pg)
            Q[shift_Q + 2] -= rho_prox
            Q[shift_Q + 3] -= rho_prox
        end
    end
    return
end

function build_qp_problem!(baseMVA, gen_mod::ProxALGeneratorModel{<:CuArray})
    nblk_gen = div(gen_mod.ngen, 32, RoundUp)
    tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_gen _build_qp_kernel!(
        gen_mod.Q, gen_mod.c, gen_mod.t_curr, gen_mod.T,
        gen_mod.Q_ref, gen_mod.c_ref,
        baseMVA, gen_mod.c2, gen_mod.c1,
        gen_mod.tau, gen_mod.rho, gen_mod.pg_ref,
        gen_mod.l_next, gen_mod.pg_next,
        gen_mod.l_prev, gen_mod.pg_prev, gen_mod.ngen,
    )
    return tgpu
end

#=
function generator_kernel_two_level(
    gen_mod::ProxALGeneratorModel{Array{Float64,1}},
    baseMVA::Float64, u, xbar, zu, lu, rho_u
)
    if gen_mod.t_curr == 1
        tcpu = @timed generator_kernel_two_level_proxal_first(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                    u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                    gen_mod.tau, gen_mod.rho, gen_mod.pg_ref, gen_mod.l_next, gen_mod.pg_next)
    elseif gen_mod.t_curr < gen_mod.T
        tcpu = @timed generator_kernel_two_level_proxal_between(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                    u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                    gen_mod.tau, gen_mod.rho, gen_mod.pg_ref, gen_mod.l_next, gen_mod.pg_next,
                    gen_mod.l_prev, gen_mod.pg_prev, gen_mod.smax, gen_mod.smin, gen_mod.s_curr)
    else
        tcpu = @timed generator_kernel_two_level_proxal_last(baseMVA, gen_mod.ngen, gen_mod.gen_start,
                    u, xbar, zu, lu, rho_u, gen_mod.pgmin, gen_mod.pgmax, gen_mod.qgmin, gen_mod.qgmax, gen_mod.c2, gen_mod.c1,
                    gen_mod.tau, gen_mod.rho, gen_mod.pg_ref, gen_mod.l_prev, gen_mod.pg_prev,
                    gen_mod.smax, gen_mod.smin, gen_mod.s_curr)
    end

    return tcpu
end
=#


text:
"""
A quick realization of MPS-qubit circuit.

For compatibility please also install in Julia Package mode:
pkg> add Yao
pkg> add https://github.com/QuantumBFS/QuAlgorithmZoo.jl.git

Reference:
    Jin-Guo Liu, Yihong Zhang, Yuan Wan and Lei Wang,
    Variational Quantum Eigensolver with Fewer Qubits(arXiv:1902.02663)
"""

module MPSCircuit
export setMPSpar
export DCbuilder
export MPSC
export MPSbuilder

using Yao, Yao.ConstGate
include("Diff.jl")

# Structure of basic parameters.
struct setMPSpar
    nBitA  # Number of lines(bits) in the Extended circuit.
    vBit   # Number of virtual bits in the MPS circuit.
    rBit   # Number of reusable bits in the MPS circuit.
    nBit   # Number of lines(bits) in the MPS circuit.
    nBlock # Number of  MPS blocks in the MPS circuit.

    function setMPSpar(nBitA::Int64, vBit::Int64, rBit::Int64)
        if (nBitA - vBit -rBit) % rBit != 0
            println("Error: nBlock is not integer!")
            return 0
        end
        nBlock = Int((nBitA - vBit) / rBit)
        nBit = rBit + vBit
        new(nBitA, vBit, rBit, nBit, nBlock)
    end
end

# Structure of elements may needed for a Differentiable circuit.
struct DCbuilder
    block   # The block for 1 depth.
    Cblocks # The whole Blocks for n depths.
    body    # The Differentiable circuit for n depths.
    head    # The "head" of the Differentiable circuit(if directly input bits).
    tail    # The "tail" of the Differentiable circuit(2 layers of rotaional gates).

    function DCbuilder(nBit::Int64, ndepth::Int64)
        block = chain(nBit, vcat([repeat(nBit, Rz(0))], [repeat(nBit, Rx(0))], [repeat(nBit, Rz(0))],
                                    [control(nBit, i, (i-1)=>X) for i=nBit:-1:2], control(nBit, 1, nBit=>X)))
        tail = chain(nBit, repeat(nBit, Rz(0)), repeat(nBit, Rx(0)))
        head = block[2:end]
        Cblocks = chain(nBit, repeat([block],ndepth))
        body = chain(nBit, vcat([head], repeat([block],ndepth), [tail]))
        new(block, Cblocks, body, head, tail)
    end
end

# Structure of elements MPSCircuit.jl provides.
struct MPSC
    circuit  # MPS circuit.
    cBlocks  # Array of all the MPS blocks in the MPS circuit.
    cExtend  # The MPS circuit extended back to where it doesn't reuse any qubit.
    cEBlocks # Array of all the MPS blocks in the Extended circuit.
    dGates   # Differentiable gates of the MPS circuit if applicable.
    nBit     # Number of lines(bits) of the MPS circuit.
    nBlock   # Number of blocks in the MPS ciruict.

    function MPSC(MPSblock, nBitA::Int64, vBit::Int64, rBit::Int64=1)
        par2nd = setMPSpar(nBitA, vBit, rBit)
        nBlock = par2nd.nBlock
        nBit = par2nd.nBit
        #println("M1\n")
        MPS = MPSbuilder(nBitA, vBit, rBit, MPSblock)
        #println("M2\n")
        circuit = MPS[1]
        #println("M3\n")
        cExtend = MPS[2]
        dGates = collect_blocks(AbstractDiff, circuit)
        new(circuit, circuit.blocks, cExtend, cExtend.blocks, dGates, nBit, nBlock)
    end

end

# Methods for creating different types of MPS circuits.
function MPSbuilder(nBitA::Int64, vBit::Int64, rBit::Int64, blockT::Tuple{String, Int64})
    par2nd = setMPSpar(nBitA, vBit, rBit)
    nBlock = par2nd.nBlock
    nBit = par2nd.nBit

    if blockT[1] == "DC" # Differentiable circuit.
        depth = blockT[2]
        # println("s1\n")
        swapV1 = chain(nBit, [put(nBit, (inBit,inBit+1)=>ConstGate.SWAP ) for irBit=rBit:-1:1 for inBit=irBit  :(vBit+irBit-1)])
        swap1V = chain(nBit, [put(nBit, (inBit+1,inBit)=>ConstGate.SWAP ) for irBit=1:   rBit for inBit=(vBit+irBit-1):-1:irBit])
        # println("s2\n")
        cBlockHead = chain(nBit, DCbuilder(nBit, depth).block, swapV1) |> diff
        dispatch!(cBlockHead, :random)
        # println("s3\n")
        cBlock =  chain(nBit, vcat([swap1V], cBlockHead.blocks))
        # println("s4\n")
        cBlocks = []
        cEBlocks = []
        for i=1:nBlock
            push!(cBlocks, cBlock)
            # println("s4_1\n")
            push!(cEBlocks, put(nBitA, Tuple((nBitA-i*rBit-vBit+1):(nBitA-(i-1)*rBit),)=>cBlockHead))
            # println("s4_2\n")
        end
        cBlocks[1] = cBlockHead
        circuit = chain(nBit, cBlocks)
        cExtend = chain(nBitA, cEBlocks)
    end

    (circuit, cExtend)
end

function MPSbuilder(nBitA::Int64, vBit::Int64, rBit::Int64, blockT::String)
    par2nd = setMPSpar(nBitA, vBit, rBit)
    nBlock = par2nd.nBlock
    nBit = par2nd.nBit

    if blockT == "CS" # MPS blocks for 1D cluster state.
        if vBit !=1 || rBit !=1
            println("Error: The nBit=vBit+rBit of cluster state MPS blocks should be 2!\n")
            return 0
        end
        cBlocks = []
        push!(cBlocks, chain(nBit, repeat(nBit,H,(2,1)), control(nBit, 1, nBit=>Z)))
        for i =2:nBlock
            push!(cBlocks, chain(nBit,put(nBit, (2,1)=>ConstGate.SWAP), put(nBit, 1=>H), control(nBit, 1, nBit=>Z)))
        end
        circuit = chain(nBit,cBlocks)
        cExtend = chain(nBitA, vcat([repeat(nBitA,H,Tuple(nBitA:-1:1,))], [control(nBitA, i-1, i=>Z) for i=nBitA:-1:2]))
    end

    (circuit, cExtend)
end

end


text:
baseURL = "https://github.com/alan-turing-institute/PDSampler.jl/blob/master/"
tests   = readdir("../test/")
for name in tests
    if occursin(Regex("^ex_"), name) && name[end-2:end]==".jl"
        #########################
        pathin  = "../test/$name"
        URL     = baseURL * "test/$name"
        mdname  = "$(name[1:end-3]).md"
        pathout = "src/examples/$mdname"
        open(pathin) do fi
            open(pathout, "w") do fo
                skip    = true
                iscode  = false
                llclose = false
                for ln in eachline(fi)
                    # Chomp markers
                    if occursin(Regex("^#@startexample"),ln)
                        write(fo, "#"*match(Regex("^#@startexample(.*)"),ln)[1]*"\n\n")
                        write(fo, "(*the code for this example can be found "*
                                  "[here]($URL), note that the doc rendered "* "here was automatically generated, if you "*
                                  "want to fix it, please do it in the "*
                                  "julia code directly*)\n\n")
                        skip = false
                    elseif occursin(Regex("^#@endexample"),ln)
                        break
                    elseif occursin(Regex("^#="),ln)
                        if iscode
                            # close the block
                            write(fo, "\n```\n")
                            iscode = false
                        end
                    elseif occursin(Regex("^=#"),ln)
                        iscode  = true
                        llclose = true
                    else
                        if iscode && llclose
                            write(fo, "\n```julia\n")
                            llclose = false
                        end
                        skip ? nothing : write(fo, ln*"\n")
                    end
                end
            end
        end
    end
end


text:
using ReachabilityBenchmarks, MathematicalSystems, Reachability, Plots

include(@relpath "fom_model.jl")
include(@relpath "fom_specifications.jl")

S = fom_model()
X0, options = fom_specification()
options = Dict(options)

# initial value problem
problem = InitialValueProblem(S, X0)

# partition with 1D blocks
partition = [i:i for i in 1:1006]

# reachability algorithms
algorithm_reach = BFFPSV18(:δ => 1e-3, :vars => [1], :partition => partition,
                           :exp_method => "lazy", :assume_sparse => true)
algorithm_check = BFFPSV18(:δ => 1e-3, :vars => 1:1006, :partition => partition,
                           :exp_method => "lazy", :assume_sparse => true,
                           :lazy_inputs_interval => -1)

# check property
# only verify a small time horizon (full time horizon does not verify)
options_check = copy(options)
options_check[:mode] = "check"
options_check[:T] = 0.001
solution = solve(problem, options_check; op=algorithm_check)
@assert solution.satisfied

# compute flowpipe
options[:mode] = "reach"
solution = solve(problem, options; op=algorithm_reach)

# project flowpipe to time and x₁
solution.options[:plot_vars] = [0, 1]
solution_proj = project(solution)

# plot projection
plot(solution_proj)
savefig("fom")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
################################################################################
#
#  NfRel/NfRel.jl : Relative number field extensions
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: 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) 2017 Claus Fieker, Tommy Hofmann
#
################################################################################

export absolute_field

################################################################################
#
#  Copy
#
################################################################################

function Base.deepcopy_internal(a::NfRelElem{T}, dict::IdDict) where T
  z = NfRelElem{T}(Base.deepcopy_internal(data(a), dict))
  z.parent = parent(a)
  return z
end

################################################################################
#
#  Comply with ring interface
#
################################################################################

elem_type(::Type{NfRel{T}}) where {T} = NfRelElem{T}

elem_type(::NfRel{T}) where {T} = NfRelElem{T}

parent_type(::Type{NfRelElem{T}}) where {T} = NfRel{T}

needs_parentheses(::NfRelElem) = true

isnegative(x::NfRelElem) = isnegative(data(x))

show_minus_one(::Type{NfRelElem{T}}) where {T} = true

function iszero(a::NfRelElem)
  reduce!(a)
  return iszero(data(a))
end

function isone(a::NfRelElem)
  reduce!(a)
  return isone(data(a))
end

zero(K::NfRel) = K(zero(parent(K.pol)))

one(K::NfRel) = K(one(parent(K.pol)))

function zero!(a::NfRelElem)
  zero!(a.data)
  return a
end

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

promote_rule(::Type{NfRelElem{T}}, ::Type{fmpz}) where {T <: RingElement} = NfRelElem{T}

promote_rule(::Type{NfRelElem{T}}, ::Type{fmpq}) where {T <: RingElement} = NfRelElem{T}

promote_rule(::Type{NfRelElem{T}}, ::Type{T}) where {T} = NfRelElem{T}

promote_rule(::Type{NfRelElem{T}}, ::Type{NfRelElem{T}}) where T <: RingElement = NfRelElem{T}

function promote_rule(::Type{NfRelElem{T}}, ::Type{U}) where {T <: RingElement, U <: RingElement}
  promote_rule(T, U) == T ? NfRelElem{T} : Union{}
end

################################################################################
#
#  Order types
#
################################################################################

order_type(K::NfRel{T}) where {T} = NfRelOrd{T, fractional_ideal_type(order_type(base_field(K)))}

order_type(::Type{NfRel{T}}) where {T} = NfRelOrd{T, fractional_ideal_type(order_type(parent_type(T)))}

################################################################################
#
#  Field access
#
################################################################################

@inline base_field(a::NfRel{T}) where {T} = a.base_ring::parent_type(T)

@inline data(a::NfRelElem) = a.data

@inline parent(a::NfRelElem{T}) where {T} = a.parent::NfRel{T}

@inline issimple(a::NfRel) = true

################################################################################
#
#  Coefficient setting and getting
#
################################################################################

@inline coeff(a::NfRelElem{T}, i::Int) where {T} = coeff(a.data, i)

@inline setcoeff!(a::NfRelElem{T}, i::Int, c::T) where {T} = setcoeff!(a.data, i, c)

# copy does not do anything (so far), this is only for compatibility with coeffs(::AbsAlgAssElem)
function coeffs(a::NfRelElem{T}; copy::Bool = false) where {T}
  return T[coeff(a, i) for i = 0:degree(parent(a)) - 1 ]
end

################################################################################
#
#  Constructors
#
################################################################################

(K::NfRel)(x::NfRelElem) = K(base_field(K)(x))

(K::NfRel)(x::nf_elem) = K(base_field(K)(x))

(K::NfRel{T})(x::NfRelElem{T}) where {T} = K(x.data)

(K::NfRel{NfRelElem{T}})(x::NfRelElem{T}) where {T} = K(parent(K.pol)(x))

(K::NfRel{nf_elem})(x::nf_elem) = K(parent(K.pol)(x))

################################################################################
#
#  Degree
#
################################################################################

@inline degree(L::Hecke.NfRel) = degree(L.pol)

################################################################################
#
#  Reduction
#
################################################################################

function reduce!(a::NfRelElem)
  a.data = mod(a.data, parent(a).pol)
  return a
end

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

function Base.show(io::IO, a::NfRel)
  print(io, "Relative number field over\n")
  print(io, a.base_ring, "\n")
  print(io, " with defining polynomial ", a.pol)
end

function _show(io::IO, x::PolyElem, S::String)
   len = length(x)
   if len == 0
      print(io, base_ring(x)(0))
   else
      for i = 1:len - 1
         c = coeff(x, len - i)
         bracket = needs_parentheses(c)
         if !iszero(c)
            if i != 1 && !isnegative(c)
               print(io, "+")
            end
            if !isone(c) && (c != -1 || show_minus_one(typeof(c)))
               if bracket
                  print(io, "(")
               end
               show(io, c)
               if bracket
                  print(io, ")")
               end
               print(io, "*")
            end
            if c == -1 && !show_minus_one(typeof(c))
               print(io, "-")
            end
            print(io, string(S))
            if len - i != 1
               print(io, "^")
               print(io, len - i)
            end
         end
      end
      c = coeff(x, 0)
      bracket = needs_parentheses(c)
      if !iszero(c)
         if len != 1 && !isnegative(c)
            print(io, "+")
         end
         if bracket
            print(io, "(")
         end
         show(io, c)
         if bracket
            print(io, ")")
         end
      end
   end
end

function Base.show(io::IO, a::NfRelElem)
  f = data(a)
  _show(io, f, string(parent(a).S))
end

################################################################################
#
#  Constructors and parent object overloading
#
################################################################################

function NumberField(f::PolyElem{T}, s::String;
                     cached::Bool = false, check::Bool = true)  where {T <: NumFieldElem}
  S = Symbol(s)
  check && !isirreducible(f) && throw(error("Polynomial must be irreducible"))
  K = NfRel{T}(f, S, cached)
  return K, K(gen(parent(f)))
end

function NumberField(f::PolyElem{<: NumFieldElem}; cached::Bool = false, check::Bool = true)
  return NumberField(f, "_\$", cached = cached, check = check)
end

function (K::NfRel{T})(a::Generic.Poly{T}) where T
  z = NfRelElem{T}(mod(a, K.pol))
  z.parent = K
  return z
end

function (K::NfRel{T})(a::T) where T
  parent(a) != base_ring(parent(K.pol)) == error("Cannot coerce")
  z = NfRelElem{T}(parent(K.pol)(a))
  z.parent = K
  return z
end

(K::NfRel)(a::Integer) = K(parent(K.pol)(a))

(K::NfRel)(a::Rational{T}) where {T <: Integer} = K(parent(K.pol)(a))

(K::NfRel)(a::fmpz) = K(parent(K.pol)(a))

(K::NfRel)(a::fmpq) = K(parent(K.pol)(a))

(K::NfRel)() = zero(K)

gen(K::NfRel) = K(gen(parent(K.pol)))

################################################################################
#
#  Unary operators
#
################################################################################

function Base.:(-)(a::NfRelElem)
  return parent(a)(-data(a))
end

################################################################################
#
#  Binary operators
#
################################################################################

function Base.:(+)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  return parent(a)(data(a) + data(b))
end

function Base.:(-)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  return parent(a)(data(a) - data(b))
end

function Base.:(*)(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  return parent(a)(data(a) * data(b))
end

function divexact(a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  b == 0 && error("Element not invertible")
  return a*inv(b)
end

Base.:(//)(a::NfRelElem{T}, b::NfRelElem{T}) where {T} = divexact(a, b)

################################################################################
#
#  Inversion
#
################################################################################

function Base.inv(a::NfRelElem)
  a == 0 && throw(error("Element not invertible"))
  g, s, _ = gcdx(data(a), parent(a).pol)
  @assert g == 1
  return parent(a)(s)
end

################################################################################
#
#  Powering
#
################################################################################

Base.:(^)(a::NfRelElem, n::UInt) = a^Int(n)

function Base.:(^)(a::NfRelElem, n::Int)
  K = parent(a)
  if iszero(a)
    return zero(K)
  end

  if n == 0
    return one(K)
  end

  if n < 0 && iszero(a)
    error("Element is not invertible")
  end

  return K(powmod(data(a), n, K.pol))
end

function Base.:(^)(a::NfRelElem, b::fmpz)
  if b < 0
    return inv(a)^(-b)
  elseif b == 0
    return parent(a)(1)
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  elseif mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

################################################################################
#
#  Comparison
#
################################################################################

function Base.:(==)(a::NfRelElem{T}, b::NfRelElem{T}) where T
  reduce!(a)
  reduce!(b)
  return data(a) == data(b)
end

if !isdefined(Nemo, :promote_rule1)
  function Base.:(==)(a::NfRelElem{T}, b::Union{Integer, Rational}) where T
    return a == parent(a)(b)
  end
end

################################################################################
#
#  Ad hoc operators
#
################################################################################

function Base.:(*)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(data(a) * b)
end

Base.:(*)(a::T, b::NfRelElem{T}) where {T <: NumFieldElem} = b * a

function Base.:(+)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(data(a) + b)
end

Base.:(+)(a::T, b::NfRelElem{T}) where {T <: NumFieldElem} = b + a

function Base.:(-)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(data(a) - b)
end

function Base.:(-)(a::T, b::NfRelElem{T}) where {T <: NumFieldElem}
  return parent(b)(a - data(b))
end

function divexact(a::NfRelElem{T}, b::T) where {T <: NumFieldElem}
  return parent(a)(divexact(data(a), b))
end

Base.:(//)(a::NfRelElem{T}, b::T) where {T <: NumFieldElem} = divexact(a, b)

for F in [fmpz, fmpq, Int]
  @eval begin
    function Base.:(*)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(data(a) * b)
    end

    Base.:(*)(a::$F, b::NfRelElem{T}) where {T <: NumFieldElem} = b * a

    function Base.:(+)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(data(a) + b)
    end

    Base.:(+)(a::$F, b::NfRelElem{T}) where {T <: NumFieldElem} = b + a

    function Base.:(-)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(data(a) - b)
    end

    function Base.:(-)(a::$F, b::NfRelElem{T}) where {T <: NumFieldElem}
      return parent(b)(a - data(b))
    end

    function divexact(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem}
      return parent(a)(divexact(data(a), b))
    end

    Base.:(//)(a::NfRelElem{T}, b::$F) where {T <: NumFieldElem} = divexact(a, b)
  end
end

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

function mul!(c::NfRelElem{T}, a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  mul!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

function mul!(c::NfRelElem{T}, a::NfRelElem{T}, b::T) where {T}
  mul!(c.data, a.data, b)
  c = reduce!(c)
  return c
end

function addeq!(b::NfRelElem{T}, a::NfRelElem{T}) where {T}
  addeq!(b.data, a.data)
  b = reduce!(b)
  return b
end

function add!(c::NfRelElem{T}, a::NfRelElem{T}, b::NfRelElem{T}) where {T}
  add!(c.data, a.data, b.data)
  c = reduce!(c)
  return c
end

################################################################################
#
#  Isomorphism to absolute field (AnticNumberField)
#
################################################################################

#@doc Markdown.doc"""
#    absolute_field(K::NfRel{nf_elem}, cached::Bool = false) -> AnticNumberField, Map, Map
#Given an extension $K/k/Q$, find an isomorphic extension of $Q$.
#"""
function absolute_field(K::NfRel{nf_elem}, cached::Bool = false)
  Ka, a, b, c = _absolute_field(K, cached)
  #return Ka, NfRelToNf(K, Ka, a, b, c), hom(base_ring(K), Ka, a, check = false)
  return Ka, NfToNfRel(Ka, K, a, b, c), hom(base_field(K), Ka, a, check = false)
end

#@doc Markdown.doc"""
#    absolute_field(K::NfRel{NfRelElem}, cached::Bool = false) -> NfRel, Map, Map
#Given an extension $E/K/k$, find an isomorphic extension of $k$.
#In a tower, only the top-most steps are collapsed.
#"""
function absolute_field(K::NfRel{NfRelElem{T}}, cached::Bool = false) where T
  Ka, a, b, c = _absolute_field(K)
  return Ka, NfRelToNfRelRel(Ka, K, a, b, c), hom(base_field(K), Ka, a, check = false)
end

#Trager: p4, Algebraic Factoring and Rational Function Integration
function _absolute_field(K::NfRel, cached::Bool = false)
  f = K.pol
  kx = parent(f)
  k = base_ring(kx)
  Qx = parent(k.pol)

  l = 0
  g = f
  N = norm(g)

  while true
    @assert degree(N) == degree(g) * degree(k)

    if !isconstant(N) && issquarefree(N)
      break
    end

    l += 1

    g = compose(f, gen(kx) - l*gen(k))
    N = norm(g)
  end

  Ka, gKa = NumberField(N, "x", cached = cached, check = false)

  KaT, T = PolynomialRing(Ka, "T", cached = false)

  # map Ka -> K: gen(Ka) -> gen(K)+ k gen(k)

  # gen(k) -> Root(gcd(g, poly(k)))  #gcd should be linear:
  # g in kx = (Q[a])[x]. Want to map x -> gen(Ka), a -> T

  gg = zero(KaT)
  for i=degree(g):-1:0
    auxp = change_ring(Qx(coeff(g, i)), KaT)
    gg = gg*gKa
    add!(gg, gg, auxp)
    #gg = gg*gKa + auxp
  end

  q = gcd(gg, change_ring(k.pol, KaT))
  @assert degree(q) == 1
  al = -trailing_coefficient(q)//lead(q)
  be = gKa - l*al
  ga = gen(K) + l*gen(k)

  #al -> gen(k) in Ka
  #be -> gen(K) in Ka
  #ga -> gen(Ka) in K
  return Ka, al, be, ga
end

function check_parent(a, b)
  return a==b
end

function hash(a::Hecke.NfRelElem{nf_elem}, b::UInt)
  return hash(a.data, b)
end

# Calls simplify on the output of absolute_field
function simplified_absolute_field(L::NfRel{nf_elem}, cached::Bool = false)
  Ka, a, b, c = _absolute_field(L, cached)
  KatoL = NfToNfRel(Ka, L, a, b, c)
  OKa = maximal_order_via_relative(Ka, KatoL)
  K, KtoKa = simplify(Ka)
  KatoK = inv(KtoKa)
  aa = KatoK(a)
  bb = KatoK(b)
  cc = KatoL(KtoKa(gen(K)))
  ktoK = hom(base_field(L), K, aa, check = false)
  KtoL = NfToNfRel(K, L, aa, bb, cc)
  return K, KtoL, ktoK
end

################################################################################
#
#  Basis & representation matrix
#
################################################################################

function basis_matrix(v::Vector{<: NfRelElem})
  K = parent(v[1])
  k = base_field(K)
  z = zero_matrix(k, length(v), degree(K))
  for i in 1:length(v)
    for j in 0:(degree(K) - 1)
      z[i, j + 1] = coeff(v[i], j)
    end
  end
  return z
end

function elem_to_mat_row!(M::Generic.Mat{T}, i::Int, a::NfRelElem{T}) where T
  for c = 1:ncols(M)
    M[i, c] = deepcopy(coeff(a, c - 1))
  end
  return nothing
end

function elem_from_mat_row(L::NfRel{T}, M::Generic.Mat{T}, i::Int) where T
  t = L(1)
  a = L()
  for c = 1:ncols(M)
    a += M[i, c]*t
    mul!(t, t, gen(L))
  end
  return a
end

function representation_matrix(a::NfRelElem)
  L = a.parent
  n = degree(L)
  M = zero_matrix(base_field(L), n, n)
  t = gen(L)
  b = deepcopy(a)
  for i = 1:n-1
    elem_to_mat_row!(M, i, b)
    mul!(b, b, t)
  end
  elem_to_mat_row!(M, n, b)
  return M
end

function norm(a::NfRelElem{nf_elem}, new::Bool = !true)
  if new && ismonic(parent(a).pol) #should be much faster - eventually
    return resultant_mod(a.data, parent(a).pol)
  end
  M = representation_matrix(a)
  return det(M)
end

function norm(a::NfRelElem, new::Bool = true)
  if new && ismonic(parent(a).pol)
    return resultant(a.data, parent(a).pol)
  end
  M = representation_matrix(a)
  return det(M)
end

function tr(a::NfRelElem)
  M = representation_matrix(a)
  return tr(M)
end

################################################################################
#
#  Minimal and characteristic polynomial
#
################################################################################

function _poly_norm_to(f, k::T) where {T}
  if base_ring(f) isa T
    @assert base_ring(f) == k
    return f
  else
    return _poly_norm_to(norm(f), k)
  end
end

#TODO: investigate charpoly/ minpoly from power_sums, aka tr(a^i) and
#      Newton identities
#TODO: cache traces of powers of the generator on the field, then
#      the trace does not need the matrix

@doc Markdown.doc"""
    charpoly(a::NfRelElem) -> PolyElem

Given an element $a$ in an extension $L/K$, this function returns the
characteristic polynomial of $a$ over $K$.
"""
function charpoly(a::NfRelElem)
  M = representation_matrix(a)
  R = PolynomialRing(base_field(parent(a)), cached = false)[1]
  return charpoly(R, M)
end

@doc Markdown.doc"""
    minpoly(a::NfRelElem) -> PolyElem

Given an element $a$ in an extension $L/K$, this function returns the minimal
polynomial of $a$ of $K$.
"""
function minpoly(a::NfRelElem)
  M = representation_matrix(a)
  R = PolynomialRing(base_field(parent(a)), cached = false)[1]
  return minpoly(R, M, false)
end

function charpoly(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})
  f = charpoly(a)
  return _poly_norm_to(f, k)
end

function absolute_charpoly(a::NfRelElem)
  return charpoly(a, FlintQQ)
end

function minpoly(a::NfRelElem, k::Union{NfRel, AnticNumberField, FlintRationalField})

  if parent(a) == k
    R, y  = PolynomialRing(k, cached = false)
    return y - a
  end

  f = minpoly(a)
  while base_ring(f) != k
    f = norm(f)
    g = gcd(f, derivative(f))
    if !isone(g)
      f = divexact(f, g)
    end
  end
  return f
end

function absolute_minpoly(a::NfRelElem)
  return minpoly(a, FlintQQ)
end

#

#

function (R::Generic.PolyRing{nf_elem})(a::NfRelElem{nf_elem})
  if base_ring(R)==base_field(parent(a))
    return a.data
  end
  error("wrong ring")
end

function (R::Generic.PolyRing{NfRelElem{T}})(a::NfRelElem{NfRelElem{T}}) where T
  if base_ring(R)==base_field(parent(a))
    return a.data
  end
  error("wrong ring")
end

################################################################################
#
#  issubfield and isisomorphic
#
################################################################################

function issubfield(K::NfRel, L::NfRel)
  @assert base_field(K) == base_field(L)
  f = K.pol
  g = L.pol
  if mod(degree(g), degree(f)) != 0
    return false, hom(K, L, zero(L), check = false)
  end
  Lx, x = PolynomialRing(L, "x", cached = false)
  fL = Lx()
  for i = 0:degree(f)
    setcoeff!(fL, i, L(coeff(f, i)))
  end
  fac = factor(fL)
  for (a, b) in fac
    if degree(a) == 1
      c1 = coeff(a, 0)
      c2 = coeff(a, 1)
      h = parent(K.pol)(-c1*inv(c2))
      return true, hom(K, L, h(gen(L)), check = false)
    end
  end
  return false, hom(K, L, zero(L), check = false)
end

function isisomorphic(K::NfRel, L::NfRel)
  @assert base_field(K) == base_field(L)
  f = K.pol
  g = L.pol
  if degree(f) != degree(g)
    return false, hom(K, L, zero(L), check = false)
  end
  return issubfield(K, L)
end

################################################################################
#
#  Normal basis
#
################################################################################

# Mostly the same as in the absolute case
function normal_basis(L::NfRel{nf_elem}, check::Bool = false)
  O = EquationOrder(L)
  K = base_field(L)
  OK = base_ring(O)
  d = discriminant(O)

  for p in PrimeIdealsSet(OK, degree(L), -1, indexdivisors = false, ramified = false)
    if valuation(d, p) != 0
      continue
    end

    # Check if p is totally split
    F, mF = ResidueField(OK, p)
    mmF = extend(mF, K)
    Ft, t = PolynomialRing(F, "t", cached = false)
    ft = nf_elem_poly_to_fq_poly(Ft, mmF, L.pol)
    pt = powmod(t, order(F), ft)

    if degree(gcd(ft, pt - t)) == degree(ft)
      # Lift an idempotent of O/pO
      immF = pseudo_inv(mmF)
      fac = factor(ft)
      gt = divexact(ft, first(keys(fac.fac)))
      g = fq_poly_to_nf_elem_poly(parent(L.pol), immF, gt)
      return L(g)
    end
  end
end

################################################################################
#
#  Linear disjointness
#
################################################################################

function islinearly_disjoint(K1::NfRel, K2::NfRel)
  if base_field(K1) != base_field(K2)
    throw(error("Number fields must have the same base field"))
  end

  if gcd(degree(K1), degree(K2)) == 1
    return true
  end

  f = change_base_ring(defining_polynomial(K1), K2)
  return isirreducible(f)
end