text:
using Test
# Unconstrained Results Constructors
n = rand(1:10)
m = rand(1:10)
N = rand(10:10:100)
r = TrajectoryOptimization.UnconstrainedVectorResults(n,m,N)
@test (length(r.X[1]),length(r.X)) == (n,N)
@test (length(r.U[1]),length(r.U)) == (m,N-1)
@test (size(r.K[1])...,length(r.K)) == (m,n,N-1)

r2 = TrajectoryOptimization.UnconstrainedVectorResults(n,m,N)
r.X[1] .= 1:n
copyto!(r2.X,r.X)
@test r2.X[1] == 1:n  # Make sure the copy worked
r2.X[1][1] = 4
@test r.X[1][1] == 1  # Make sure the copies aren't linked

# Constrained Results
p = rand(1:5)
p_N = rand(1:5)
r = ConstrainedVectorResults(n,m,p,N,p_N)
@test (length(r.C[1]),length(r.C)) == (p,N)
@test (size(r.Iμ[1])...,length(r.Iμ)) == (p,p,N)

# Test init_results
N = 10
model, obj = Dynamics.dubinscar
obj_con = ConstrainedObjective(obj, u_min=-10, u_max=10)
solver = Solver(model, obj_con, N=N)
solver.state.infeasible = true
n,m = get_sizes(solver)
p, = get_num_constraints(solver)
p_N, = TrajectoryOptimization.get_num_terminal_constraints(solver)
X = rand(n,N)
U = rand(m,N-1)

results = init_results(solver, X, U)
@test isapprox(to_array(results.X),X)
@test isapprox(to_array(results.U)[1:m,:],U)
@test results.λ[1] == zeros(p)

# Test warm start (partial λs)
λ = [i != N ? rand(p-n) : rand(p_N) for i = 1:N]
results = init_results(solver, X, U, λ=λ)
@test isapprox(to_array(results.X),X)
@test isapprox(to_array(results.U)[1:m,:],U)
@test results.λ[1] == [λ[1]; zeros(n)]
λ[1][1] = 10
@test results.λ[1][1] != 10

# Test warm start (all λs)
λ = [rand(p) for i = 1:N]
push!(λ,rand(n))
results = init_results(solver, X, U, λ=λ)
@test isapprox(to_array(results.X),X)
@test isapprox(to_array(results.U)[1:m,:],U)
@test results.λ[1] == λ[1]
λ[1][1] = 10
@test results.λ[1][1] != 10

λ = [rand(p-1) for i = 1:N-1]
push!(λ,rand(n))
@test_throws ArgumentError init_results(solver, X, U, λ=λ)


text:
using Statistics
using FileIO

mutable struct ObjectParam
	Coord::Array{Int,1}
	Layer::AbstractString
	Img_map::Array{Float64,3}
	Img_ID::Array{Float64,3}
	Source::Bool
	DispPrio::Int
	function ObjectParam(A::Array{Int,1},B::AbstractString,C::Array{Float64,3},D::Array{Float64,3},E::Bool,F::Int)
		if (size(C)[1:2] != size(D)[1:2])
			error(":Img_map: need to have the same size as :Img_ID:")
		end
		new(A,B,C,D,E,F)
	end
end

mutable struct LightParam
	Coord::Array{Int,1}
	Layer::AbstractString
	Intensity::Float64
	GlowAmount::Float64
	Color::Array{Float64,1}
	function LightParam(A::Array{Int,1},B::AbstractString,C::Float64,D::Float64,E::Array{Float64,1})
		if (C<0)
			error("Need :Intensity: cannot be inferior to 0.0")
		end
		if (D<0)
			error("Need :GlowAmount: cannot be inferior to 0.0")
		end
		if (E[1]<0||E[1]>1||E[2]<0||E[2]>1||E[3]<0||E[3]>1)
			error("Need :Color: value between 1.0 and 0.0")
		end
		new(A,B,C,D,E)
	end
end

mutable struct LayerParam
	LayerDist::Float64
	function LayerParam(A::Float64)
		new(A)
	end
end

mutable struct SceneLights
	Light::Dict{AbstractString,LightParam}
	Object::Dict{AbstractString,ObjectParam}
	Layer::Dict{AbstractString,LayerParam}
	function SceneLights()
		new(Dict(),Dict(),Dict("central"=>LayerParam(0.0)))
	end
end

#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function askSceneLights(SL::SceneLights)
	printstyled("\nLights :\n",color=:green)
	for (k,i) in SL.Light
		println("\t['",k,"'] : @",i.Coord,"/layer-> ",i.Layer,".\n\t|> Color : [",i.Color[1],";",i.Color[2],";",i.Color[3],"] ",i.Intensity*100,"%"," |> Amount of glow :",i.GlowAmount*100,"%")
	end
	printstyled("\nObjects :\n",color=:green)
	for (k,i) in SL.Object
		println("\t['",k,"'] : @",i.Coord,"/layer-> ",i.Layer,".\n\t|> Light Source ? : ",i.Source, " |> Display Pritority : ",i.DispPrio)
	end
	printstyled("\nLayers :\n",color=:green)
	for (k,i) in SL.Layer
		println("\tLayer[",k,"] @",i.LayerDist,"m from central layer.")
	end
	println()
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
"""
	\tLoadImgFloat(fileIMG::String)\n
	\t\tReturn (Array{Float64,2},Array{Float64,2},Array{Float64,2})
Load an image , extract the R G B matrices and convert them to Float64 (btw 0.0 and 1.0).\n
Return the (Float) R G B matrices\n
"""
function ImgtoFloat(fileIMG::String)
	img = load(fileIMG)
	(h,l) = size(img)
	r = Array{Float64,3}(undef,h,l,4)
	@inbounds for i = 1:h , j = 1:l
		r[i,j,1] = img[i,j].r
		r[i,j,2] = img[i,j].g
		r[i,j,3] = img[i,j].b
		r[i,j,4] = img[i,j].alpha
	end
	return r
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function IDtoNormal(fileNORM::String)
	norm = load(fileNORM)
	(h,l) = size(norm)
	r = Array{Float64,3}(undef,h,l,3)
	@inbounds for i = 1:h , j = 1:l
		if norm[i,j].alpha == 0.0
			r[i,j,1] = false
			r[i,j,2] = false
			r[i,j,3] = false
		else
		(r[i,j,1],r[i,j,2],r[i,j,3]) = norm_angle(Float64(norm[i,j].r),Float64(norm[i,j].g),Float64(norm[i,j].b))
		end
	end
	return r
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function spwnLight!(List::SceneLights,Name::AbstractString,Coord::Array{Int,1},Layer::AbstractString,Intensity::Float64,Glow::Float64,Color::Array{Float64,1})
	temp = LightParam(Coord,Layer,Intensity,Glow,Color)
	if length(List.Layer) > 0
		if length(filter(x -> (x == Layer), keys(List.Layer))) == 0
			printstyled("\n\t * Layer -",temp.Layer,"- doesn't exists !"," Light *",Name,"* ,not created\n",color=:red)
			return false;
		else
			List.Light[Name] = temp
			return true;
		end
	else
		printstyled("\n\t * Layers not initialized !\n",color=:red)
		return false;
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function spwnObject!(List::SceneLights,Name::AbstractString,Coord::Array{Int,1},Layer::AbstractString,Img::Array{Float64,3},ID::Array{Float64,3};Source=false,DispPrio=0)
	temp = ObjectParam(Coord,Layer,Img,ID,Source,DispPrio)
	if length(List.Layer) > 0
		if length(filter(x -> (x == Layer), keys(List.Layer))) == 0
			printstyled("\n\t * Layer -",temp.Layer,"- doesn't exists !"," Object *",Name,"* ,not created\n",color=:red)
			return false;
		else
			List.Object[Name] = temp
			return true;
		end
	else
		printstyled("\n\t * Layers not initialized !\n",color=:red)
		return false;
	end
end
function spwnObject!(List::SceneLights,Name::AbstractString,Coord::Array{Int,1},Layer::AbstractString,fImg::AbstractString,fID::AbstractString;Source=false,DispPrio=0)
	tImg = ImgtoFloat(fImg)
	tID = IDtoNormal(fID)
	temp = ObjectParam(Coord,Layer,tImg,tID,Source,DispPrio)
	if length(List.Layer) > 0
		if length(filter(x -> (x == Layer), keys(List.Layer))) == 0
			printstyled("\n\t * Layer -",temp.Layer,"- doesn't exists !"," Object *",Name,"* ,not created\n",color=:red)
			return false;
		else
			List.Object[Name] = temp
			return true;
		end
	else
		printstyled("\n\t * Layers not initialized !\n",color=:red)
		return false;
	end
end
function spwnObject!(List::SceneLights,Name::AbstractString,Coord::Array{Int,1},Layer::AbstractString,Img::Array{Float64,3};DispPrio=0)
	tID = Array{Float64,3}(undef,1,1,1)
	temp = ObjectParam(Coord,Layer,Img,tID,true,DispPrio)
	if length(List.Layer) > 0
		if length(filter(x -> (x == Layer), keys(List.Layer))) == 0
			printstyled("\n\t * Layer -",temp.Layer,"- doesn't exists !"," Object *",Name,"* ,not created\n",color=:red)
			return false;
		else
			List.Object[Name] = temp
			return true;
		end
	else
		printstyled("\n\t * Layers not initialized !\n",color=:red)
		return false;
	end
end
function spwnObject!(List::SceneLights,Name::AbstractString,Coord::Array{Int,1},Layer::AbstractString,fImg::AbstractString;DispPrio=0)
	tImg = ImgtoFloat(fImg)
	tID = Array{Float64,3}(undef,1,1,1)
	temp = ObjectParam(Coord,Layer,tImg,tID,true,DispPrio)
	if length(List.Layer) > 0
		if length(filter(x -> (x == Layer), keys(List.Layer))) == 0
			printstyled("\n\t * Layer -",temp.Layer,"- doesn't exists !"," Object *",Name,"* ,not created\n",color=:red)
			return false;
		else
			List.Object[Name] = temp
			return true;
		end
	else
		printstyled("\n\t * Layers not initialized !\n",color=:red)
		return false;
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function spwnLayer!(List::SceneLights,Layer::AbstractString,Dist::Float64)
	temp = (Layer => LayerParam(Dist))
	if length(List.Layer) > 0
		if length(filter(x -> (x == Layer), keys(List.Layer))) > 0
			printstyled("\n\t * Layer ",temp.Layer," already exists !\n",color=:red)
			return false;
		else
			push!(List.Layer,temp)
			return true;
		end
	else
		push!(List.Layer,temp)
		return true;
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function killLight!(List::SceneLights,Name::AbstractString)
	if haskey(List.Light,Name)
		delete!(List.Light,Name)
		return true
	else
		printstyled("Light(s) *",Name,"* not found\n",color=:red)
		return false
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function killObject!(List::SceneLights,Name::AbstractString)
	if haskey(List.Object,Name)
		delete!(List.Object,Name)
		return true
	else
		printstyled("Object(s) *",Name,"* not found\n",color=:red)
		return false
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function killLayer!(List::SceneLights,Layer::Int)
	#Can't kill layer 0
	if Layer == "central"
		printstyled("Can't kill central layer \n",color=:red)
		return false
	else
		temp = filter(x -> (x != Layer), keys(List.Layer))
		if(length(temp) == length(List.Layer))
			printstyled("Can't kill a layer that doesn't exists !\n",color=:red)
			return false
		else
			List.Layer = temp
			return true
		end
	end

end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function tuneLight!(List::SceneLights,Name::AbstractString;
	Coord::Union{String,Array{Int,1}}="x",
	Layer::AbstractString="x",
	Intensity::Union{String,Float64}="x",
	GlowAmount::Union{String,Float64}="x",
	Color::Union{String,Array{Float64,1}}="x"
	)
	if haskey(List.Light,Name)
		if typeof(Coord) == Array{Int,1}
			List.Light[Name].Coord = Coord
		end
		if Layer != "x"
			List.Light[Name].Layer = Layer
		end
		if typeof(Intensity) == Float64
			List.Light[Name].Intensity = Intensity
		end
		if typeof(Intensity) == Float64
			List.Light[Name].GlowAmount = GlowAmount
		end
		if typeof(Color) == Array{Float64,1}
		    List.Light[Name].Color = Color
		end
		return true
	else
		printstyled("Light(s) *",Name,"* not found\n",color=:red)
		return false
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function tuneLayer!(List::SceneLights,Layer::AbstractString;
	LayerDist::Union{String,Float64}="x"
	)

	if length(filter(x -> (x == Layer), keys(List.Layer))) != 0
		if typeof(LayerDist) == Float64
			List.Layer[Layer].LayerDist = LayerDist
		end
		return true
	else
		printstyled("Layer *",Num,"* not found\n",color=:red)
		return false
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function tuneObject!(List::SceneLights,Name::AbstractString;
	Coord::Union{String,Array{Int,1}}="x",
	Layer::AbstractString="x",
	Img_map::Union{String,Array{Float64,3}}="x",
	Img_ID::Union{String,Array{Float64,3}}="x",
	Source::Union{String,Bool}="x",
	DispPrio::Union{String,Int}="x"
	)
	if haskey(List.Object,Name)
		if typeof(Coord) == Array{Int,1}
			List.Object[Name].Coord = Coord
		end
		if Layer != "x"
			List.Object[Name].Layer = Layer
		end
		if typeof(Img_map) == Array{Float64,3}
			List.Object[Name].Img_map = Img_map
		end
		if typeof(Img_ID) == Array{Float64,3}
		    List.Object[Name].Img_ID = Img_ID
		end
		if typeof(Source) == Bool
		    List.Object[Name].Source = Source
		end
		if typeof(DispPrio) == Int
		    List.Object[Name].DispPrio = DispPrio
		end
		return true
	else
		printstyled("Object(s) *",Name,"* not found\n",color=:red)
		return false
	end
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function norm_angle(r::Float64,g::Float64,b::Float64)
	tr = (r != 0 && r < 1.0)-(r==1.0)
	tg = (g != 0 && g < 1.0)-(g==1.0)
	tb = (b==1.0)-(b != 0 && b < 1.0)
	return tr, tg, tb
end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
# Intensity -> Distance in pixel for the Lenticular Light Limit (decrease of 60%) | Glow -> Gaussian form factor
light_distance(intensity,distance,glow) = exp(-(distance/intensity)^(1+glow))
rd(x,y,z) = sqrt(x^2 + y^2 + z^2)
tetha(x,y,z) = acos(z/rd(x,y,z))
function noneg(x)
	if x < 0
		return 0
	else
		return x
	end
end
function phi(x,y)
	if y >= 0
		return acos(x/sqrt(x^2+y^2))
	else
		return 2pi - acos(x/sqrt(x^2+y^2))
	end
end
function n_light(x,y,z,idx,idy,idz)
	if (idx==0) && (idy==0)
		part1 = 0
	else
		part1 = phi(x,y)-phi(idx,idy)
	end
	return part1 + (tetha(x,y,z)-tetha(idx,idy,idz))
end
#n_light(x,y,z,idx,idy,idz) = -cos(phi(x,y) - phi(idx,idy))
#n_light(x,y,z,idx,idy,idz) = -cos(tetha(idx,idy,idz)-tetha(x,y,z))
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function Ray_ID(List::SceneLights,Name::AbstractString)
	#if object is Source, no shading needed
	if List.Object[Name].Source
		return List.Object[Name].Img_map
	end

	(h,l) = size(List.Object[Name].Img_map)

	#black object w/ same opacity features
	current  = copy(List.Object[Name].Img_map)
	current[:,:,1:3] .= 0.0

	debug1 = zeros(size(current)[1],size(current)[2],2)

	#for each light
	for (k,n) in enumerate(keys(sort(List.Light)))
		z = List.Layer[List.Object[Name].Layer].LayerDist - List.Layer[List.Light[n].Layer].LayerDist

		#for each pixels
		for i=1:h , j=1:l
			#ray vector
			x = (List.Object[Name].Coord[1]-List.Light[n].Coord[1])+(j-1)
			y = (List.Object[Name].Coord[2]-List.Light[n].Coord[2])+(i-1)

			#distance between light and pixel
			r = rd(x,y,z)

			#for each color
			if (List.Object[Name].Img_map[i,j,4] != 0)
				#light vs normals
				normal_weight = n_light(x,y,z,List.Object[Name].Img_ID[i,j,1],List.Object[Name].Img_ID[i,j,2],List.Object[Name].Img_ID[i,j,3])

				debug1[i,j,k] = normal_weight

				if ( normal_weight > 0 )
					#light vs distance
					li = ((light_distance(List.Light[n].Intensity,r,List.Light[n].GlowAmount).*List.Light[n].Color) .* List.Object[Name].Img_map[i,j,1:3])
					#applying normals
					li = li .* normal_weight
					#adding light to pixel
					current[i,j,1:3] = current[i,j,1:3] .+ li
				end
			end
		end
	end
	return current,debug1
end

function Ray_ID(List::SceneLights)

end
#::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
function Render_ID()

end

ext: jl


text:
function fevd(estimate::VARRepresentation, H)
	Psi(estimate, H)
	FEVD = [sum(estimate.Psi[i,j,:].^2) for i=1:estimate.vars, j=1:estimate.vars]
	FEVD = hcat([FEVD[:,i]/sum(FEVD[:,i]) for i=1:estimate.vars]...)
	return FEVD
end

function genFEVD(estimate::VARRepresentation, H::Int; nocorr = false)
	Phi(estimate, H)

	A = estimate.Phi
	Σ = deepcopy(estimate.Σ)
	K = estimate.vars

	den = zeros(size(A)[1])
	num = zeros(A[:,:,1])

	if nocorr
		Σ = diagm(diag(deepcopy(estimate.Σ)))
	end

	for i=1:H
		t = A[:,:,i]
		z = t'*Σ
		den += diag(z*t)
		num += z.^2
	end

	θ = [num[i,j]/(den[i]*Σ[j,j]) for i=1:K, j=1:K]
	for i=1:K
		θ[i,:] = θ[i,:]/sum(θ[i,:])
	end
	θ = convert(Array{Float64}, θ)

	return θ
end

# Fourier decomposition of the FEVDs

function fftFEVD(estimate::VARRepresentation, H; nocorr = false, range = nothing)
	k = estimate.vars

	Phi(estimate, H)

	if nocorr
		Σ = diagm(diag(deepcopy(estimate.Σ)))
	else
		Σ = estimate.Σ
	end

	P = chol(Σ)

	if !(typeof(range) <: Tuple)
		range[1] = 1
		range[2] = H
	end

	ft = mapslices(fft, estimate.Phi, [3])
	decomp = [(abs(ft[:,i,z]'*P'[:,j]).^2)[1] / H for i=1:k, j = 1:k, z = 1:H]
	denom = mapslices(sum, decomp[:,:,range[1]:range[2]], [3, 2])

	return  mapslices((x) -> x./denom[:,:,1], decomp, [1,2])
end

function fftGenFEVD(estimate::VARRepresentation, H; nocorr = false, range = nothing)
	k = estimate.vars

	Phi(estimate, H)

	if nocorr
		Σ = diagm(diag(estimate.Σ))
	else
		Σ = estimate.Σ
	end

	P = convert(Array{Float64,2}, chol(Σ))
	ft = mapslices(fft, estimate.Phi, [3])
	decompNew = zeros(k,k,H)
	for i=1:k
		temp = ft[:,i,:]
		for j=1:k
			temp2 = Σ[j,:]
			for h=1:H
				decompNew[i,j,h] = (abs(temp[:,h]'*temp2).^2)[1] / H
			end
		end
	end

	decomp = zeros(k,k,H)

	for i=1:k
		for j=1:k
			for h=1:H
				decomp[i,j,h] = (abs(ft[:,i,h]'*P'[:,j]).^2)[1] / H
			end
		end
	end

	if !(typeof(range) <: Tuple)
		range = [1,H]
	end

	denom = vec(mapslices(sum, decomp[:,:,range[1]:range[2]], [3, 2]))

	θ = [decompNew[i,j,h]/(denom[i]*Σ[j,j]) for i=1:k, j = 1:k, h = 1:H]

	for i=1:k
		s = sum(θ[i,:,range[1]:range[2]])
		for j=1:k
			for h=1:H
				θ[i,j,h] = θ[i,j,h]/s
			end
		end
	end

	return θ
end


text:
@testset "clock" begin
    @testset "basic" begin
        @system SClock(Controller)
        s = instance(SClock)
        @test s.context.clock.time' == 0u"hr"
        @test s.context.clock.tick' === 0
        update!(s)
        @test s.context.clock.time' == 1u"hr"
        @test s.context.clock.tick' === 1
        update!(s)
        @test s.context.clock.time' == 2u"hr"
        @test s.context.clock.tick' === 2
    end

    @testset "config" begin
        @system SClockConfig(Controller)
        o = :Clock => (#=:init => 5,=# :step => 10)
        s = instance(SClockConfig; config=o)
        @test s.context.clock.time' == 0u"hr"
        @test s.context.clock.tick' === 0
        update!(s)
        @test s.context.clock.time' == 10u"hr"
        @test s.context.clock.tick' === 1
        update!(s)
        @test s.context.clock.time' == 20u"hr"
        @test s.context.clock.tick' === 2
    end

    @testset "daily" begin
        @system SClockDaily{Context => Cropbox.DailyContext}(Controller) begin
            a => 1 ~ accumulate
        end
        s = instance(SClockDaily)
        @test s.context isa Cropbox.DailyContext
        @test s.context.clock isa Cropbox.DailyClock
        @test s.context.clock.time' == 0u"d"
        @test s.context.clock.tick' === 0
        update!(s)
        @test s.context.clock.time' == 1u"d"
        @test s.context.clock.tick' === 1
        @test s.a' == 1
    end
end


text:
using DReal

b = ForallVar(Float64, -10.0, 10.0)
c = ForallVar(Float64,  -10.0, 10.0)
d = Var(Float64,-10.,10.)

f(x,y) = (x + y)*d + 6

add!(f(b,c) > 5.0)
@show is_satisfiable()
@show model(d)


text:
function en(an,anm1,Ln)
	return sqrt((1-abs2(an))/((one(an)-abs2(anm1))*(one(an)-abs2(Ln))))
end

function ln(z::Vector{Complex128},a::Vector{Complex128},phinm1::Vector{Complex128})
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(Complex128)
        phinm1star = phinm1
        f1 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
        f2 = [((one(Complex128)-conj(anm1)*z[i])/(one(Complex128)-conj(an)*z[i]))*phinm1star[i] for i=1:length(z)]
        num = sum(f1)
        den = sum(f2)
        r = -num/den
    else
        an = copy(a[end])
        anm1 = copy(a[end-1])
        bProd = length(a)>=3 ? conjBlaschkeProd(z,[a[1:end-2]]) : ones(Complex128,length(z))
        f1 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
        f2 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*phinm1[i]*bProd[i] for i=1:length(z)]
        num = sum(f1)
        den = sum(f2)
        r = -num/den
    end
    return r
end

function phiNext(z::Vector{Complex128},a::Vector{Complex128},phinm1::Vector{Complex128})
    phin = zeros(Complex128,length(phinm1))
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(Complex128)
        L = ln(z,a,phinm1)
        c = en(an,anm1,L)
        for k=1:length(z)
            q = (z[k]-anm1)
            p = (one(Complex128)-conj(anm1)*z[k])
            C = c/(z[k]-an)
            PHI = p*phinm1[k]
            PHISTAR = q*phinm1[k]
            phin[k] = C*(PHI+conj(L)*PHISTAR)
        end
    else
        an = a[end]
        anm1 = a[end-1]
        L = ln(z,a,phinm1)
        c = en(an,anm1,L)
        bProd = length(a)>=3 ? blaschkeProd(z,[a[1:end-2]]) : ones(Complex128,length(z))
        for k=1:length(z)
            q = (z[k]-an)
            p = (one(Complex128)-conj(anm1)*z[k])
            C = c*p/q
            phin[k] = C*(phinm1[k]+conj(L)*bProd[k]*conj(phinm1[k]))
        end
    end
    return (phin, L)
end

function phiNext(z::Vector{Complex128},a::Vector{Complex128},phinm1::Vector{Complex128},L::Complex128)
    phin = zeros(Complex128,length(phinm1))
    if length(a)==1
        an = copy(a[end])
        anm1 = zero(Complex128)
        c = en(an,anm1,L)
        for k=1:length(z)
            q = (z[k]-anm1)
            p = (one(Complex128)-conj(anm1)*z[k])
            C = c/(z[k]-an)
            PHI = p*phinm1[k]
            PHISTAR = q*phinm1[k]
            phin[k] = C*(PHI+conj(L)*PHISTAR)
        end
    else
        an = a[end]
        anm1 = a[end-1]
        c = en(an,anm1,L)
        bProd = length(a)>=3 ? blaschkeProd(z,[a[1:end-2]]) : ones(Complex128,length(z))
        for k=1:length(z)
            q = (z[k]-an)
            p = (one(Complex128)-conj(anm1)*z[k])
            C = c*p/q
            phin[k] = C*(phinm1[k]+conj(L)*bProd[k]*conj(phinm1[k]))
        end
    end
    return phin
end

function blaschkeProd(z::Vector{Complex128},aj::Vector{Complex128})
    r = ones(Complex128,length(z))
    for i=1:length(aj)
        t = [(one(Complex128)-conj(aj[i])*z[k])./(z[k]-aj[i]) for k=1:length(z)]
        r = r.*t
    end
    return r
end

function conjBlaschkeProd(z::Vector{Complex128},aj::Vector{Complex128})
    r = ones(Complex128,length(z))
    for i=1:length(aj)
        t = [(z[k]-aj[i])/(one(Complex128)-conj(aj[i])*z[k]) for k=1:length(z)]
        r = r.*t
    end
    return r
end

function linspace2(s,e,n::Int)
    #this produces a linear spacing between
    #s and e, such that n total elements are
    #in the vector, but the last element
    #is e - spacing, rather than e, like in
    #the normal linspace function
    return [s:((e-s)/(n)):(e-((e-s)/(n)))]
end

function linspace2upsample(s,e,n::Int,m::Int)
    #this function upsamples what linspace2 would
    #return, such that the two vectors share the same
    #points at integer values of the spacing defined
    #by linspace2: (e-s)/m
    return [s:((e-s)/(n*m)):(e-((e-s)/(n*m)))]
end

function genZ(N::Int)
	return exp(im*linspace2(0.0,2*pi,N))
end

function genZArc(N::Int,Div::Int)
    return exp(im*linspace2(0.0,2*pi,N))
end

function genZUpSample(N::Int,Mult::Int)
    return exp(im*linspace2upsample(0.0,2*pi,N,Mult))
end

function genBasis(z::Vector,a::Vector{Complex128})
    P = zeros(Complex128,(length(z),length(a)+1))
    P[:,1] = ones(Complex128,length(z))/sqrt(length(z));
    L = zeros(Complex128,length(a))
    for k=1:length(a)
        (P[:,k+1],L[k]) = phiNext(z,[a[1:k]],P[:,k])
    end
    return (P,L)
end

function genBasisGivenL(z::Vector,a::Vector{Complex128},L::Vector{Complex128},N::Int)
    P = zeros(Complex128,(length(z),length(a)+1))
    P[:,1] = ones(Complex128,length(z))/sqrt(N);
    for k=1:length(a)
        P[:,k+1] = phiNext(z,[a[1:k]],P[:,k],L[k])
    end
    return P
end


text:
function fast_non_dominated_sort!(pop::Vector{T}) where {T}
    F = T[]
    for p in pop
        empty!(p.dom_list)
        p.dom_count = 0
        for q in pop
            if p ⋖ q
                push!(p.dom_list, q)
            elseif q ⋖ p
                p.dom_count += 1
            end
        end
        if p.dom_count == 0
            p.rank = 1
            push!(F, p)
        end
    end
    res = Vector{T}[]
    i = 2
    while !isempty(F)
        Q = T[]
        for p in F
            for q in p.dom_list
                q.dom_count -= 1
                if q.dom_count == 0
                    q.rank = i
                    push!(Q, q)
                end
            end
        end
        i += 1
        push!(res, F)
        F = Q
    end
    res
end

function niching_based_selection(p1, p2, H)

    if p1.CV != p2.CV
        return p1.CV < p2.CV ? p1 : p2
    end

    if p1.ref_point != p2.ref_point
        return rand(Bool) ? p1 : p2
    end

    if p1.rank != p2.rank
        return p1.rank < p2.rank ? p1 : p2
    end

    return p1.distance_ref < p2.distance_ref ? p1 : p2
end

distance(x,h) = norm(-x-dot(-x,h)*h)

function DennisDas(nbobj, p)
    vals = [i for i = 0:1/p:1]
    vcurrent = [[v] for v in vals]
    for i = 2:nbobj
        vnext = Vector{Vector{Float64}}()
        for x in vcurrent
            for v in vals
                if sum(x) + v <= 1.
                    push!(vnext, vcat(x, v))
                end
            end
        end
        vcurrent = vnext
    end
    filter(x -> sum(x) == 1., vcurrent)
end

function normalize_pop!(pop::Vector{T}) where T
    maximums = [pop[1].y...]
    minimums = [pop[1].y...]

    for indiv in pop
        for obj = 1:length(maximums)
            maximums[obj] = max(maximums[obj], indiv.y[obj])
            minimums[obj] = min(minimums[obj], indiv.y[obj])
        end
    end

    for i = 1:length(maximums)
        if maximums[i] == minimums[i]
            minimums[i] -= 1
        end
    end

    for indiv in pop
        indiv.normalized_y = [((maximums[i] - indiv.y[i]) / (maximums[i] - minimums[i])) for i in 1:length(maximums)]
    end

end

function associate_references!(pop, references)
    for ind in pop
        ref = 1
        best_distance = distance(ind.normalized_y, references[1])
        for j_ref = 2:length(references)
            j_distance = distance(ind.normalized_y, references[j_ref])
            if j_distance < best_distance
                best_distance = j_distance
                ref = j_ref
            end
        end
        ind.ref_point = ref
        ind.distance_ref = best_distance
    end
end

function niche_count(pop, references)
    res = zeros(Int, length(references))
    for ind in pop
        res[ind.ref_point] += 1
    end
    res
end

function niching!(P, T, n_count)

    ref = indmin(n_count)

    points = find(x -> x.ref_point == ref, T)

    if isempty(points)
        n_count[ref] = maximum(n_count) + 1
    else
        n_count[ref] += 1
        ipoint = points[1]
        for j = 2:length(points)
            if T[points[j]].distance_ref < T[ipoint].distance_ref
                ipoint = points[j]
            end
        end
        push!(P, T[ipoint])
        deleteat!(T, ipoint)
    end
end


text:
function step!(crr::CorefResolver, params::StepParameters)
    _step!(crr, params)
    push!(crr.crspace.param_log, params)
end

function _step!(crr::CorefResolver, params::InitParams)
   empty!(crr.crspace.param_log)
   init!(crr, params)
end

function _step!(crr::CorefResolver, params::ExtendParams)
    extend!(crr, params)
    crr |> switch!
end

function _step!(crr::CorefResolver, params::MatchParams)
    match!(crr, params)
    crr |> switch!
end

function _step!(crr::CorefResolver, params::IntegrateParams)
    integrate!(crr)
end


text:
validate_devices(partition, ::CPU, devices) = nothing
validate_devices(p, ::CPU, ::Nothing) = nothing

# If no device is specified on the GPU, use only the default device
validate_devices(p, ::GPU, ::Nothing) = 1

function validate_devices(partition, ::GPU, devices)
    @assert length(unique(devices)) <= length(CUDA.devices())
    @assert maximum(devices) <= length(CUDA.devices())
    @assert length(devices) <= length(partition)
    return devices
end

function validate_devices(partition, ::GPU, devices::Number)
    @assert devices <= length(CUDA.devices())
    @assert devices <= length(partition)
    return devices
end
assign_devices(p, ::Nothing) = Tuple(CPU() for i in 1:length(p))

function assign_devices(p::AbstractPartition, dev::Number)
    part     = length(p)
    repeat   = part ÷ dev
    leftover = mod(part, dev)
    devices  = []

    for i in 1:dev
        CUDA.device!(i-1)
        for j in 1:repeat
            push!(devices, CUDA.device())
        end
        if i <= leftover
            push!(devices, CUDA.device())
        end
    end
    return Tuple(devices)
end

function assign_devices(p::AbstractPartition, dev::Tuple)
    part     = length(p)
    repeat   = part ÷ length(dev)
    leftover = mod(part, length(dev))
    devices  = []

    for i in 1:length(dev)
        CUDA.device!(dev[i])
        for j in 1:repeat
            push!(devices, CUDA.device())
        end
        if i <= leftover
            push!(devices, CUDA.device())
        end
    end
    return Tuple(devices)
end


text:
using Random, Statistics

@testset "Assortativity" begin
    # Test definition of assortativity as Pearson correlation coefficient
    # between excess of degrees
    @testset "Small graphs" for n = 5:10
        @test @inferred assortativity(wheel_graph(n)) ≈ -1/3
    end
    @testset "Directed ($seed)" for seed in [1, 2, 3], (n, ne) in [(14, 18), (10, 22), (7, 16)]
        g = erdos_renyi(n, ne; is_directed=true, seed=seed)
        assort = assortativity(g)
        x = [outdegree(g, src(d))-1 for d in edges(g)]
        y = [indegree(g, dst(d))-1 for d in edges(g)]
        @test @inferred assort ≈ cor(x, y)
    end
    @testset "Undirected ($seed)" for seed in [1, 2, 3], (n, ne) in [(14, 18), (10, 22), (7, 16)]
        g = erdos_renyi(n, ne; is_directed=false, seed=seed)
        assort = assortativity(g)
        x = [outdegree(g, src(d))-1 for d in edges(g)]
        y = [indegree(g, dst(d))-1 for d in edges(g)]
        @test @inferred assort ≈ cor([x;y], [y;x])
    end
end


text:
module  TestAlleleFreqs

using PopGen
using DataFrames
using Test

cats = @nancycats;
sharks = @gulfsharks;

@testset "PopGen.allele_freq" begin
    @test length(PopGen.allele_freq(cats)) == 9
    @test eltype(PopGen.allele_freq(cats)) == Dict{Int16,Float64}
    @test typeof(PopGen.allele_freq(cats.loci.genotype)) == Dict{Int16,Float64}
    df_cats = DataFrames.combine(
        groupby(cats.loci, :population),
        :genotype => PopGen.allele_freq => :genos
    )
    @test size(df_cats) == (17,2)
    @test eltype(df_cats.genos) == Dict{Int16,Float64}
    @test length(PopGen.allele_freq(cats, "fca8")) == 16
    @test typeof(PopGen.allele_freq(cats, "fca8")) == Dict{Int16,Float64}

    @test length(PopGen.allele_freq(sharks)) == 2209
    @test eltype(PopGen.allele_freq(sharks)) == Dict{Int8,Float64}
    @test typeof(PopGen.allele_freq(sharks.loci.genotype)) == Dict{Int8,Float64}
    df_sharks = DataFrames.combine(
        groupby(sharks.loci, :population),
        :genotype => PopGen.allele_freq => :genos
    )
    @test size(df_sharks) == (7,2)
    @test eltype(df_sharks.genos) == Dict{Int8,Float64}
    @test length(PopGen.allele_freq(sharks, "contig_2784")) == 2
    @test typeof(PopGen.allele_freq(sharks, "contig_2784")) == Dict{Int8,Float64}
end

@testset "other PopGen.allele_freqs" begin
    @test PopGen.allele_freq_vec([(1,1), (2,2), (2,1)]) == [0.5, 0.5]
    @test PopGen.allele_freq_vec(missing) === missing
end

@testset "geno counts" begin
    @test length(PopGen.geno_count_observed(cats.loci.genotype)) == 295
    @test typeof(PopGen.geno_count_observed(cats.loci.genotype)) == Dict{Tuple,Float32}
    @test length(PopGen.geno_count_expected(cats.loci.genotype)) == 6241
    @test typeof(PopGen.geno_count_expected(cats.loci.genotype)) == Dict{Tuple,Float64}

    @test length(PopGen.geno_count_observed(sharks.loci.genotype)) == 12
    @test typeof(PopGen.geno_count_observed(sharks.loci.genotype)) == Dict{Tuple,Float32}
    @test length(PopGen.geno_count_expected(sharks.loci.genotype)) == 25
    @test typeof(PopGen.geno_count_expected(sharks.loci.genotype)) == Dict{Tuple,Float64}
end

@testset "geno freqs" begin
@test length(PopGen.geno_freq(cats.loci.genotype)) == 295
@test typeof(PopGen.geno_count_observed(cats.loci.genotype)) == Dict{Tuple,Float32}
@test length(PopGen.geno_freq(cats, "fca8")) == 51
@test typeof(PopGen.geno_freq(cats, "fca8")) == Dict{Tuple,Float32}
@test typeof(PopGen.geno_freq(cats, "fca8", population = true)) == DataFrame
@test size(PopGen.geno_freq(cats, "fca8", population = true)) == (17,2)

@test length(PopGen.geno_freq_expected(cats.loci.genotype)) == 6241
@test typeof(PopGen.geno_freq_expected(cats.loci.genotype)) == Dict{Tuple,Float64}
@test length(PopGen.geno_freq_expected(cats, "fca8")) == 256
@test typeof(PopGen.geno_freq_expected(cats, "fca8")) == Dict{Tuple,Float64}
@test typeof(PopGen.geno_freq_expected(cats, "fca8", population = true)) == DataFrame
@test size(PopGen.geno_freq_expected(cats, "fca8", population = true)) == (17,2)

@test length(PopGen.geno_freq(sharks.loci.genotype)) == 12
@test typeof(PopGen.geno_count_observed(sharks.loci.genotype)) == Dict{Tuple,Float32}
@test length(PopGen.geno_freq(sharks, "contig_2784")) == 2
@test typeof(PopGen.geno_freq(sharks, "contig_2784")) == Dict{Tuple,Float32}
@test typeof(PopGen.geno_freq(sharks, "contig_2784", population = true)) == DataFrame
@test size(PopGen.geno_freq(sharks, "contig_2784", population = true)) == (7,2)

@test length(PopGen.geno_freq_expected(sharks.loci.genotype)) == 25
@test typeof(PopGen.geno_freq_expected(sharks.loci.genotype)) == Dict{Tuple,Float64}
@test length(PopGen.geno_freq_expected(sharks, "contig_2784")) == 4
@test typeof(PopGen.geno_freq_expected(sharks, "contig_2784")) == Dict{Tuple,Float64}
@test typeof(PopGen.geno_freq_expected(sharks, "contig_2784", population = true)) == DataFrame
@test size(PopGen.geno_freq_expected(sharks, "contig_2784", population = true)) == (7,2)
end

end # module


text:
export qft_circuit, simple_state_prep_circuit, ghz_circuit

"""
    qft_circuit(n::Integer)

Creates a QFT circuit
"""
function qft_circuit(n::Integer)
    qc = Gates.QuantumCircuit(n)
    for i in 1:n
        push!(qc.ops, Gates.GateOperator([i], Gates.h))
        for j in i:n-1
            push!(qc.ops, Gates.GateOperator([j+1 , i], Gates.rk(j-i+2)))
        end
    end
    for i = 1:convert(Integer, floor(n//2))
        push!(qc.ops, Gates.GateOperator([i , n + 1 - i], Gates.swap))
    end
    qc
end

"""
    function simple_state_prep_circuit(qubits::Integer, depth::Integer)

Create a simple state preparation circuit with depth layers
"""
function simple_state_prep_circuit(qubits::Integer, depth::Integer)
    qc = Gates.QuantumCircuit(qubits)
    for d in 1:depth
        angles = rand(qubits, 3) * 2π
        for i in 1:qubits
            gate = Gates.u3(angles[i, :]...)
            push!(qc.ops, Gates.GateOperator([i], gate))
        end
        for i in ((d % 2)+1):2:qubits-1
            push!(qc.ops, Gates.GateOperator([i , i+1], Gates.cx))
        end
    end
    qc
end

"""
    ghz_circuit(n::Integer)

Creates a GHZ preparation circuit
"""
function ghz_circuit(n::Integer)
    qc = Gates.QuantumCircuit(n)
    push!(qc.ops, Gates.GateOperator([1], Gates.h))
    for i in 1:n-1
        push!(qc.ops, Gates.GateOperator([i, i+1], Gates.cx))
    end
    qc
end


text:
using FractionalDiffEq, Plots

a=1; mu=4

function Brusselator(du, u, p, t)
    du[1] = a-(mu+1)*u[1]+u[1]^2*u[2]
    du[2] = mu*u[1]-u[1]^2*u[2]
    #du
end
alpha=[0.8; 0.8]
t0=0; tfinal=50; y0=[0.2; 0.03]
h=0.01; tspan = (t0, tfinal)
prob = FODESystem(Brusselator, alpha, y0, tspan)
#(t, y) = solve(prob, h, FLMMBDF())
#(t, y) = solve(prob, h, FLMMNewtonGregory())
(t, y) = solve(prob, h, FLMMTrap())

plot(y[1, :], y[2, :])


text:
include("../gcc_common.jl")
build_and_upload_gcc(v"10.2.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
import ..UncertainValues: CertainValue
using RecipesBase

@recipe f(::Type{CertainValue{T}}, x::CertainValue{T}) where {T} = [x.value]

@recipe function f(certainvals::Vector{CertainValue})
    @series begin
        [val.value for val in certainvals]
    end
end

@recipe function f(certainvals::Vararg{CertainValue,N}) where {N}
    @series begin
        [val.value for val in certainvals]
    end
end


text:
mutable struct Clock <: AbstractClock
  player::AbstractPlayer

  last_move::Union{DateTime, Nullable{DateTime}}

  drop::Union{UUIDs.UUID, Nullable{UUIDs.UUID}}
  lock::Union{UUIDs.UUID, Nullable{UUIDs.UUID}}
  hold::Union{UUIDs.UUID, Nullable{UUIDs.UUID}}
  drag::Union{UUIDs.UUID, Nullable{UUIDs.UUID}}

  push::Union{UUIDs.UUID, Nullable{UUIDs.UUID}}
end

function Clock(cur_player::AbstractPlayer=Player())
  cur_clock = Clock(
    cur_player,
    Nullable{DateTime}(),
    Nullable{UUIDs.UUID}(),
    Nullable{UUIDs.UUID}(),
    Nullable{UUIDs.UUID}(),
    Nullable{UUIDs.UUID}(),
    Nullable{UUIDs.UUID}()
  )

  cur_player.clock = cur_clock

  cur_clock
end


text:
module UCI

using ..Chess

export BestMoveInfo, BoundType, Engine, Option, OptionType, OptionValue, Score, SearchInfo

export mpvsearch,
    newgame,
    parsebestmove,
    parsesearchinfo,
    quit,
    runengine,
    search,
    sendcommand,
    sendisready,
    setboard,
    setoption,
    touci

"""
    OptionType

Type representing an UCI option type. This is an enum with values corresponding
to the option types defined in the UCI protocol: `check`, `spin`, `combo`,
`button` and `string`.
"""
@enum OptionType begin
    check
    spin
    combo
    button
    string
end

"""
    OptionValue

Type representing the value of a UCI option. This is a union type containing
the types `Nothing` (for options of type `button`), `Bool` (for options of type
`check`, `Int` (for options of type `spin`) and `String` (for options of type
`combo` or `string`).)
"""
const OptionValue = Union{Nothing,Bool,Int,String}

"""
    Option

Type representing a UCI option. This is a mutable struct with the following slots:

- `name`: A `String`, the name of the option, as provided by the engine.
- `type`: An `OptionType`, the type of the option, as provided by the engine.
- `defaultValue`: An `OptionValue`, as provided by the engine.
- `value`: An `OptionValue`, the current value of the option.
- `min`: The minimum value of this option. Only used for options of type `spin`.
- `max`: The maximum value of this option. Only used for options of type `spin`.
- `combovals`: Vector of the possible values for this option. Only used for
  options of type `Combo`.
"""
mutable struct Option
    name::String
    type::OptionType
    defaultvalue::OptionValue
    value::OptionValue
    min::Int
    max::Int
    combovals::Vector{String}
end

function parseoptionname(s::String)::Tuple{String,String}
    @assert startswith(s, "name ")
    result = ""
    for c in s[6:end]
        result *= c
        endswith(result, " type") && break
    end
    (result[1:end-5], s[length(result)+2:end])
end

function parseoptiontype(s::String)::Tuple{OptionType,String}
    @assert startswith(s, "type ")
    result = ""
    for c in s[6:end]
        c == ' ' && break
        result *= c
    end
    optionsbyname = Dict(
        "check" => check,
        "spin" => spin,
        "combo" => combo,
        "button" => button,
        "string" => string,
    )
    (optionsbyname[result], s[length(result)+6:end])
end

function parseoptiondefault(ot::OptionType, s::String)::Tuple{OptionValue,String}
    if ot == button
        (nothing, s)
    else
        @assert(startswith(s, " default "))
        s = s[10:end]
        if ot == check
            (startswith(s, "true"), s)
        elseif ot == spin
            result = ""
            for c in s
                c == ' ' && break
                result *= c
            end
            (parse(Int, result), s[length(result)+1:end])
        elseif ot == combo
            result = ""
            for c in s
                endswith(result, " var ") && break
                result *= c
            end
            (result[1:end-5], s[length(result)-3:end])
        elseif ot == string
            (s, s)
        end
    end
end

function parsespinminmax(s::String)::Tuple{Int,Int}
    @assert startswith(s, " min ")
    tokens = split(s[2:end], r"\s+")
    @assert length(tokens) >= 4
    @assert tokens[1] == "min" && tokens[3] == "max"
    (parse(Int, tokens[2]), parse(Int, tokens[4]))
end

function parsecombovar(s::String)::Tuple{String,String}
    @assert startswith(s, "var ")
    result = ""
    for c in s[5:end]
        endswith(result, " var ") && break
        result *= c
    end
    if endswith(result, " var ")
        (result[1:end-5], s[length(result)+1:end])
    else
        (result, "")
    end
end

function parsecombovars(s::String)::Vector{String}
    result = String[]
    while startswith(s, "var ")
        (next, rest) = parsecombovar(s)
        push!(result, next)
        s = rest
    end
    result
end

function parseoption(s::String)::Option
    @assert startswith(s, "option name")
    s = s[8:end]
    (name, s) = parseoptionname(s)
    (otype, s) = parseoptiontype(s)
    (default, s) = parseoptiondefault(otype, s)
    (min, max) = otype == spin ? parsespinminmax(s) : (0, 0)
    combovals = otype == combo ? parsecombovars(s) : String[]
    Option(name, otype, default, default, min, max, combovals)
end

"""
    Engine

Type representing a UCI chess engine.

This is a struct with the following slots:

- `name`: The engine name, as provided by the engine in response to the `uci`
  command.
- `author`: The engine author name, as provided by the engine in response to
  the `uci` command.
- `options`: The UCI options for this engine. This is a dictionary mapping
  option names (`String`s) to options (instances of the `Option` type).
- `io`: A `Base.Process` object used to communicate with the engine.

Engines are created by calling the `runengine` function, which takes a pathname
for an UCI engine as input, runs the engine, and returns an `Engine` object.

# Examples

The below is a typical interaction with a UCI engine. The example assumes that
you have a UCI engine with the file name `stockfish` somewhere in your `PATH`.

```julia-repl
julia> sf = runengine("stockfish");

julia> setoption(sf, "Hash", 128)

julia> setboard(sf, fromfen("1kbr3r/pp6/8/P1n2ppq/2N3n1/R3Q1P1/3B1P2/2R2BK1 w - -"))

julia> search(sf, "go depth 18", infoaction=println)
info depth 1 seldepth 1 multipv 1 score cp -842 nodes 88 nps 88000 tbhits 0 time 1 pv f1g2 g4e3 d2e3
info depth 2 seldepth 2 multipv 1 score cp -842 nodes 207 nps 103500 tbhits 0 time 2 pv f1g2 g4e3
info depth 3 seldepth 3 multipv 1 score cp -844 nodes 270 nps 135000 tbhits 0 time 2 pv f1g2 g4e3 d2e3 d8d1 c1d1 h5d1 g2f1
info depth 4 seldepth 5 multipv 1 score cp -844 nodes 367 nps 183500 tbhits 0 time 2 pv f1g2 g4e3 d2e3 d8d1 c1d1
info depth 5 seldepth 7 multipv 1 score cp -953 nodes 866 nps 433000 tbhits 0 time 2 pv f1g2 h5h2 g1f1 g4e3 d2e3 c5d3 e3g5
info depth 6 seldepth 8 multipv 1 score cp -1060 nodes 1507 nps 753500 tbhits 0 time 2 pv f1g2 g4e3 d2e3 c5d3 e3d2 d3c1
info depth 7 seldepth 11 multipv 1 score cp -876 nodes 1995 nps 665000 tbhits 0 time 3 pv f1g2 g4e3
info depth 8 seldepth 10 multipv 1 score cp -882 nodes 2771 nps 923666 tbhits 0 time 3 pv f1g2 g4e3 d2e3 c8e6 e3c5 d8d1 c1d1 h5d1 g2f1 e6c4
info depth 9 seldepth 15 multipv 1 score cp -1068 nodes 12059 nps 1339888 tbhits 0 time 9 pv f1g2 g4e3 a3e3 c8e6 e3e6 c5e6 d2e1 d8d1 c1d1 h5d1
info depth 10 seldepth 15 multipv 1 score cp -1050 nodes 17558 nps 1463166 tbhits 0 time 12 pv f1g2 g4e3 a3e3 c8e6 g1f1 h5g4 e3e6 c5e6 c4e3 g4a4 d2c3
info depth 11 seldepth 17 multipv 1 score cp -1013 nodes 23543 nps 1569533 tbhits 0 time 15 pv f1g2 g4e3 a3e3 c8e6 g1f1 h5g4 e3e6 c5e6 c4e3 g4a4 d2c3
info depth 12 seldepth 20 multipv 1 score cp -995 nodes 48497 nps 1672310 tbhits 0 time 29 pv f1g2 g4e3 a3e3 c8e6 g1f1 c5e4 d2e1 d8c8 g3g4 h5g4 c4e5 g4h4 c1c8 h8c8
info depth 13 seldepth 29 multipv 1 score cp -1116 nodes 205726 nps 1959295 tbhits 0 time 105 pv f1g2 g4e3 f2e3 c5e4 a3a2 h5h2 g1f1 e4g3 f1f2 g3e4 f2f1 e4d2 c4d2 h2g3 d2f3 d8d3 f1g1 d3e3 f3d4 e3e8 a2e2
info depth 14 seldepth 28 multipv 1 score cp -148 nodes 247247 nps 1993927 tbhits 0 time 124 pv e3f4 g5f4 d2f4 b8a8 c4b6 a7b6 a5b6 c5a6
info depth 15 seldepth 16 multipv 1 score cp 120 nodes 248911 nps 1991288 tbhits 0 time 125 pv e3f4 g5f4 d2f4 g4e5 f4e5 d8d6 e5d6 b8a8 f1g2 c8e6 d6c5
info depth 16 seldepth 24 multipv 1 score cp 1117 nodes 251829 nps 1982905 tbhits 0 time 127 pv e3f4 g4e5 f4e5 d8d6 e5d6 b8a8 f1g2 h5h2 g1f1 c5e4 g2e4 f5e4 d2g5 c8h3 f1e2
info depth 17 seldepth 22 multipv 1 score cp 1500 nodes 258707 nps 1990053 tbhits 0 time 130 pv e3f4 g4e5 f4e5 d8d6 e5d6 b8a8 f1g2 h5h2 g1f1 c5e4 g2e4 h2h3 e4g2
info depth 18 seldepth 22 multipv 1 score mate 11 nodes 281736 nps 1984056 tbhits 0 time 142 pv e3f4 g5f4 d2f4 g4e5 f4e5 d8d6 e5d6 b8a8 c4b6 a7b6 a5b6 c5a6 c1c8 h8c8
BestMoveInfo (best=e3f4, ponder=g5f4)
```
"""
mutable struct Engine
    name::String
    author::String
    options::Dict{String,Option}
    io::Base.Process
end

Base.show(io::IO, engine::Engine) = print(io, "Engine: $(engine.name)")

"""
    function runengine(path::String)::Engine

Runs the engine at the specified path, returning an `Engine`.
"""
function runengine(path::String)::Engine
    process = open(`$path`, "r+")
    options = Dict{String,Option}()
    println(process, "uci")
    name = path
    author = ""
    while true
        line = readline(process)
        line == "uciok" && break
        if startswith(line, "id name ")
            name = line[9:end]
        elseif startswith(line, "id author ")
            author = line[11:end]
        elseif startswith(line, "option name")
            o = parseoption(line)
            options[o.name] = o
        end
    end
    Engine(name, author, options, process)
end

function isvalidvalueforoption(o::Option, v::OptionValue)::Bool
    if o.type == string
        typeof(v) == String
    elseif o.type == combo
        typeof(v) == String && any(cv -> cv == v, o.combovals)
    elseif o.type == check
        typeof(v) == Bool
    elseif o.type == spin
        typeof(v) == Int && v >= o.min && v <= o.max
    elseif o.type == button
        isnothing(v)
    else
        false
    end
end

"""
    sendcommand(e::Engine, cmd::String)

Sends the UCI command `cmd` to the engine `e`.
"""
function sendcommand(e::Engine, cmd::String)
    println(e.io, cmd)
end

"""
    sendisready(e::Engine)::Bool

Sends the engine `e` an "isready" command and waits for the "readyok" response.

Returns `true` on success, `false` on failure (i.e. if the engine replies with
anything other than "readyok").
"""
function sendisready(e::Engine)::Bool
    sendcommand(e, "isready")
    readline(e.io) == "readyok"
end

"""
    newgame(e::Engine)

Instructs the engine that a new game is about to begin.
"""
function newgame(e::Engine)
    sendcommand(e, "ucinewgame")
    sendisready(e)
end

"""
    function setoption(e::Engine, name::String, value::OptionValue = nothing)

Sets the UCI option named `name` to the new value `value`.

Throws an error if the engine `e` does not have an option with the provided
name, or if the value is incompatible with the type of the option.
"""
function setoption(e::Engine, name::String, value::OptionValue = nothing)
    if !haskey(e.options, name)
        throw("Engine $(e.name) has no option named $name")
    elseif !isvalidvalueforoption(e.options[name], value)
        throw("Invalid value $value for option $name")
    else
        e.options[name].value = value
        if isnothing(value)
            sendcommand(e, "setoption name $name")
        else
            sendcommand(e, "setoption name $name value $value")
        end
    end
    sendisready(e)
end

"""
    function quit(e::Engine)

Sends the UCI engine `e` the `"quit"` command.
"""
function quit(e::Engine)
    sendcommand(e, "quit")
end

"""
    BoundType

An enum type representing the score bound types `lower`, `upper` and `exact`.
"""
@enum BoundType begin
    lower = 1
    upper = 2
    exact = 3
end

"""
    Score

A struct type representing a score returned by a UCI engine.

The struct has the following slots:

- `value`: An `Int` representing the score value.
- `ismate`: A `Bool` that tells whether this is a mate score.
- `bound`: A `BoundType`, indicating whether this score is a lower bound, an
  upper bound, or an exact value.
"""
struct Score
    value::Int
    ismate::Bool
    bound::BoundType
end

"""
    BestMoveInfo

A struct representing the contents of a UCI engine's `bestmove` output.

Contains the following slots:

- `bestmove`: The `Move` returned by the engine as the best move.
- `ponder`: The engine's expected reply to the best move, a `Move` or
  `nothing`.
"""
struct BestMoveInfo
    bestmove::Move
    ponder::Union{Move,Nothing}
end

function Base.show(io::IO, bmi::BestMoveInfo)
    if isnothing(bmi.ponder)
        print(io, "BestMoveInfo (best=$(tostring(bmi.bestmove)))")
    else
        print(
            io,
            "BestMoveInfo (best=$(tostring(bmi.bestmove)), ponder=$(tostring(bmi.ponder)))",
        )
    end
end

"""
    parsebestmove(line::String)::BestMoveInfo

Parses a `bestmove` line printed by a UCI engine to a `BestMoveInfo` object.
"""
function parsebestmove(line::String)::BestMoveInfo
    @assert startswith(line, "bestmove")
    tokens = split(line, r"\s+")
    bestmove = movefromstring(String(tokens[2]))
    ponder = if length(tokens) >= 4 && tokens[3] == "ponder"
        movefromstring(String(tokens[4]))
    else
        nothing
    end
    BestMoveInfo(bestmove, ponder)
end

"""
    SearchInfo

A struct representing the contents of a UCI engine's `info` output.

Contains the following slots, all of which can be `nothing` for a given line
of search output:

- `depth`: The current search depth.
- `seldepth`: The current selective search depth.
- `time`: The time spent searching so far, in milliseconds.
- `nodes`: The number of nodes searched so far.
- `pv`: The main line, as a vector of `Move` values.
- `multipv`: The multipv index of the line currently printed.
- `score`: The score, a value of type `Score`.
- `currmove`: The move currently searched, a value of type `Move`.
- `currmovenumber`: The index of the move currently searched in the move list.
- `hashfull`: A number in the range 0--100, indicating the transposition table
  saturation percentage.
- `nps`: Nodes/second count.
- `tbhits`: Number of tablebase hits.
- `cpuload`: CPU load percentage.
- `string`: An arbitrary string sent by the engine.
"""
mutable struct SearchInfo
    depth::Union{Nothing,Int}
    seldepth::Union{Nothing,Int}
    time::Union{Nothing,Int}
    nodes::Union{Nothing,Int}
    pv::Union{Nothing,Vector{Move}}
    multipv::Union{Nothing,Int}
    score::Union{Nothing,Score}
    currmove::Union{Nothing,Move}
    currmovenumber::Union{Nothing,Int}
    hashfull::Union{Nothing,Int}
    nps::Union{Nothing,Int}
    tbhits::Union{Nothing,Int}
    cpuload::Union{Nothing,Int}
    string::Union{Nothing,String}
end

function SearchInfo()
    SearchInfo(
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
    )
end

function Base.show(io::IO, si::SearchInfo)
    println(io, "SearchInfo:")
    if !isnothing(si.depth)
        println(io, " depth: $(si.depth)")
    end
    if !isnothing(si.seldepth)
        println(io, " seldepth: $(si.seldepth)")
    end
    if !isnothing(si.time)
        println(io, " time: $(si.time)")
    end
    if !isnothing(si.nodes)
        println(io, " nodes: $(si.nodes)")
    end
    if !isnothing(si.nps)
        println(io, " nps: $(si.nps)")
    end
    if !isnothing(si.score)
        println(io, " score: $(si.score)")
    end
    if !isnothing(si.currmove)
        println(io, " currmove: $(tostring(si.currmove))")
    end
    if !isnothing(si.currmovenumber)
        println(io, " currmovenumber: $(si.currmovenumber)")
    end
    if !isnothing(si.hashfull)
        println(io, " hashfull: $(si.hashfull)")
    end
    if !isnothing(si.tbhits)
        println(io, " tbhits: $(si.tbhits)")
    end
    if !isnothing(si.cpuload)
        println(io, " cpuload: $(si.cpuload)")
    end
    if !isnothing(si.string)
        println(io, " string: $(si.string)")
    end
    if !isnothing(si.multipv)
        println(io, " multipv: $(si.multipv)")
    end
    if !isnothing(si.pv)
        print(io, " pv:")
        for m ∈ si.pv
            print(io, " $(tostring(m))")
        end
        print(io, "\n")
    end
end

const TokenList = Vector{SubString{String}}

function parseinfoint(tokens::TokenList)::Tuple{Int,TokenList}
    (parse(Int, tokens[2]), tokens[3:end])
end

function parseinfocurrmove(tokens::TokenList)::Tuple{Move,TokenList}
    (movefromstring(String(tokens[2])), tokens[3:end])
end

function parseinfoscore(tokens::TokenList)::Tuple{Score,TokenList}
    ismate = tokens[2] == "mate"
    value = parse(Int, tokens[3])
    bound = exact
    if length(tokens) >= 4
        if tokens[4] == "lowerbound"
            bound = lower
        elseif tokens[4] == "upperbound"
            bound = upper
        end
    end
    (Score(value, ismate, bound), bound == exact ? tokens[4:end] : tokens[5:end])
end

function parseinfopv(tokens::TokenList)::Vector{Move}
    result = Move[]
    for t in tokens
        if !isempty(t)
            push!(result, movefromstring(String(t)))
        end
    end
    result
end

"""
    parsesearchinfo(line::String)::SearchInfo

Parses an `info` line printed by a UCI engine to a `SearchInfo` object.

See the documentation for `SearchInfo` for information about how to inspect
and use the return value.
"""
function parsesearchinfo(line::String)::SearchInfo
    @assert startswith(line, "info")
    result = SearchInfo()
    tokens = split(line, r"\s+")[2:end]
    while !isempty(tokens)
        if tokens[1] == "depth"
            (result.depth, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "seldepth"
            (result.seldepth, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "time"
            (result.time, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "nodes"
            (result.nodes, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "multipv"
            (result.multipv, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "currmove"
            (result.currmove, tokens) = parseinfocurrmove(tokens)
        elseif tokens[1] == "currmovenumber"
            (result.currmovenumber, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "hashfull"
            (result.hashfull, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "nps"
            (result.nps, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "tbhits"
            (result.tbhits, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "cpuload"
            (result.cpuload, tokens) = parseinfoint(tokens)
        elseif tokens[1] == "score"
            (result.score, tokens) = parseinfoscore(tokens)
        elseif tokens[1] == "pv"
            result.pv = parseinfopv(tokens[2:end])
            break
        else
            tokens = tokens[2:end]
        end
    end
    result
end

"""
    function search(e::Engine, gocmd::String; infoaction = nothing)

Tells a UCI engine to start searching.

The parameter `gocmd` is the actual command you want to send to the engine;
e.g. `"go movetime 10000"` or `"go infinite"`. The named parameter `infoaction`
is a function accepting the output of the engine's `"info"` commands and doing
something with the output. Usually, it will be some function making internal use
of `parsesearchinfo()`.

The return value is of type `BestMoveInfo`.
"""
function search(e::Engine, gocmd::String; infoaction = nothing)::BestMoveInfo
    sendcommand(e, gocmd)
    while true
        line = readline(e.io)
        if startswith(line, "bestmove")
            return parsebestmove(line)
        elseif !isnothing(infoaction) && startswith(line, "info")
            infoaction(line)
        end
    end
end

"""
    setboard(e::Engine, b::Board)
    setboard(e::Engine, g::SimpleGame)
    setboard(e::Engine, g::Game)

Set the engine's current board position to the given board/game state.
"""
function setboard(e::Engine, b::Board)
    sendcommand(e, touci(b))
end

function setboard(e::Engine, g::SimpleGame)
    sendcommand(e, touci(g))
end

function setboard(e::Engine, g::Game)
    sendcommand(e, touci(g))
end

"""
    touci(b::Board)
    touci(g::SimpleGame)
    touci(g::Game)

Create a UCI string representation of a board or a game.

# Examples
```julia-repl
julia> touci(startboard())
"position fen rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq -"

julia> sg = SimpleGame(); domove!(sg, "e4"); domove!(sg, "c5"); domove!(sg, "Nf3");

julia> touci(sg)
"position fen rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - moves e2e4 c7c5 g1f3"

julia> g = Game(); domove!(g, "d4"); domove!(g, "Nf6"); domove!(g, "c4");

julia> touci(g)
"position fen rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - moves d2d4 g8f6 c2c4"
```
"""
function touci(b::Board)::String
    "position fen $(fen(b))"
end

function touci(g::SimpleGame)::String
    result = "position fen $(fen(g.startboard))"
    if !isatbeginning(g)
        result *= " moves"
        for ply = 1:(g.ply-1)
            result *= " " * tostring(g.history[ply].move)
        end
    end
    result
end

function touci(g::Game)::String
    result = "position fen $(fen(g.root.board))"
    if !isatbeginning(g)
        result *= " moves"
        moves = String[]
        n = g.node
        while !isnothing(n.parent)
            push!(moves, tostring(lastmove(n.board)))
            n = n.parent
        end
        for m in reverse(moves)
            result *= " " * m
        end
    end
    result
end

"""
    mpvsearch(game, engine; nodes, depth, pvs)::Vector{SearchInfo}

Performs a multi-PV search and returns the result as a vector of `SearchInfo`.

Parameters:

- `game`: A `Game`, a `SimpleGame` or a `Board`.
- `engine`: An `Engine`.
- `nodes`: A named parameter instructing the engine to search to the desired
  tree size.
- `depth`: A named parameter instructing the engine to search to the given
  depth.
- `pvs`: The number of desired lines. Analysis for the `pvs` best moves is
  returned. If `pvs` is greater than the number of legal moves, analysis for
  all legal moves is returned.

At least one of `nodes` and `depth` must be supplied. If both are supplied,
the function will use `depth` and ignore `nodes`.

The function returns a vector of `SearchInfo` values, one for each of the `pvs`
best moves.
"""
function mpvsearch(
    g::Union{Board,SimpleGame,Game},
    e::Engine;
    nodes = nothing,
    depth = nothing,
    pvs = 100,
)::Vector{SearchInfo}

    result = SearchInfo[]

    function infoaction(info::String)
        info = parsesearchinfo(info)
        if !isnothing(info.multipv)
            info.multipv == 1 && empty!(result)
            push!(result, info)
        end
    end

    setoption(e, "MultiPV", pvs)
    setboard(e, g)
    if !isnothing(nodes)
        search(e, "go nodes $nodes", infoaction = infoaction)
    elseif !isnothing(depth)
        search(e, "go depth $depth", infoaction = infoaction)
    else
        Error("You must supply either the `nodes` or the `depth` parameter")
    end

    result
end

end


text:
using KernelAbstractions: Event, MultiEvent

using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceananigans.Architectures: device
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ExplicitFreeSurface, PrescribedVelocityFields

import Oceananigans.Utils: launch!

get_face(obj, face_index) = obj
get_face(t::Tuple, face_index) = Tuple(get_face(t_elem, face_index) for t_elem in t)
get_face(nt::NamedTuple, face_index) = NamedTuple{keys(nt)}(get_face(nt_elem, face_index) for nt_elem in nt)

get_face(grid::ConformalCubedSphereGrid, face_index) = grid.faces[face_index]
get_face(faces::CubedSphereFaces, face_index) = faces[face_index]

get_face(free_surface::ExplicitFreeSurface, face_index) =
    ExplicitFreeSurface(get_face(free_surface.η, face_index), free_surface.gravitational_acceleration)

get_face(velocities::PrescribedVelocityFields, face_index) =
    PrescribedVelocityFields(get_face(velocities.u, face_index),
                             get_face(velocities.v, face_index),
                             get_face(velocities.w, face_index),
                             velocities.parameters)

function get_face(op::KernelFunctionOperation, face_index)
    LX, LY, LZ = location(op)
    computed_dependencies = get_face(op.computed_dependencies, face_index)
    parameters = get_face(op.parameters, face_index)
    face_grid = get_face(op.grid, face_index)
    return KernelFunctionOperation{LX, LY, LZ}(op.kernel_function,
                                               computed_dependencies,
                                               parameters,
                                               face_grid)
end

function launch!(arch, grid::ConformalCubedSphereGrid, dims, kernel!, args...; kwargs...)

    events = []

    for (face_index, face_grid) in enumerate(grid.faces)
        face_args = Tuple(get_face(arg, face_index) for arg in args)
        event = launch!(arch, face_grid, dims, kernel!, face_args...; kwargs...)
        push!(events, event)
    end

    events = filter(e -> e isa Event, events)

    return MultiEvent(Tuple(events))
end

@inline launch!(arch, grid::ConformalCubedSphereGrid, ::Val{dims}, args...; kwargs...) where dims = launch!(arch, grid, dims, args...; kwargs...)


text:
using Revise
using CairoMakie, ZigZagBoomerang, SparseArrays, LinearAlgebra
using Random

function zz_sticky_events(ts0, xs0, i)
    first = true
    y = []
    t = Float64[]
    for (s, x) in zip(ts0, xs0)
        if x[i] == 0.0  && first == true
            push!(y, x)
            push!(t, s)
            first = false
        elseif x[i] == 0.0 && first == false
            continue
        else
            first = true
        end
    end
    return t, y
end

Random.seed!(10)
function ϕ(x, i, μ)
    x[i] - μ[i]
end
κ = .5
n = 2
μ = zeros(n)
x0 = [-1.0, -0.35]
θ0 = [1.0, -1.0]

T = 10.0
c = 0.1*ones(n)
@time trace0, _ = ZigZagBoomerang.sspdmp(ϕ, 0.0, x0, θ0, T, c, ZigZag(sparse(1.0I,n,n), zeros(n)), [κ,κ], μ)
ts0, xs0 = splitpairs(trace0)
# AbstractPlotting.available_marker_symbols()

fig  = Figure(resolution = (1000, 500))
p1 = fig[1, 1] = ax1 = Axis(fig, title = "ZigZag")

co = [(:black, 0.5) for i in 1:length(xs0)-1]
#co[8] = co[4] = co[3] = (:red, 0.5)
#co[12] = co[4] = (:blue, 0.5)
fr1 = findall(x->x[1]==0, diff(xs0))
fr2 = findall(x->x[2]==0, diff(xs0))
co[fr1] .=  [(:red, 0.5)]
co[fr2] .=  [(:blue, 0.5)]

dupl(x) = [x[1:end-1]'; x[2:end]'][:]
segs = dupl(xs0)
#lines!(getindex.(xs0,1), getindex.(xs0,2), color=co, markersize=0.1)
linesegments!(first.(segs), last.(segs), linewidth=1.5,color=co)
arrows!(first.(xs0)[1:1],last.(xs0)[1:1],first.(xs0)[2:2].-first.(xs0)[1:1],last.(xs0)[2:2].-last.(xs0)[1:1], arrowsize = 10.0, lengthscale = 0.5)

t1, y1 = zz_sticky_events(ts0, xs0, 1)
t2, y2 = zz_sticky_events(ts0, xs0, 2)
scatter!(getindex.(y1,1), getindex.(y1,2),
    color = (:red, 0.4), strokecolor = (:red, 0.6), markersize = 15, marker = :star4)
scatter!(getindex.(y2,1), getindex.(y2,2),
    color = (:blue, 0.4), strokecolor = (:blue, 0.6), markersize = 15, marker = :star4)
scatter!([xs0[1][1]], [xs0[1][2]], color = (:black, 1.0), strokecolor = (:black, 1.0), markersize = 5.5)
#text!("x(0), y(0)", textsize = 0.25, position = (xs0[1][1]- 0.05, xs0[1][2] + 0.05))
p1.xlabel = "x"
p1.ylabel = "y"
hidespines!(p1)
#p1.aspect = DataAspect()

p2 = fig[1,2] = ax2 = Axis(fig, title = "Coordinates")
linkyaxes!(ax1, ax2)

co1 = [(:black, 0.5) for i in 1:length(xs0)-1]
co2 = [(:black, 0.5) for i in 1:length(xs0)-1]

fr1 = findall(x->x[1]==0, diff(xs0))
fr2 = findall(x->x[2]==0, diff(xs0))
co1[fr1] .= [(:red, 0.5)]
co2[fr2] .= [(:blue, 0.5)]
#p2.aspect = DataAspect()

times = dupl(ts0)
linesegments!(times, getindex.(segs,1), color = co1, label = "x1", leg = true)
linesegments!(times, getindex.(segs,2), color = co2, linestyle = :dash, label = "x2")

scatter!(t1, getindex.(y1,1),
    color = (:red, 0.4), strokecolor = (:red, 0.6), markersize = 15, marker = :star4)
scatter!(t2, getindex.(y2,2),
    color = (:blue, 0.4), strokecolor = (:blue, 0.6), markersize = 15, marker = :star4)
hidespines!(p2)
p2.xlabel = "t"
p2.ylabel = "x,y"
scatter!(ts0[1:1], first.(xs0[1:1]),
    color = (:red, 0.4), strokecolor = (:red, 1.0), markersize = 5.5)
scatter!(ts0[1:1], getindex.(xs0[1:1],2),
    color = (:blue, 0.4), strokecolor = (:blue, 1.0), markersize = 5.5)

save("./figures/2d_sticky_samplers.png", fig)
fig


text:
using ModelingToolkit, DiffEqBase, LinearAlgebra

# Define some variables
@parameters t x y
@variables u(..)

ConstrainedEquation([x ~ 0, y < 1 / 2], u(t, x, y) ~ x + y^2)


text:
## Metsa-Oy Forest Products Supply Chain
## Edward J. Xu <edxu96@outlook.com>
## Jan 24th, 2020

"Get parameters for model 3."
function get_para()
  alpha_t = [190, 150, 120, 180, 150, 150]
  beta_t = [1, 0.5, 3, 0.2, 0.3, 0.2]
  a_i = [2, 2, 2, 2.8, 2.8]
  b_i = [0.8, 0.8, 0.8, 1.6, 1.6]
  e_i = [0.2, 0.2, 0.2, 0.2, 0.2]
  c_i = [550, 500, 450, 2500, 2600]
  r_saw = 200
  r_plywood = 90
  a_j = [4.8, 4.2]
  c_j = [820, 800]
  r_j = [220, 180]
  a_paper_t = [0, 1, 0]
  b_paper_j = [0.2, 0.2]
  c_paper = 1700
  r_paper = 80
  gamma_ik = [
    1600 1300 1400 1500;
    1400 1200 1300 1400;
    1300 1400 1500 1600;
    4400 3800 3600 3500;
    4300 4100 3900 3800
    ]
  delta_ik = [
    4 5 12 15;
    2 10 12 15;
    14 20 24 25;
    4 10 12 18;
    4 7 12 15
    ]
  gamma_jk = [
    2300 2500 2300 2600;
    2500 2800 2300 2500
    ]
  delta_jk = [
    1 4 5 6;
    3 2 6 7
    ]
  gamma_paper_k = [4500, 4700, 4300, 4800]
  delta_paper_k = [7, 10, 12, 15]
  p_fuel = 40

  i_t1 = [[1], [2, 4], [3, 5]]
  i_t2 = [[1], [2, 4], [3, 5]]
  j_t2 = [[1], [], [2]]

  ## Extra parameters in model 2
  sigma = 0.95
  omega_i = [1.010, 1.008, 1.015, 1.015, 1.020]
  omega_j = [1.025, 1.030]
  omega_paper = 1.035
  nu_saw = 1.5
  nu_plywood = 1.5
  nu_j = [2, 2]
  nu_paper = 2
  o_saw = 100
  o_plywood = 300
  o_j = [500, 500]
  o_paper = 700

  ## Extra parameters in model 3
  pi_n = [0.25, 0.25, 0.25, 0.25]
  rho_nm = [
    1 1.05 1.07;
    1 1.05 0.95;
    1 0.95 1.05;
    1 0.95 0.93
    ]

  return alpha_t, beta_t, a_i, b_i, e_i, c_i, r_saw, r_plywood, a_j, c_j,
    r_j, a_paper_t, b_paper_j, c_paper, r_paper, gamma_ik, delta_ik, gamma_jk,
    delta_jk, gamma_paper_k, delta_paper_k, p_fuel, sigma, omega_i, omega_j,
    omega_paper, nu_saw, nu_plywood, nu_j, nu_paper, o_saw, o_plywood, o_j,
    o_paper, i_t1, i_t2, j_t2, pi_n, rho_nm
end


text:
module Collatz

using Base.Checked

export collatz

function collatz(n::Integer)
    k = 0
    while n > 1
        n = iseven(n) ? n >> 1 : checked_add(checked_mul(3,n), 1)
        k += 1
    end
    return k
end

end


text:
using Nemo, LinearAlgebra, DelimitedFiles

# using Revise
using SymbolicCRN

M = [63   90   36  -45  -18
     -7  -10   -4    5    2
     28   40   16  -20   -8
     -56  -80  -32   40   16
     7   10    4   -5   -2]

K = [2  5  0  0  0  0
     0  0  1  0  0  1
     0  0  0  5  1  0
     0  7  2  4  0  0
     7  0  0  0  2  5]

R, vrs = PolynomialRing(Nemo.ZZ, vcat(["k$i" for i in 1:5], ["x$i" for i in 1:4]));
F = [vrs[1] * vrs[2] * vrs[6] - vrs[8]; vrs[3] * vrs[9] + 2 * vrs[7]];

@test Jacobian(R, F, 6:9) == [vrs[1] * vrs[2]  0  -1   0
                              0  2   0  vrs[3]]
@test Jacobian(R, F) == [vrs[2] * vrs[6]  vrs[1] * vrs[6]   0  0  0  vrs[1] * vrs[2]  0  -1   0
                         0      0  vrs[9]  0  0      0  2   0  vrs[3]]


text:
using MatrixCompletion

@testset "$(format("Tracker: IndexTracker[convert]"))" begin
  let
    data = [1, 2, 3, 4]
    tc = convert(Array{<:CartesianIndex}, data)
    @test typeof(tc) == Array{CartesianIndex{1}, 1}
  end
  let
    data = rand(5, 5)
    tc = convert(Array{<:CartesianIndex}, findall(x -> x < 0.5 , data))
    @test typeof(tc) == Array{CartesianIndex{2}, 1}
  end
end

@testset "$(format("Tracker: IndexTracker[data stream initialization]"))" begin
  let
    data1 = [:Gaussian :Bernoulli;
             :Gaussian :Bernoulli]
    data2 = [:Observed :Observed;
             :Missing  :Missing]
    tc = IndexTracker{Symbol}(data1, data2)
  end
  let
    data1 = [:Gaussian :Bernoulli;
             :Gaussian :Bernoulli]
    data2 = [:Observed :Observed :Observed;
             :Missing  :Missing :Missing]
    @test_throws DimensionMismatch  tc = IndexTracker{Symbol}(data1, data2)
  end
end

@testset "$(format("Tracker: IndexTracker[disjoint join]"))" begin
  let
    tc = IndexTracker{Symbol}([:a, :b])
    @test_throws DimensionMismatch disjoint_join(tc, [:a])
    @test_throws MethodError disjoint_join(tc, [:a, :c])
  end
  let
    tc = IndexTracker{Symbol}([:a :b; :c :d])
  end
end

@testset "$(format("Tracker: IndexTracker[getindex]"))" begin
  let
    data = [:a,:b]
    tc = IndexTracker{Symbol}(data)
    @test data[tc[:a]] == data[findall(x -> x == :a, data)]
  end
  let
    data = [:a :b;
            :c :d]
    tc = IndexTracker{Symbol}(data)
    @test data[tc[:a]] == data[findall(x -> x == :a, data)]
    @test data[tc[:b]] == data[findall(x -> x == :b, data)]
    @test data[tc[:c]] == data[findall(x -> x == :c, data)]
    @test data[tc[:d]] == data[findall(x -> x == :d, data)]
  end
end

@testset "$(format("Tracker: IndexTracker[symbol]"))" begin
  let
    tc1 = IndexTracker{Symbol}([:a, :b])
    @test typeof(tc1) == IndexTracker{Symbol}
    @test size(tc1) == (2, )
    @test tc1.dimension == (2, )
  end
end

@testset "$(format("Tracker: IndexTracker[groupby]"))" begin
  let
    data = Array{Symbol}(undef, 20, 20)
    data[:, 1:10]  .= :Gaussian
    data[:, 11:20] .= :Binomial
    data2 = Array{Symbol}(undef, 20, 20)
    data2[:, 1:5]   .= :Observed
    data2[:, 10:15] .= :Observed
    data2[:, 6:10]  .= :Missing
    data2[:, 16:20] .= :Missing
    tc = IndexTracker{Symbol}(data)
    # @show(size(Iterators.flatten(data[:, 1:10])))
    @test tc[:Gaussian] == findall(x -> x == :Gaussian, data)
    @test tc[:Binomial] == findall(x -> x == :Binomial, data)
    disjoint_join(tc, data2)
    @test tc[:Observed] == findall(x -> x == :Observed, data2)
    @test tc[:Missing] == findall(x -> x == :Missing, data2)
    tc2 = groupby(tc, [:Observed, :Missing])
    @test tc2[:Gaussian][:Observed] == intersect(findall(x -> x == :Gaussian, data),
                                                 findall(x -> x == :Observed, data2))
    @test tc2[:Gaussian][:Missing] == intersect(findall(x -> x == :Gaussian, data),
                                                findall(x -> x == :Missing, data2))
    @test tc2[:Binomial][:Observed] == intersect(findall(x -> x == :Binomial, data),
                                                 findall(x -> x == :Observed, data2))
    @test tc2[:Binomial][:Missing] == intersect(findall(x -> x == :Binomial, data),
                                                findall(x -> x == :Missing, data2))
  end
end


text:
struct ImageTexture{T} <: Texture where T <: Real
    width::Integer
    height::Integer
    img::Array{RGBA{T}, 2}

    function ImageTexture(img::Array{RGBA{T}, 2}) where T <: Real
        img_size::Tuple{Integer, Integer} = size(img)
        new{T}(img_size[2], img_size[1], img)
    end
end

function value(it::ImageTexture, u::AbstractFloat, v::AbstractFloat, point::AbstractVec)::AbstractVec
    i::Integer = convert(Integer, floor(u * it.width)) + 1
    j::Integer = convert(Integer, floor((1.0 - v) * it.height)) + 1
    i = max(1, i)
    j = max(1, j)
    i = min(it.width, i)
    j = min(it.height, j)
    pixel::RGBA = it.img[j, i]
    r::AbstractFloat = pixel.r
    g::AbstractFloat = pixel.g
    b::AbstractFloat = pixel.b
    a::AbstractFloat = pixel.alpha
    Vec(r * a, g * a, b * a)
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
################################################################################
#
#  Two element generated ideals
#
################################################################################

@doc Markdown.doc"""
    has_2_elem(A::NfAbsOrdIdl) -> Bool

Returns whether $A$ is generated by two elements.
"""
function has_2_elem(A::NfAbsOrdIdl)
  return isdefined(A, :gen_two)
end

@doc Markdown.doc"""
    has_weakly_normal(A::NfAbsOrdIdl) -> Bool

Returns whether $A$ has weakly normal two element generators.
"""
function has_weakly_normal(A::NfAbsOrdIdl)
  return (isdefined(A, :gens_weakly_normal) &&
        A.gens_weakly_normal == true) || has_2_elem_normal(A)
end

@doc Markdown.doc"""
    has_2_elem_normal(A::NfAbsOrdIdl) -> Bool

Returns whether $A$ has normal two element generators.
"""
function has_2_elem_normal(A::NfAbsOrdIdl)
  #the one ideal <1, ?> is automatomatically normal>
  return isdefined(A, :gens_normal) && (A.gen_one == 1 || A.gens_normal > 1)
end

################################################################################
#
#  Predicates
#
################################################################################

function isconsistent(A::NfAbsOrdIdl)
  if has_2_elem_normal(A) && !iszero(A.gen_two)
    if !defines_2_normal(A)
      return false
    end
  end
  if has_norm(A)
    b = basis_mat(A, copy = false)
    if det(b) != A.norm
      return false
    end
  end
  return true
end

# check if gen_one,gen_two is a P(gen_one)-normal presentation
# see Pohst-Zassenhaus p. 404
function defines_2_normal(A::NfAbsOrdIdl)
  m = A.gen_one
  gen = A.gen_two
  mg = denominator(inv(gen.elem_in_nf), order(A))
  # the minimum of ideal generated by g
  g = gcd(m,mg)
  return gcd(m, div(mg, g)) == 1
end

###########################################################################################
#
#  2-element normal presentation
#
###########################################################################################

# The following makes sure that A has a weakly normal presentation
# Recall that (x,y) are a weakly normal presentation for A
# if and only if norm(A) = gcd(norm(x), norm(y))
#
# Maybe we should allow an optional paramter (an fmpz),
# which should be the first generator.
# So far, the algorithm just samples (lifts of) random elements of A/m^2,
# where m is the minimum of A.

function _assure_weakly_normal_presentation(A::NfAbsOrdIdl)
  if has_2_elem(A) && has_weakly_normal(A)
    return nothing
  end

  if isdefined(A, :princ_gen)
    x = A.princ_gen
    b = x.elem_in_nf

    bi = inv(b)

    A.gen_one = denominator(bi, order(A))
    A.minimum = A.gen_one
    A.gen_two = x
    if !isdefined(A, :norm)
      A.norm = abs(_normmod(A.gen_one^degree(order(A)), b))
    end
    @hassert :NfOrd 1 gcd(A.gen_one^degree(order(A)),
                    FlintZZ(norm(A.gen_two))) == A.norm

    if A.gen_one == 1
      A.gens_normal = 2*A.gen_one
    else
      A.gens_normal = A.gen_one
    end
    A.gens_weakly_normal = 1
    return nothing
  end

  O = order(A)
  if minimum(A) == 0
    A.gen_one = minimum(A)
    A.gen_two = zero(O)
    A.gens_weakly_normal = true
    return nothing
  end

  if minimum(A) == 1
    A.gen_one = minimum(A)
    A.gen_two = one(O)
    A.gens_weakly_normal = true
    A.gens_normal = fmpz(2)
    return nothing
  end

  if !issimple(nf(order(A))) || !isdefining_polynomial_nice(nf(order(A)))
    B = Array{fmpz}(undef, degree(O))
    Amin2 = minimum(A)^2
    Amind = minimum(A)^degree(O)
    BB = basis(A, copy = false)
    r = -Amin2:Amin2

    while true
      rand!(B, r)
      gen = dot(B, BB)
      if norm(A) == gcd(Amind, norm(gen))
        A.gen_one = minimum(A)
        A.gen_two = O(gen, false)
        A.gens_weakly_normal = true
        return nothing
      end
    end
  end

  Amin2 = minimum(A)^2
  Amind = gcd(minimum(A)^degree(O), minimum(A)*norm(A))

  B = Array{fmpz}(undef, degree(O))

  gen = O()

  r = -Amin2:Amin2

  m = zero_matrix(FlintZZ, 1, degree(O))

  cnt = 0
  while true
    cnt += 1

    if cnt > 100 && is2_normal_difficult(A)
      assure_2_normal_difficult(A)
      return nothing
    end

    #if cnt > 1000
    #  println("Having a hard time find weak generators for $A")
    #end

    rand!(B, r)

    # Put the entries of B into the (1 x d)-Matrix m
    for i in 1:degree(O)
      s = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz}, (Ref{fmpz_mat}, Int, Int), m, 0, i - 1)
      ccall((:fmpz_set, :libflint), Nothing, (Ptr{fmpz}, Ref{fmpz}), s, B[i])
    end

    if iszero(m)
      continue
    end

    mul!(m, m, basis_mat(A, copy = false))
    d = denominator(basis_mat(O, copy = false))
    mul!(m, m, basis_mat(O, copy = false).num)
    gen = elem_from_mat_row(nf(O), m, 1, d)
    d = denominator(gen)
    f, e = ppio(d, minimum(A, copy = false))
    gen = mod(numerator(gen), f*minimum(A)^2)//f
    if iszero(gen)
      continue
    end

    # the following should be done inplace
    #gen = dot(reshape(Array(mm), degree(O)), basis(O))
    if norm(A) == _normmod(Amind, O(gen, false))#gcd(Amind, numerator(norm(gen)))
      A.gen_one = minimum(A)
      A.gen_two = O(gen, false)
      A.gens_weakly_normal = true
      return nothing
    end
  end
end
#
# Here some random guesses for the difficult 2-element generators
# degree  | d
#   < 7   | 1
#  8 - 12 | 2 * 3
# 13 - 20 | 2 * 3 * 5
#  >= 21  | 2 * 3 * 5 * 7

function is2_normal_difficult(A::NfAbsOrdIdl)
  d = 2
  m = minimum(A)
  ZK = order(A)
  n = degree(ZK)

  if n < 7
    return false
  end

  if n < 12 && isone(gcd(m, 2 * 3))
    return false
  end

  if n < 20 && isone(gcd(m, 2 * 3 * 5))
    return false
  end

  if isone(gcd(m, 2 * 3 * 5 * 7))
    return false
  end

  return true
end

function valuation_of_ideal_difficult(A::NfAbsOrdIdl, P::NfAbsOrdIdl)
  B = basis(A, copy = false)
  res = valuation(B[1], P)
  for i = 2:length(B)
    if iszero(res)
      break
    end
    res = min(res, valuation(B[i], P))
  end
  return res
end

function assure_2_normal_difficult(A::NfAbsOrdIdl)

  m = minimum(A)
  ZK = order(A)
  n = degree(ZK)

  if !is2_normal_difficult(A)
    assure_2_normal(A)
    return nothing
  end

  if n < 12
    d = Int[2, 3]
  elseif n < 20
    d = Int[2, 3, 5]
  else
    d = Int[2, 3, 5, 7]
  end

  I = ideal(ZK, 1)
  for i = 1:length(d)
    m1, m = ppio(m, fmpz(d[i]))
    if isone(m1)
      continue
    end
    A1 = gcd(A, m1)
    lp = prime_decomposition(ZK, d[i])
    v = Int[valuation_of_ideal_difficult(A1, p[1]) for p = lp]
    for i = 1:length(v)
      if v[i] > 0
        I *= lp[i][1]^v[i]
      end
    end
  end
  #Now the ideal I has a 2 normal presentation
  if isone(m)
    A.gen_one = I.gen_one
    A.gen_two = I.gen_two
    A.gens_normal = I.gens_normal
    A.gens_weakly_normal = I.gens_weakly_normal
    A.gens_short = I.gens_short
    return nothing
  end
  J = gcd(A, m)
  assure_2_normal(J)
  C = I * J
  A.gen_one = C.gen_one
  A.gen_two = C.gen_two
  A.gens_normal = C.gens_normal
  A.gens_weakly_normal = C.gens_weakly_normal
  A.gens_short = C.gens_short
  return nothing
end

################################################################################
#
#  Assure 2 normal
#
################################################################################

function assure_2_normal(A::NfAbsOrdIdl)
  if has_2_elem(A) && has_2_elem_normal(A)
    @hassert :NfOrd 1 isconsistent(A)
    return nothing
  end
  O = order(A)
  K = nf(O)
  n = degree(K)

  if norm(A) == 1
    A.gen_one = fmpz(1)
    A.gen_two = one(O)
    A.gens_normal = fmpz(1)
    return
  end

  if has_2_elem(A)
    m = minimum(A)
    bas = basis(O)
    # Magic constants
    if m > 1000
      r = -500:500
    else
      r = -div(Int(m)+1,2):div(Int(m)+1,2)
    end
    #gen = K()
    #s = K()
    gen = zero(O)
    s = O()
    cnt = 0
    while true
      cnt += 1
      if cnt > 100 && is2_normal_difficult(A)
        assure_2_normal_difficult(A)
        @hassert :NfOrd 1 isconsistent(A)
        return
      end
      if cnt > 1000
        error("Having a hard time making generators normal for $A")
      end
      rand!(s, O, r)
      mul!(s, s, A.gen_two)
      add!(gen, rand(r)*A.gen_one, gen)
      add!(gen, s, gen)
      gen = mod(gen, m^2)

      if iszero(gen)
        continue
      end

      mg = _minmod(m^2, gen)
      g = gcd(m, mg)
      if gcd(m, div(mg, g)) == 1
        if gcd(m^n, norm(gen)) != norm(A)
          @vprint :NfOrd 2 "\n\noffending ideal $A \ngen is $gen\nWrong ideal\n"
          cnt += 10
          continue
        end
        break
      end
    end
    @vprint :NfOrd 2 "used $cnt attempts\n"
    A.gen_one = m
    A.gen_two = gen
    A.gens_normal = m
    @hassert :NfOrd 1 isconsistent(A)
    return
  end

  m = minimum(A)
  bas = nf_elem[ elem_in_nf(b) for b in basis(A, copy = false)]
  # Magic constants
  if m > 1000
    r = -500:500
  else
    r = -div(Int(m)+1,2):div(Int(m)+1,2)
  end
  cnt = 0
  while true
    cnt += 1
    if cnt > 100 && is2_normal_difficult(A)
      assure_2_normal_difficult(A)
      @hassert :NfOrd 1 isconsistent(A)
      return
    end
    if cnt > 1000
      error("Having a hard time making generators normal for $A")
    end
    gen = O(rand(bas, r))
    gen = mod(gen, m^2)

    if iszero(gen)
      continue
    end

    mg = _minmod(m^2, gen)
    if gcd(m, div(mg, gcd(mg, m))) == 1
      if gcd(m^n, norm(gen)) != norm(A)
        @vprint :NfOrd 1 "\n\noffending ideal $A \ngen is $gen\nWrong ideal\n"
        cnt += 10
        continue
      end
      break
    end
  end
  @vprint :NfOrd 2 "used $cnt attempts\n"
  A.gen_one = m
  A.gen_two = gen
  A.gens_normal = m
  @hassert :NfOrd 1 isconsistent(A)
  return nothing
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Calcium_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Calcium")
JLLWrappers.@generate_main_file("Calcium", UUID("fcfa6d1b-d8ce-59d5-8c0a-c0d7f69e4f40"))
end  # module Calcium_jll


text:
#
# Example of a slightly larger graphene calculation.
#src tags: long
#

using DFTK
setup_threading()

kgrid = [12, 12, 4]
Ecut = 15

lattice = [4.66 -2.33 0.00;
           0.00  4.04 0.00
           0.00  0.00 15.12]
C = ElementPsp(:C, psp=load_psp("hgh/pbe/c-q4"))
atoms     = [C, C]
positions = [[0.0, 0.0, 0.0], [1//3, 2//3, 0.0]]

model = model_PBE(lattice, atoms, positions; temperature=1e-3, smearing=Smearing.Gaussian())
basis = PlaneWaveBasis(model; Ecut, kgrid)

# Run SCF
scfres = self_consistent_field(basis)

println("Carbon forces:")
println(compute_forces(scfres)[1])

# Print obtained energies
println()
println(scfres.energies)


text:
function mkf2e(t,elemtype,nd)

ne = size(t,2);

# get local faces of an element
face = getelemface(nd,elemtype);
nvf,nfe = size(face);

# sort faces for all elements to make two identical faces have the same numbering order
N = nfe*ne;
tf = reshape(t[face[:],:],nvf, N);
tf = sort(tf, dims=1);

tf, jx = sortcolumns(tf); # do this so that two identical faces are next each other
jx = Int.(jx);
dx = tf[:,2:end]-tf[:,1:end-1]; # do this to find two identical faces
dx = sum(dx.*dx,dims=1);           # identical faces corresponds to zero entries
in1 = findall(dx[:].==0);           # interior global faces connected to 1st elements (1st identical faces)
in2 = in1.+1;                 # interior global faces connected to 2nd elements (2nd identical faces)

in0 = setdiff(collect(1:length(jx)), unique([in1[:]; in2[:]])); # boundary global faces connected to 1st elements

nf = length(in0)+length(in1);
f2e = zeros(Int,4,nf);   # allocate memory
e2e = zeros(Int,nfe,ne); # allocate memory

# interior faces
e1 = Int.(ceil.(Float64.(jx[in1])./Float64(nfe)));      # 1st elements
l1 = jx[in1] .- (e1.-1).*nfe;   # 1st local faces
e2 = Int.(ceil.(Float64.(jx[in2])./Float64(nfe)));      # 2nd elements
l2 = jx[in2] .- (e2.-1).*nfe;   # 2nd local faces
g = 1:length(in1);           # indices for interior faces
f2e[1,g] = e1;
f2e[2,g] = l1;
f2e[3,g] = e2;
f2e[4,g] = l2;
for i = 1:length(e1)
    e2e[l1[i],e1[i]] = e2[i];
    e2e[l2[i],e2[i]] = e1[i];
end

# boundary faces
e1 = Int.(ceil.(Float64.(jx[in0])./Float64(nfe)));      # 1st elements
l1 = jx[in0] .- (e1.-1).*nfe;   # 1st local faces
g = (length(in1)+1):nf;      # indices for boundary faces
f2e[1,g] = e1;
f2e[2,g] = l1;

return f2e, e2e

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
export
    IterativeSampler

"""
    struct IterativeSampler <: IterativeAlgorithm

An iterative algorithm that uses a sampler to accumulate more samples on each refinement.
"""
struct IterativeSampler <: IterativeAlgorithm
    base_algorithm :: InstantAlgorithm
end

function prepare(alg::IterativeSampler, runtime::InstantRuntime,
    evidence::Dict{Symbol, <:Score} = Dict{Symbol, Score}(),
    interventions::Dict{Symbol, <:Dist} = Dict{Symbol, Dist}(),
    placeholder_beliefs = get_placeholder_beliefs(runtime))
    net = get_network(runtime)
    for (n,e) in evidence
        v = get_node(net, n)
        inst = current_instance(runtime, v)
        post_evidence!(runtime, inst, e)
    end
    for (n,i) in interventions
        v = get_node(net, n)
        inst = current_instance(runtime, v)
        post_intervention!(runtime, inst, i)
    end
    for (n,b) in placeholder_beliefs
        p = get_node(net, n)
        inst = current_instance(runtime, p)
        post_belief!(runtime, inst, b)
    end
    set_state!(runtime, :particles, Particles(Sample[], Float64[]))
end

function refine(alg::IterativeSampler, runtime::InstantRuntime)
    current_particles = get_state(runtime, :particles)
    infer(alg.base_algorithm, runtime)
    new_particles = get_state(runtime, :particles)
    all_samples = copy(current_particles.samples)
    append!(all_samples, new_particles.samples)
    all_log_weights = copy(current_particles.log_weights)
    append!(all_log_weights, new_particles.log_weights)
    set_state!(runtime, :particles, Particles(all_samples, all_log_weights))
end

answer(q::Query, alg::IterativeSampler, run::Runtime, inst::VariableInstance) = answer(q, alg.base_algorithm, run, inst)

answer(q::Query, alg::IterativeSampler, run::Runtime, insts::Vector{VariableInstance}) = answer(q, alg.base_algorithm, run, insts)


text:
module SampleTissueMap

using CSV
using DataFrames

function getTissueMap(file)
    #df = CSV.read(file, delim="\t", weakrefstrings=false, nullable=false)
    df = CSV.read(file, delim="\t")
    sampleTissueMap = Dict()
    for row in eachrow(df)
        sampleTissueMap[Symbol(row[Symbol("sample")])] = row[Symbol("tissue_type")]
    end

    return sampleTissueMap
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath