text:
using Plots: plot3d, plot, Animation, frame, gif, GridLayout, EmptyLayout, create_grid
using StaticArrays

"""
    Body(m, x, y, z, v)

Create a celestial body object with mass `m` kilograms, position `x` `y` `z`,
and 3-component velocity vector `v` (m/s).
"""
mutable struct Body{T <: Float64}
	# The mass of the body
	m::T

	# Its coords
	x::T
	y::T
	z::T

	# Its velocity as a vector of x y z values
	v::MVector{3, T}
end

"""
    getforce(b1, b2) -> float

Return the magnitude of the force between Body objects `b1` and `b2`.

Uses the standard gravitational equation.

See also: [`getforcevector`](@ref)
"""
function getforce(b1::Body, b2::Body)::Float64
	a = b1.x - b2.x
	b = b1.y - b2.y
	c = b1.z - b2.z

	r2 = a^2 + b^2 + c^2
	if r2 == 0; error("Cannot divide by 0. The positions of two bodies cannot be the same."); end

	# G \frac{m_1 m_2}{r^2}
	(6.674e-11 * b1.m * b2.m) / r2
end

"""
    getforcevector(b1, b2) -> [float, float, float]

Return the 3D vector of the force from Body object `b1` towards `b2`.

Expressed as a 3-component vector. Magnitude can be computed
by `getforce(b1, b2)`. Angles must be computed separately.

See also: [`getforce`](@ref)
"""
function getforcevector(b1::Body, b2::Body)::MVector{3, Float64}
	# azimuth is the angle on the x-y plane
	Δx = b2.x - b1.x
	Δy = b2.y - b1.y

	azimuth = atan(Δy, Δx)

	# altitude is the angle from the x-y plane in the z axis
	hyp = sqrt(Δx^2 + Δy^2)
	Δz = b2.z - b1.z
	altitude = atan(Δz, hyp)

	F = getforce(b1, b2)

	# We use F to compute the z component and the new hypotenuse to use
	# for computing the x and y components
	vec_z = F * sin(altitude)
	xy_plane_hyp = F * cos(altitude)

	# We then do some more trig to compute the x and y components
	vec_y = xy_plane_hyp * sin(azimuth)
	vec_x = xy_plane_hyp * cos(azimuth)

	MVector{3, Float64}(vec_x, vec_y, vec_z)
end

"""
    step!(bodies, Δt)

Apply forces between all Body objects in `bodies` over time step `Δt` seconds.
"""
function step!(bodies::Vector{Body}, Δt::Float64)
	# Calculate the force on each body and update its velocity accordingly
	for i in 1:length(bodies)
		b = bodies[i]

		forceonbody = copy(b.v)
		# We loop over every body and if it's not the same one, we add the add the force
		for j in 1:length(bodies)
			if j != i
				# We add the body's velocity to the force applied by gravity
				# to get the resultant force
				forceonbody .+= getforcevector(b, bodies[j])
			end

		end

		# Update the velocity of the body

		# a = (v - u)/t and F = ma ⟹  F = m * (v - u)/t
		# ⟹  v = Ft/m + u

		# Since u is the initial value, we can just use +=
		# We actually need .+= here to broadcast over the MVector
		b.v .+= ((forceonbody * Δt) / b.m)
	end

	# We need to wait until we've calculated all the forces before we start to change the positions,
	# otherwise we would change a position and then calculate a force based on the new position,
	# rather than the old one
	for i in 1:length(bodies)
		b = bodies[i]

		# Then we update the position of the body using its velocity
		b.x += Δt * b.v[1]
		b.y += Δt * b.v[2]
		b.z += Δt * b.v[3]
	end
end

"""
    drawframes(bodies, framecount, Δt)

Return frame data for the list of Body objects `bodies` with time step `Δt`.

This frame data is a Vector{Vector{SVector{3, Float64}}}. It can be thought of as a list of frames,
where each frame has a list of bodies, and each body has a list of x y z coordinates.

See also: [`drawgif`](@ref), [`creategif`](@ref)
"""
function drawframes(bodies::Vector{Body}, framecount::Int, Δt::Float64)::Vector{Vector{SVector{3, Float64}}}
	frames::Vector{Vector{SVector{3, Float64}}} = [[SVector{3, Float64}(bodies[i].x, bodies[i].y, bodies[i].z) for i in 1:length(bodies)]]

	for _ in 1:framecount
		step!(bodies, Δt)

		# We build the position data for this frame
		thisframe::Vector{SVector{3, Float64}} = [SVector{3, Float64}(bodies[i].x, bodies[i].y, bodies[i].z) for i in 1:length(bodies)]

		push!(frames, thisframe)
	end

	return frames
end

"""
    drawgif(positions, cube=false, bounds=nothing, filename=nothing)

Draw a GIF using the position data.

`positions` is a Vector{Vector{SVector{3, Float64}}}. It can be thought of as a list of frames,
where each frame has a list of bodies, and each body has a list of x y z coordinates.

`cube` is a bool for whether to draw the plot bounds as a cube.

`bounds` is a Vector{Tuple{Float64, Float64}} - a list of xlimits, ylimits, and zlimits.
If it's nothing (by default), then the bounds will be auto-generated.

`filename` is an optional filename. If none is provided, one will be programmatically generated,
following the pattern of `out.gif`, `out_1.gif`, `out_2.gif`, etc.

See also: [`drawframes`](@ref), [`creategif`](@ref)
"""
function drawgif(positions::Vector{Vector{SVector{3, Float64}}}, cube::Bool, bounds::Union{Nothing, Vector{Tuple{Float64, Float64}}}=nothing, filename::Union{String, Nothing}=nothing)
	# This is a multiplier to make the bounding box just a bit bigger than strictly necessary,
	# just to make it look a bit nicer
	multiplier = 1.05

	if cube && !isnothing(bounds)
		# We get the maximum absolute value in the bounds and use that to build the cube
		maxbound = multiplier * max(abs.(Iterators.flatten(bounds))...)

		xlimits = (-1 * maxbound, maxbound)
		ylimits = (-1 * maxbound, maxbound)
		zlimits = (-1 * maxbound, maxbound)

	elseif !cube && !isnothing(bounds)
		# We just use the bounds given, multiplied by the multiplier
		xlimits = multiplier .* bounds[1]
		ylimits = multiplier .* bounds[2]
		zlimits = multiplier .* bounds[3]

	elseif cube && isnothing(bounds)
		# This line gets absolute value of every float in positions and then gets the max
		# Because we're dealing with a nested vector, we have to flatten it, and then
		# flatten the result, and then splat it for max
		maxbound = multiplier * max(abs.(Iterators.flatten(Iterators.flatten(positions)))...)

		xlimits = (-1 * maxbound, maxbound)
		ylimits = (-1 * maxbound, maxbound)
		zlimits = (-1 * maxbound, maxbound)

	# Not cube or explicit bounds
	else
		# These list comprehensions loop through the list to get the coords
		# i is every frame in the positions list, and j is every body in each frame
		# We only need to check length(positions[1]) to get j because the length of
		# every sub-list is the same because every frame has the same number of bodies
		xs = [positions[i][j][1] for i in 1:length(positions) for j in 1:length(positions[1])]
		ys = [positions[i][j][2] for i in 1:length(positions) for j in 1:length(positions[1])]
		zs = [positions[i][j][3] for i in 1:length(positions) for j in 1:length(positions[1])]

		xlimits = multiplier .* extrema(xs)
		ylimits = multiplier .* extrema(ys)
		zlimits = multiplier .* extrema(zs)
	end

	# This is just the number of bodies
	n = length(positions[1])

	# We have 4 different subplots. xyz is the 3D perspective camera angle. The others are the orthogonal views.
	# We then put these subplots into a full plot to show everything comprehensively
	xyz = plot3d(n; legend=false, lw=1, title="$n Body Gravity Simulation", xlim=xlimits, ylim=ylimits, zlim=zlimits)
	# We're using standard 2D line plots here to avoid issues with cameras in 3D space
	xy = plot(n; legend=false, lw=0.5, title="Plan", titlefontsize=8, xlim=xlimits, ylim=ylimits, tickfontsize=4)
	xz = plot(n; legend=false, lw=0.5, title="Front", titlefontsize=8, xlim=xlimits, ylim=zlimits, tickfontsize=4)
	yz = plot(n; legend=false, lw=0.5, title="Side", titlefontsize=8, xlim=ylimits, ylim=zlimits, tickfontsize=4)

	# This is the expanded form of the @gif macro over a for loop
	# We're using the expanded form rather than the macro itself because
	# that lets us control the filename
	anim = Animation()
	counter = 1
	# For each frame
	for i = 1:length(positions)
		# For each body in this frame
		for j in 1:length(positions[i])
			# We add the position data of every body to its respective series in the plots
			bodydata = positions[i][j]
			push!(xyz, j, bodydata[1], bodydata[2], bodydata[3])
			push!(xy, j, bodydata[1], bodydata[2])
			push!(xz, j, bodydata[1], bodydata[3])
			push!(yz, j, bodydata[2], bodydata[3])
		end

		# Only if this counter is a multiple of 10, do we add this frame to the animation
		if counter % 10 == 0
			# I can't directly import the @layout macro, so we do it this way
			layout = eval(create_grid(:([a{0.7h}; b c d])))
			frame(anim, plot(xyz, xy, xz, yz; layout = layout, size = (1000, 1000)))
		end

		counter += 1
	end

	if isnothing(filename)
		if !in("out.gif", readdir())
			filename = "out.gif"
		else
			i = 1
			while in("out_$i.gif", readdir())
				i += 1
			end
			filename = "out_$i.gif"
		end
	end

	gif(anim, filename)
end

"""
    creategif(bodies, framecount, Δt, cube=false, bounds=nothing, filename=nothing)

Generate a GIF from the list of Body objects `bodies`, with `framecount` frames and a time step of `Δt` seconds.

`cube` is a bool for whether to draw the plot bounds as a cube.

`bounds` is a Vector{Tuple{Float64, Float64}} - a list of xlimits, ylimits, and zlimits.
If it's nothing (by default), then the plot bounds will be auto-generated.

`filename` is an optional filename. If none is provided, one will be programmatically generated,
following the pattern of `out.gif`, `out_1.gif`, `out_2.gif`, etc.

See also: [`drawgif`](@ref), [`drawframes`](@ref)
"""
creategif(bodies::Vector{Body}, framecount::Int, Δt::Float64, cube::Bool=false, bounds::Union{Nothing, Vector{Tuple{Float64, Float64}}}=nothing, filename::Union{String, Nothing}=nothing) = drawgif(drawframes(bodies, framecount, Δt), cube, bounds, filename)


text:
@testset "Init Buoyancy" begin

# Resolution [even integer]
nx, ny  = 64, 64
angle_grid = π/4
omega = 2π/(24*60*60)
f0 = 2 * omega * sin( angle_grid )

grid = Grid( 1e6, nx, 1e6, ny )

# Gather parameters in the structure model
model = SQG( grid, f0 )

# Initial condition for the buoyancy
init_buoyancy!(model)

@test true

end

ext: jl


text:
using Test

using BirkhoffDecomposition

# Generate a random doubly stochastic matrix (n is the dimension)
n  = 3;
x  = randomDoublyStochasticMatrix(n);

# Find decomposition
P, w = birkdecomp(x);


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

const CC = Core.Compiler
import Core: MethodInstance, CodeInstance
import .CC: WorldRange, WorldView

# define new `AbstractInterpreter` that satisfies the minimum interface requirements
# while managing its cache independently
macro newinterp(name)
    cachename = gensym(string(name, "Cache"))
    name = esc(name)
    quote
        struct $cachename
            dict::IdDict{MethodInstance,CodeInstance}
        end
        struct $name <: CC.AbstractInterpreter
            interp::CC.NativeInterpreter
            cache::$cachename
            $name(world = Base.get_world_counter();
                interp = CC.NativeInterpreter(world),
                cache = $cachename(IdDict{MethodInstance,CodeInstance}())
                ) = new(interp, cache)
        end
        CC.InferenceParams(interp::$name) = CC.InferenceParams(interp.interp)
        CC.OptimizationParams(interp::$name) = CC.OptimizationParams(interp.interp)
        CC.get_world_counter(interp::$name) = CC.get_world_counter(interp.interp)
        CC.get_inference_cache(interp::$name) = CC.get_inference_cache(interp.interp)
        CC.code_cache(interp::$name) = WorldView(interp.cache, WorldRange(CC.get_world_counter(interp)))
        CC.get(wvc::WorldView{<:$cachename}, mi::MethodInstance, default) = get(wvc.cache.dict, mi, default)
        CC.getindex(wvc::WorldView{<:$cachename}, mi::MethodInstance) = getindex(wvc.cache.dict, mi)
        CC.haskey(wvc::WorldView{<:$cachename}, mi::MethodInstance) = haskey(wvc.cache.dict, mi)
        CC.setindex!(wvc::WorldView{<:$cachename}, ci::CodeInstance, mi::MethodInstance) = setindex!(wvc.cache.dict, ci, mi)
    end
end

# `OverlayMethodTable`
# --------------------
import Base.Experimental: @MethodTable, @overlay

@newinterp MTOverlayInterp
@MethodTable(OverlayedMT)
CC.method_table(interp::MTOverlayInterp) = CC.OverlayMethodTable(CC.get_world_counter(interp), OverlayedMT)

@overlay OverlayedMT sin(x::Float64) = 1
@test Base.return_types((Int,), MTOverlayInterp()) do x
    sin(x)
end == Any[Int]

max_stars_repo_name: shngt/julia

max_issues_repo_name: shngt/julia

max_forks_repo_name: shngt/julia


text:
################################################################################
# PML update in 2D - see Coefficients in yee.jl if you are confused
################################################################################
###
# pml[z][x] is the effect on x field from pml at z boundary!
# the nodes are grouped differently at the lower boundary

# ΦE[y][x] is backwards from the notation in books - ΦExy
function updatePML(Φ::PMLAux{T}, F::VecArray{T, 2}, C::Coefficients, comp::Int, N::Int) where T
    #allocate arrays
    maxxx =  zeros(T, N, C.s.y.N); maxyy = zeros(T, C.s.x.N, N);
    minxx =  zeros(T, N, C.s.y.N); minyy = zeros(T, C.s.x.N, N);

    x,y,z = (1,2,3)
    maxyx = cover2d(C.bu[comp], y).*Φ.max[y][x]
                    .+ cover2d(C.cu[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, (size(F[z])[y]-N):size(F[z])[y]), y)), y, 2:N+1)
    minyx = cover2d(C.bl[comp], y).*Φ.min[y][x]
                    .+ cover2d(C.cl[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, 1:N+1), y)), y, 1:N)

    x,y,z = (2,3,1)

    maxxy = cover2d(C.bu[comp], z).*Φ.max[z][x]
                    .+ cover2d(C.cu[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, (size(F[y])[z]-N):size(F[y])[z]), z)), z, 2:N+1)
    minxy = cover2d(C.bl[comp], z).*Φ.min[z][x]
                    .+ cover2d(C.cl[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, 1:N+1), z)), z, 1:N)
    x,y,z = (3,1,2)
    maxxz = cover2d(C.bu[comp], y).*Φ.max[y][x]
                    .+ cover2d(C.cu[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, (size(F[z])[y]-N):size(F[z])[y]), y)), y, 2:N+1)
    minxz = cover2d(C.bl[comp], y).*Φ.min[y][x]
                    .+ cover2d(C.cl[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, 1:N+1), y)), y, 1:N)

    maxyz = cover2d(C.bu[comp], z).*Φ.max[z][x]
                    .+ cover2d(C.cu[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, (size(F[y])[z]-N):size(F[y])[z]), z)), z, 2:N+1)
    minyz = cover2d(C.bl[comp], z).*Φ.min[z][x]
                    .+ cover2d(C.cl[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, 1:N+1), z)), z, 1:N)

    return PMLAux([VecField3(maxxx, maxxy, maxxz), VecField3(maxyx, maxyy, maxyz)],
                      [VecField3(minxx, minxy, minxz), VecField3(minyx, minyy, minyz)])
end

function applyPML(F::VecArray{T,2}, Φ::PMLAux{T}, C::Coefficients, comp::Int, N::Int) where T
    Fx = deepcopy(F[1]); Fy = deepcopy(F[2]); Fz = deepcopy(F[3])

    x,y,z = (1,2,3)
    selectdim(Fx, y, 1:N) .= selectdim(F[x], y, 1:N)
                                    .+ C.s.x.Δ .* selectdim(C[comp][x+3], y, 1:N).*Φ.min[y][x]
    selectdim(Fx, y, (size(F[x],y)-N+1):size(F[x],y)) .= selectdim(F[x], y, (size(F[x],y)-N+1):size(F[x],y))
                                                            .+  C.s.x.Δ .* selectdim(C[comp][x+3], y, (size(F[x],y)-N+1):size(F[x],y)).*Φ.max[y][x]

   x,y,z = (2,3,1)

   selectdim(Fy, z, 1:N) .= selectdim(F[x], z, 1:N)
                               .-  C.s.x.Δ .* selectdim(C[comp][x+3], z, 1:N).*Φ.min[z][x]
   selectdim(Fy, z, (size(F[x],z)-N+1):size(F[x],z)) .= selectdim(F[x], z, (size(F[x],z)-N+1):size(F[x],z))
                                                           .-  C.s.x.Δ .*selectdim(C[comp][x+3], z, (size(F[x],z)-N+1):size(F[x],z)).*Φ.max[z][x]
   x,y,z = (3,1,2)
   selectdim(Fz, y, 1:N) .= selectdim(F[x], y, 1:N)
                                   .+ C.s.x.Δ .* selectdim(C[comp][x+3], y, 1:N).*Φ.min[y][x]
   selectdim(Fz, y, (size(F[x],y)-N+1):size(F[x],y)) .= selectdim(F[x], y, (size(F[x],y)-N+1):size(F[x],y))
                                                           .+  C.s.x.Δ .* selectdim(C[comp][x+3], y, (size(F[x],y)-N+1):size(F[x],y)).*Φ.max[y][x]

   selectdim(Fz, z, 1:N) .= selectdim(F[x], z, 1:N)
                               .-  C.s.x.Δ .* selectdim(C[comp][x+3], z, 1:N).*Φ.min[z][x]
   selectdim(Fz, z, (size(F[x],z)-N+1):size(F[x],z)) .= selectdim(F[x], z, (size(F[x],z)-N+1):size(F[x],z))
                                                           .-  C.s.x.Δ .*selectdim(C[comp][x+3], z, (size(F[x],z)-N+1):size(F[x],z)).*Φ.max[z][x]
    return (Fx, Fy, Fz)
end

function updatePML(Φ::PMLAux{T}, F::VecArray{T,3}, C::Coefficients, comp::Int, N::Int) where T
   #allocate arrays
   maxxx =  zeros(T, N, C.s.y.N, C.s.z.N); maxyy = zeros(T, C.s.x.N, N, C.s.z.N); maxzz = zeros(T, C.s.x.N, C.s.y.N, N);
   minxx =  zeros(T, N, C.s.y.N, C.s.z.N); minyy = zeros(T, C.s.x.N, N, C.s.z.N); minzz = zeros(T, C.s.x.N, C.s.y.N, N);

   x,y,z = (1,2,3)
   maxyx = cover3d(C.bu[comp], y).*Φ.max[y][x]
                   .+ cover3d(C.cu[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, (size(F[z])[y]-N):size(F[z])[y]), y)), y, 2:N+1)
   minyx = cover3d(C.bl[comp], y).*Φ.min[y][x]
                   .+ cover3d(C.cl[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, 1:N+1), y)), y, 1:N)

   maxzx = cover3d(C.bu[comp], z).*Φ.max[z][x]
                   .+ cover3d(C.cu[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, (size(F[y])[z]-N):size(F[y])[z]), z)), z, 2:N+1)
   minzx = cover3d(C.bl[comp], z).*Φ.min[z][x]
                   .+ cover3d(C.cl[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, 1:N+1), z)), z, 1:N)

   x,y,z = (2,3,1)
   maxzy = cover3d(C.bu[comp], y).*Φ.max[y][x]
                   .+ cover3d(C.cu[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, (size(F[z])[y]-N):size(F[z])[y]), y)), y, 2:N+1)
   minzy = cover3d(C.bl[comp], y).*Φ.min[y][x]
                   .+ cover3d(C.cl[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, 1:N+1), y)), y, 1:N)

   maxxy = cover3d(C.bu[comp], z).*Φ.max[z][x]
                   .+ cover3d(C.cu[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, (size(F[y])[z]-N):size(F[y])[z]), z)), z, 2:N+1)
   minxy = cover3d(C.bl[comp], z).*Φ.min[z][x]
                   .+ cover3d(C.cl[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, 1:N+1), z)), z, 1:N)
   x,y,z = (3,1,2)
   maxxz = cover3d(C.bu[comp], y).*Φ.max[y][x]
                   .+ cover3d(C.cu[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, (size(F[z])[y]-N):size(F[z])[y]), y)), y, 2:N+1)
   minxz = cover3d(C.bl[comp], y).*Φ.min[y][x]
                   .+ cover3d(C.cl[comp], y) .*selectdim(( ⊗(C.δ[comp], selectdim(F[z], y, 1:N+1), y)), y, 1:N)

   maxyz = cover3d(C.bu[comp], z).*Φ.max[z][x]
                   .+ cover3d(C.cu[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, (size(F[y])[z]-N):size(F[y])[z]), z)), z, 2:N+1)
   minyz = cover3d(C.bl[comp], z).*Φ.min[z][x]
                   .+ cover3d(C.cl[comp], z) .*selectdim(( ⊗(C.δ[comp], selectdim(F[y], z, 1:N+1), z)), z, 1:N)

   return PMLAux{T}([VecField3(maxxx, maxxy, maxxz), VecField3(maxyx, maxyy, maxyz), VecField3(maxzx, maxzy, maxzz)],
                     [VecField3(minxx, minxy, minxz), VecField3(minyx, minyy, minyz), VecField3(minzx, minzy, minzz)])
end
################################################################################
# PML update in 3d - see Coefficients in yee.jl if you are confused
################################################################################
###
# pml[z][x] is the effect on x field from pml at z boundary!
# the nodes are grouped differently at the lower boundary

# ΦE[y][x] is backwards from the notation in books - ΦExy

function applyPML(F::VecArray{T,3}, Φ::PMLAux{T}, C::Coefficients, comp::Int, N::Int) where T
   Fx = deepcopy(F[1]); Fy = deepcopy(F[2]); Fz = deepcopy(F[3])

   x,y,z = (1,2,3)
   selectdim(Fx, y, 1:N) .= selectdim(F[x], y, 1:N)
                                   .+ C.s.x.Δ .* selectdim(C[comp][x+3], y, 1:N).*Φ.min[y][x]
   selectdim(Fx, y, (size(F[x],y)-N+1):size(F[x],y)) .= selectdim(F[x], y, (size(F[x],y)-N+1):size(F[x],y))
                                                           .+  C.s.x.Δ .* selectdim(C[comp][x+3], y, (size(F[x],y)-N+1):size(F[x],y)).*Φ.max[y][x]

   selectdim(Fx, z, 1:N) .= selectdim(F[x], z, 1:N)
                               .-  C.s.x.Δ .* selectdim(C[comp][x+3], z, 1:N).*Φ.min[z][x]
   selectdim(Fx, z, (size(F[x],z)-N+1):size(F[x],z)) .= selectdim(F[x], z, (size(F[x],z)-N+1):size(F[x],z))
                                                           .-  C.s.x.Δ .*selectdim(C[comp][x+3], z, (size(F[x],z)-N+1):size(F[x],z)).*Φ.max[z][x]

  x,y,z = (2,3,1)
  selectdim(Fy, y, 1:N) .= selectdim(F[x], y, 1:N)
                                  .+ C.s.x.Δ .* selectdim(C[comp][x+3], y, 1:N).*Φ.min[y][x]
  selectdim(Fy, y, (size(F[x],y)-N+1):size(F[x],y)) .= selectdim(F[x], y, (size(F[x],y)-N+1):size(F[x],y))
                                                          .+  C.s.x.Δ .* selectdim(C[comp][x+3], y, (size(F[x],y)-N+1):size(F[x],y)).*Φ.max[y][x]

  selectdim(Fy, z, 1:N) .= selectdim(F[x], z, 1:N)
                              .-  C.s.x.Δ .* selectdim(C[comp][x+3], z, 1:N).*Φ.min[z][x]
  selectdim(Fy, z, (size(F[x],z)-N+1):size(F[x],z)) .= selectdim(F[x], z, (size(F[x],z)-N+1):size(F[x],z))
                                                          .-  C.s.x.Δ .*selectdim(C[comp][x+3], z, (size(F[x],z)-N+1):size(F[x],z)).*Φ.max[z][x]
  x,y,z = (3,1,2)
  selectdim(Fz, y, 1:N) .= selectdim(F[x], y, 1:N)
                                  .+ C.s.x.Δ .* selectdim(C[comp][x+3], y, 1:N).*Φ.min[y][x]
  selectdim(Fz, y, (size(F[x],y)-N+1):size(F[x],y)) .= selectdim(F[x], y, (size(F[x],y)-N+1):size(F[x],y))
                                                          .+  C.s.x.Δ .* selectdim(C[comp][x+3], y, (size(F[x],y)-N+1):size(F[x],y)).*Φ.max[y][x]

  selectdim(Fz, z, 1:N) .= selectdim(F[x], z, 1:N)
                              .-  C.s.x.Δ .* selectdim(C[comp][x+3], z, 1:N).*Φ.min[z][x]
  selectdim(Fz, z, (size(F[x],z)-N+1):size(F[x],z)) .= selectdim(F[x], z, (size(F[x],z)-N+1):size(F[x],z))
                                                          .-  C.s.x.Δ .*selectdim(C[comp][x+3], z, (size(F[x],z)-N+1):size(F[x],z)).*Φ.max[z][x]
   return (Fx, Fy, Fz)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

abstract type AbstractConfig end

struct Config <: AbstractConfig
    cursor::Char
    up_arrow::Char
    down_arrow::Char
    scroll_wrap::Bool
    ctrl_c_interrupt::Bool
end

struct MultiSelectConfig <: AbstractConfig
    config::Config
    checked::String
    unchecked::String
end

"""
    Config(; scroll_wrap=false, ctrl_c_interrupt=true, charset=:ascii, cursor::Char, up_arrow::Char, down_arrow::Char)

Configure behavior for selection menus via keyword arguments:

- `scroll_wrap`, if `true`, causes the menu to wrap around when scrolling above the first
  or below the last entry
- `ctrl_c_interrupt`, if `true`, throws an `InterruptException` if the user hits Ctrl-C
  during menu selection. If `false`, [`TerminalMenus.request`](@ref) will return the
  default result from [`TerminalMenus.selected`](@ref).
- `charset` affects the default values for `cursor`, `up_arrow`, and `down_arrow`,
  and can be `:ascii` or `:unicode`
- `cursor` is the character printed to indicate the option that will be chosen by
  hitting "Enter." Defaults are '>' or '→', depending on `charset`.
- `up_arrow` is the character printed when the display does not include the first entry.
  Defaults are '^' or '↑', depending on `charset`.
- `down_arrow` is the character printed when the display does not include the last entry.
  Defaults are 'v' or '↓', depending on `charset`.

Subtypes of `ConfiguredMenu` will print `cursor`, `up_arrow`, and `down_arrow` automatically
as needed, your `writeline` method should not print them.

!!! compat Julia 1.6
    `Config` is available as of Julia 1.6. On older releases use the global `CONFIG`.
"""
function Config(;
                charset::Symbol = :ascii,
                cursor::Char = '\0',
                up_arrow::Char = '\0',
                down_arrow::Char = '\0',
                scroll_wrap::Bool = false,
                ctrl_c_interrupt::Bool = true)
    charset === :ascii || charset === :unicode ||
        throw(ArgumentError("charset should be :ascii or :unicode, received $charset"))
    if cursor == '\0'
        cursor = charset === :ascii ? '>' : '→'
    end
    if up_arrow == '\0'
        up_arrow = charset === :ascii ? '^' : '↑'
    end
    if down_arrow == '\0'
        down_arrow = charset === :ascii ? 'v' : '↓'
    end
    return Config(cursor, up_arrow, down_arrow, scroll_wrap, ctrl_c_interrupt)
end

"""
    MultiSelectConfig(; charset=:ascii, checked::String, unchecked::String, kwargs...)

Configure behavior for a multiple-selection menu via keyword arguments:

- `checked` is the string to print when an option has been selected.
  Defaults are "[X]" or "✓", depending on `charset`.
- `unchecked` is the string to print when an option has not been selected.
  Defaults are "[ ]" or "⬚", depending on `charset`.

All other keyword arguments are as described for [`TerminalMenus.Config`](@ref).
`checked` and `unchecked` are not printed automatically, and should be printed by
your `writeline` method.

!!! compat Julia 1.6
    `MultiSelectConfig` is available as of Julia 1.6. On older releases use the global `CONFIG`.
"""
function MultiSelectConfig(;
                           charset::Symbol = :ascii,
                           checked::String = "",
                           unchecked::String = "",
                           kwargs...)
    charset === :ascii || charset === :unicode || throw(ArgumentError("charset should be :ascii or :unicode, received $charset"))
    if isempty(checked)
        checked = charset === :ascii ? "[X]" : "✓"
    end
    if isempty(unchecked)
        unchecked = charset === :ascii ? "[ ]" : "⬚"
    end
    return MultiSelectConfig(Config(; charset=charset, kwargs...), checked, unchecked)
end

## Below is the old-style CONFIG interface, kept for backwards compatibility.
## Not recommended for any new menu types.

"""
    CONFIG

Global menu configuration parameters

!!! compat Julia 1.6
    `CONFIG` is deprecated, instead configure menus via their constructors.
"""
const CONFIG = Dict{Symbol,Union{Char,String,Bool}}()

"""
    config( <see arguments> )

Keyword-only function to configure global menu parameters

# Arguments
 - `charset::Symbol=:na`: ui characters to use (`:ascii` or `:unicode`); overridden by other arguments
 - `cursor::Char='>'|'→'`: character to use for cursor
 - `up_arrow::Char='^'|'↑'`: character to use for up arrow
 - `down_arrow::Char='v'|'↓'`: character to use for down arrow
 - `checked::String="[X]"|"✓"`: string to use for checked
 - `unchecked::String="[ ]"|"⬚")`: string to use for unchecked
 - `scroll::Symbol=:nowrap`: If `:wrap` wrap cursor around top and bottom, if :`nowrap` do not wrap cursor
 - `supress_output::Bool=false`: Ignored legacy argument, pass `suppress_output` as a keyword argument to `request` instead.
 - `ctrl_c_interrupt::Bool=true`: If `false`, return empty on ^C, if `true` throw InterruptException() on ^C

!!! compat Julia 1.6
    As of Julia 1.6, `config` is deprecated. Use `Config` or `MultiSelectConfig` instead.
"""
function config(;charset::Symbol = :na,
                scroll::Symbol = :na,
                cursor::Char = '\0',
                up_arrow::Char = '\0',
                down_arrow::Char = '\0',
                checked::String = "",
                unchecked::String = "",
                supress_output::Union{Nothing, Bool}=nothing,   # typo was documented, unfortunately
                ctrl_c_interrupt::Union{Nothing, Bool}=nothing)

    if charset === :ascii
        cursor     = '>'
        up_arrow   = '^'
        down_arrow = 'v'
        checked    = "[X]"
        unchecked  = "[ ]"
    elseif charset === :unicode
        cursor     = '→'
        up_arrow   = '↑'
        down_arrow = '↓'
        checked    = "✓"
        unchecked  = "⬚"
    elseif charset === :na
    else
        throw(ArgumentError("charset should be :ascii or :unicode, received $charset"))
    end

    scroll ∉ [:na, :wrap, :nowrap] && throw(ArgumentError("scroll must be :wrap or :nowrap, received $scroll"))
    scroll === :wrap   && (CONFIG[:scroll_wrap] = true)
    scroll === :nowrap && (CONFIG[:scroll_wrap] = false)

    cursor     != '\0' && (CONFIG[:cursor]     = cursor)
    up_arrow   != '\0' && (CONFIG[:up_arrow]   = up_arrow)
    down_arrow != '\0' && (CONFIG[:down_arrow] = down_arrow)
    checked    != ""   && (CONFIG[:checked]    = checked)
    unchecked  != ""   && (CONFIG[:unchecked]  = unchecked)
    supress_output isa Bool   && (CONFIG[:supress_output] = supress_output)
    ctrl_c_interrupt isa Bool && (CONFIG[:ctrl_c_interrupt] = ctrl_c_interrupt)

    return nothing
end

# Set up defaults
config(charset=:ascii, scroll=:nowrap, supress_output=false, ctrl_c_interrupt=true)

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
################################################################################

"Feed on hidden Markov model object"
function feed!(self::HMM, frame::Int64, d::Array{Float64, 2}; args::HMMParams)
  for ix in eachindex(self.tbM)
    plus = euclDist(self.dataM[ix], d[frame, :])
    for jx in eachindex(self.tbM)
      lpen = copy(args.pen)
      if (jx == ix) lpen = 0 end
      if self.tbM[jx][frame + 1] < 0. || self.tbM[ix][frame] + plus + lpen < self.tbM[jx][frame + 1]
        self.tbM[jx][frame + 1] = self.tbM[ix][frame] + plus + lpen
      end
    end
  end
end

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

max_stars_repo_name: grabherr/HiddenMarkovModelReaders

max_issues_repo_name: grabherr/HiddenMarkovModelReaders

max_forks_repo_name: grabherr/HiddenMarkovModelReaders


text:
function ar1_sample_mean(N, beta, alpha, s)
    sm = 0.0
    x = beta / (1 - alpha)
    for i in 1:N
        sm += x
        x = beta + alpha * x + s * randn()
    end
    return sm / N
end

N = 1e7
beta = 1.0
alpha = 0.9
s = 1.0
tic()
result = ar1_sample_mean(N, beta, alpha, s)
println("mean = $result")
toc()


text:
mymean(x)=sum(x)/length(x)   # 定义一个函数，求平均值
a1=[1,2,3,4]                 # Vector{Int64}
@time mymean(a1)             # 0.012165 seconds，其中包含 JIT 编译时间
a2 = 2 * a1                  # Vector{Int64}
@time mymean(a2)             # 0.000005 seconds，由于 a2 和 a1 同为 Vector{Int64}
                             # 此次调用不需再次编译
a3=Float64[1,2,3,4]          # Vector{Float64}，不同的类型
@time mymean(a3)             # 0.030347 seconds，针对新的类型，编译出特定的版本
@time mymean(a3)             # 0.000005 seconds，不需编译第二次

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
using DotaBot
using HTTP
using ArgParse
using Sockets
using JSON
using NEAT

s = ArgParseSettings()
@add_arg_table! s begin
    "--breezyIP"
    help = "breezy server IP adress"
    default = "127.0.0.1"
    "--breezyPort"
    help = "breezy server port number"
    default = "8085"
    "--agentIP"
    help = "agent server IP adress"
    default = "127.0.0.1"
    "--agentPort"
    help = "agent server port number"
    default = "8086"
    "--ind"
    help = "NEAT individual DNA file"
    default = "models/best_neat.dna"
    "--cfg"
    help = "NEAT config file"
    default = "config/neat.yaml"
end

args = parse_args(ARGS, s)
cfg = get_config(args["cfg"])
dna = JSON.parsefile(args["ind"])
individual = NEATInd(cfg, JSON.json(dna))

@async HTTP.serve(parse(IPAddr, args["agentIP"]), parse(Int, args["agentPort"])) do request::HTTP.Request
   @show request
   try
       f_act(features) = get_indiv_action(features, individual)
       response = server_handler(request, args, f_act)
       @show response
       return response
   catch e
       println(e)
       return HTTP.Response(200, JSON.json(Dict("done"=>false)))
   end
end

start_game(args)
while true
    sleep(1)
end


text:
import Base: size, eltype, similar, convert, show, isempty

"""AbstractCvMat{T,N} represents `N`-dimentional arrays in OpenCV (cv::Mat,
cv::UMat, etc), which element type are bound to `T`.
"""
abstract type AbstractCvMat{T,N} <: AbstractArray{T,N} end

# NOTE: subtypes of AbstractCvMat should have `handle` as a member.
handle(m::AbstractCvMat) = m.handle

### Types and methods for C++ types ###

const cvMatExpr = cxxt"cv::MatExpr"
function _size(expr::cvMatExpr)
    s::cvSize = icxx"$expr.size();"
    height(s), width(s)
end
_type(m::cvMatExpr) = icxx"$m.type();"
depth(m::cvMatExpr) = mat_depth(_type(m))
channels(m::cvMatExpr) = mat_channel(_type(m))
eltype(m::cvMatExpr) = jltype(depth(m))

const cvMat = cxxt"cv::Mat"
const cvUMat = cxxt"cv::UMat"
const cvMatVariants = Union{cvMat, cvUMat}

@inline depth(m::cvMatVariants) = convert(Int, icxx"$m.depth();")
@inline channels(m::cvMatVariants) = convert(Int, icxx"$m.channels();")
@inline eltype(m::cvMatVariants) = jltype(depth(m))
@inline dims(m::cvMatVariants) = convert(Int, icxx"$m.dims;")

@inline data(m::cvMat) = icxx"$m.data;"

import Cxx: CppEnum

const UMatUsageFlags = CppEnum{Symbol("cv::UMatUsageFlags"),Int32}

const USAGE_DEFAULT = UMatUsageFlags(0)
const USAGE_ALLOCATE_HOST_MEMORY = UMatUsageFlags(1 << 0)
const USAGE_ALLOCATE_DEVICE_MEMORY = UMatUsageFlags(1 << 1)
const USAGE_ALLOCATE_SHARED_MEMORY = UMatUsageFlags(1 << 2)
const __UMAT_USAGE_FLAGS_32BIT = UMatUsageFlags(0x7fffffff)

cvUMat(usage_flags::UMatUsageFlags=USAGE_DEFAULT) = icxx"cv::UMat($usage_flags);"
cvUMat(m::cvUMat) = icxx"cv::UMat($m);"
usageFlags(m::cvUMat) = icxx"$m.usageFlags;"::UMatUsageFlags
elemSize1(m::cvUMat) = convert(Int, icxx"$m.elemSize1();")

### Methods for AbstractCvMat ###

flags(m::AbstractCvMat) = icxx"$(m.handle).flags;"
dims(m::AbstractCvMat) = icxx"$(m.handle).dims;"
rows(m::AbstractCvMat) = convert(Int, icxx"$(m.handle).rows;")
cols(m::AbstractCvMat) = convert(Int, icxx"$(m.handle).cols;")
_size(m::AbstractCvMat) = (Int(rows(m)), Int(cols(m)))

function size(m::Union{AbstractCvMat, cvMatExpr})
    cn = channels(m)
    if cn == 1
        return _size(m)
    else
        return (_size(m)..., cn)
    end
end

clone(m::AbstractCvMat) = Mat(icxx"$(m.handle).clone();")
total(m::AbstractCvMat) = convert(Int, icxx"$(m.handle).total();")
isContinuous(m::AbstractCvMat) = icxx"$(m.handle).isContinuous();"
elemSize(m::AbstractCvMat) = convert(Int, icxx"$(m.handle).elemSize();")
depth(m::AbstractCvMat) = convert(Int, icxx"$(m.handle).depth();")
channels(m::AbstractCvMat) = channels(handle(m))
empty(m::AbstractCvMat) = icxx"$(m.handle).empty();"
isempty(m::AbstractCvMat) = empty(m)

### MatExpr{T,N} ###

"""MatExpr{T,N} represents cv::MatExpr with encoded type information

`T` and `N` represents the element type and the dimension of Mat,
respectively.

TODO: should consder wherther I make this a subtype of AbstractCvMat{T,N}
"""
mutable struct MatExpr{T,N}
    handle::cvMatExpr
end

handle(m::MatExpr) = m.handle
channels(m::MatExpr) = channels(handle(m))
size(m::MatExpr) = size(handle(m))

function MatExpr(handle::cvMatExpr)
    # determin type parameters by value
    T = eltype(handle)
    N = length(size(handle))
    MatExpr{T,N}(handle)
end

"""Mat{T,N} represents cv::Mat with encoded type information

Mat{T,N} keeps cv::Mat instance with: element type `T` and dimention `N`.
Hence, in fact it behaves like cv::Mat_<T>. Note that Mat stores its
internal data in column-major order, while Julia's arrays are in row-major.

NOTE: Mat{T,N} supports multi-channel 2-dimentional matrices and
single-channel 2-dimentional matrices for now. Should be extended for
N-dimentional cases.
"""
mutable struct Mat{T,N} <: AbstractCvMat{T,N}
    handle::cvMat
end

### Constructors ###

"""Generic constructor"""
function Mat(handle::cvMat)
    # Determine dimention and element type by value and encode eit into type
    T = eltype(handle)
    cn = channels(handle)
    dims = Int(icxx"$handle.dims;")
    N = (cn == 1) ? dims : dims + 1
    Mat{T,N}(handle)
end

"""Empty mat constructor"""
function Mat{T,N}() where {T,N}
    handle = icxx"cv::Mat();"
    Mat{T,N}(handle)
end
function Mat{T}() where T
    handle = icxx"cv::Mat();"
    Mat{T,2}(handle)
end

"""Single-channel 2-dimentional mat constructor"""
function Mat{T}(rows::Int, cols::Int) where T
    typ = maketype(cvdepth(T), 1)
    handle = icxx"cv::Mat($rows, $cols, $typ);"
    Mat{T,2}(handle)
end

function Mat{T}(rows::Int, cols::Int, s::AbstractCvScalar) where T
    typ = maketype(cvdepth(T), 1)
    handle = icxx"cv::Mat($rows, $cols, $typ, $s);"
    Mat{T,2}(handle)
end

"""Multi-chanel 2-dimentional mat constructor"""
function Mat{T}(rows::Int, cols::Int, cn::Int) where T
    typ = maketype(cvdepth(T), cn)
    handle = icxx"cv::Mat($rows, $cols, $typ);"
    Mat{T,3}(handle)
end

function Mat{T}(rows::Int, cols::Int, cn::Int, s::AbstractCvScalar) where T
    typ = maketype(cvdepth(T), cn)
    handle = icxx"cv::Mat($rows, $cols, $typ, $s);"
    Mat{T,3}(handle)
end

"""Single-channel 2-dimentaionl mat constructor with user provided data"""
function Mat{T}(rows::Int, cols::Int, data::Ptr{T}, step=0) where T
    typ = maketype(cvdepth(T), 1)
    handle = icxx"cv::Mat($rows, $cols, $typ, $data, $step);"
    Mat{T,2}(handle)
end

"""Multi-channel 2-dimentaionl mat constructor with user provided data"""
function Mat{T}(rows::Int, cols::Int, cn::Int, data::Ptr{T},
        step=0) where T
    typ = maketype(cvdepth(T), cn)
    handle = icxx"cv::Mat($rows, $cols, $typ, $data, $step);"
    Mat{T,3}(handle)
end

Mat(m::Mat{T,N}) where {T,N} = Mat{T,N}(m.handle)

### Fallback show method

function Base.show(io::IO, ::MIME{Symbol("text/plain")}, m::AbstractCvMat)
    print(io, string(typeof(m)))
    print(io, "\n")
    Base.showarray(io, convert(Mat, m)::Mat, false; header=false)
end

### Mat-specific methods ###

Base.IndexStyle(m::Mat) = Base.IndexLinear()

similar(m::Mat{T}) where {T} = Mat{T}(size(m)...)
similar_empty(m::Mat) = similar(m)

Base.show(io::IO, m::Mat) = Base.showarray(io, m, false)

# Note that cv::UMat doesn't have `data` in members.
@inline data(m::Mat) = data(m.handle)

import Base: getindex, setindex!

getindex(m::Mat, i::Int) =
    convert(eltype(m), icxx"$(m.handle).at<$(eltype(m))>($i-1);")
getindex(m::Mat, i::Int, j::Int) =
    convert(eltype(m), icxx"$(m.handle).at<$(eltype(m))>($i-1, $j-1);")
function getindex(m::Mat{T}, i::Int, j::Int, k::Int) where T
    cn = Val{channels(m)}
    convert(eltype(m), icxx"$(m.handle).at<cv::Vec<$T,$cn>>($i-1, $j-1)[$k-1];")
end

setindex!(m::Mat, v, i::Int) =
    icxx"$(m.handle).at<$(eltype(m))>($i-1) = $v;"
setindex!(m::Mat, v, i::Int, j::Int) =
    icxx"$(m.handle).at<$(eltype(m))>($i-1, $j-1) = $v;"
function setindex!(m::Mat{T}, v, i::Int, j::Int, k::Int) where T
    cn = Val{channels(m)}
    icxx"$(m.handle).at<cv::Vec<$T,$cn>>($i-1, $j-1)[$k-1] = $v;"
end

### UMat{T,N} ###

"""UMat{T,N} represents cv::UMat with encoded type information

`T` and `N` represents the element type and the dimension of Mat,
respectively.
"""
mutable struct UMat{T,N} <: AbstractCvMat{T,N}
    handle::cvUMat
end

"""Generic constructor"""
function UMat(handle::UMat)
    # Determine dimention and element type by value and encode eit into type
    T = eltype(handle)
    cn = channels(handle)
    dims = Int(icxx"$handle.dims;")
    N = (cn == 1) ? dims : dims + 1
    UMat{T,N}(handle)
end

# """Empty mat constructor"""
function UMat{T}(;usage_flags::UMatUsageFlags=USAGE_DEFAULT) where T
    handle = icxx"cv::UMat($usage_flags);"
    UMat{T,0}(handle)
end

# """Single-channel 2-dimentional mat constructor"""
function UMat{T}(rows::Int, cols::Int;
             usage_flags::UMatUsageFlags=USAGE_DEFAULT) where T
    typ = maketype(cvdepth(T), 1)
    handle = icxx"cv::UMat($rows, $cols, $typ, $usage_flags);"
    UMat{T,2}(handle)
end

# """Multi-chanel 2-dimentional mat constructor"""
function UMat{T}(rows::Int, cols::Int, cn::Int;
             usage_flags::UMatUsageFlags=USAGE_DEFAULT) where T
    typ = maketype(cvdepth(T), cn)
    handle = icxx"cv::UMat($rows, $cols, $typ, $usage_flags);"
    UMat{T,3}(handle)
end

# UMat(m::UMat{T,N}) where {T,N} = UMat{T,N}(m.handle)
function UMat(m::Mat{T,N}, flags=ACCESS_READ) where {T,N}
    UMat{T,N}(
        icxx"$(m.handle).getUMat($flags);")
end

similar(m::UMat{T}) where {T} = UMat{T}(size(m)...)
similar_empty(m::UMat) = similar(m)

### MatExpr{T,N} to Mat{T,N} conversion ###

convert(::Type{Mat}, m::MatExpr) = Mat(
    icxx"return cv::Mat($(m.handle));")
function show(io::IO, m::MatExpr)
    print(io, string(typeof(m)))
    print(io, "\n")
    Base.showarray(io, convert(Mat, m)::Mat, false; header=false)
end

### AbstractCvMat{T,N} to MatExpr{T,N}

convert(::Type{MatExpr}, m::AbstractCvMat) = MatExpr(
    icxx"return cv::MatExpr($(m.handle));")
convert(::Type{MatExpr{T,N}}, m::AbstractCvMat{T,N}) where {T,N} = MatExpr{T,N}(
    icxx"return cv::MatExpr($(m.handle));")
convert(::Type{MatExpr}, m::UMat) = convert(MatExpr, convert(Mat, m))
convert(::Type{MatExpr{T,N}}, m::UMat{T,N}) where {T,N} = MatExpr{T,N}(
    icxx"return cv::MatExpr($(m.handle));")

### UMat{T,N} to Mat{T,N} conversion ###

convert(::Type{Mat}, m::UMat{T,N}, flags=ACCESS_READ) where {T,N} = Mat{T,N}(
    icxx"$(m.handle).getMat($flags);")

### Array{T,N} to Mat conversion ###

function convert(::Type{Mat}, arr::Array{T,N}) where {T,N}
    if N != 2 && N != 3
        error("Not supported conversion for now")
    end
    Mat{T}(reverse(size(arr))..., pointer(arr))
end

### Mat to Array{T,N} conversion

@generated function convert(::Type{Array}, m::Mat{T,N}) where {T,N}
    pt = Ptr{T}
    quote
        # TODO: not sure if this is a julia bug
        p = icxx"$(m.handle).data;"
        p = reinterpret($pt, icxx"$(m.handle).data;")
        unsafe_wrap(Array, p, reverse(size(m)))
    end
end

### Matrix operations ###

import Base: +, -, *, /, transpose, ctranspose
import Base: promote_rule, broadcast

promote_rule(::Type{Mat{T,N}}, ::Type{MatExpr{T,N}}) where {T,N} = MatExpr{T,N}
promote_rule(::Type{UMat{T,N}}, ::Type{MatExpr{T,N}}) where {T,N} = MatExpr{T,N}

# For convenience
@inline handle(x::Number) = x

# to avoid method ambiguity warnings
for op in [:+, :-, :*]
    @eval begin
        $op(x::AbstractCvMat{Bool}, y::Bool) = error("not supported")
        $op(x::Bool, y::AbstractCvMat{Bool}) = error("not supported")
    end
end

+(x::MatExpr, y::MatExpr) = MatExpr(icxx"$(x.handle) + $(y.handle);")
+(x::MatExpr, y::Number) = MatExpr(icxx"$(x.handle) + $(y);")
+(x::Number, y::MatExpr) = MatExpr(icxx"$(x) + $(y.handle);")
+(x::MatExpr) = x
+(x::AbstractCvMat) = MatExpr(x)

-(x::MatExpr, y::MatExpr) = MatExpr(icxx"$(x.handle) - $(y.handle);")
-(x::MatExpr, y::Number) = MatExpr(icxx"$(x.handle) - $(y);")
-(x::Number, y::MatExpr) = MatExpr(icxx"$(x) - $(y.handle);")
-(x::MatExpr) = MatExpr(
    icxx"0 - $(x.handle);")
-(x::AbstractCvMat) = -(MatExpr(x))

function broadcast(::typeof(*), x::MatExpr, y::MatExpr)
    MatExpr(icxx"$(x.handle).mul($(y.handle));")
end

broadcast(::typeof(*), x::MatExpr, y::Number) = MatExpr(icxx"$(x.handle) * $(y);")
broadcast(::typeof(*), x::Number, y::MatExpr) = MatExpr(icxx"$(x) * $(y.handle);")

broadcast(::typeof(/), x::MatExpr, y::MatExpr) = MatExpr(icxx"$(x.handle) / $(y.handle);")
broadcast(::typeof(/), x::MatExpr, y::Number) = MatExpr(icxx"$(x.handle) / $(y);")
broadcast(::typeof(/), x::Number, y::MatExpr) = MatExpr(icxx"$(x) / $(y.handle);")

*(x::MatExpr, y::MatExpr) = MatExpr(icxx"$(x.handle) * $(y.handle);")
*(x::MatExpr, y::Number) = broadcast(*, x, y)
*(x::Number, y::MatExpr) = broadcast(*, x, y)
/(x::MatExpr, y::Number) = broadcast(/, x, y)
/(x::Number, y::MatExpr) = broadcast(/, x, y)

# For AbstractCvMats
for op in [:+, :-, :*]
    @eval begin
        @inline $op(x::AbstractCvMat, y::AbstractCvMat) =
            $op(MatExpr(x), MatExpr(y))
        @inline $op(x::MatExpr, y::AbstractCvMat) = $op(promote(x, y)...)
        @inline $op(x::AbstractCvMat, y::MatExpr) = $op(promote(x, y)...)
    end
end
for op in [:/, :*]
    @eval begin
        broadcast(::typeof($op), x::AbstractCvMat, y::AbstractCvMat) =
            $op(MatExpr(x), MatExpr(y))
        broadcast(::typeof($op), x::MatExpr, y::AbstractCvMat) = $op(promote(x, y)...)
        broadcast(::typeof($op), x::AbstractCvMat, y::MatExpr) = $op(promote(x, y)...)
    end
end

# Mat and scalars
for op in [:+, :-, :*]
    @eval begin
        @inline $op(x::AbstractCvMat, y::Number) = $op(MatExpr(x), y)
        @inline $op(x::Number, y::AbstractCvMat) = $op(x, MatExpr(y))
    end
end
for op in [:/, :*]
    @eval begin
        broadcast(::typeof($op), x::AbstractCvMat, y::Number) = $op(MatExpr(x), y)
        broadcast(::typeof($op), x::Number, y::AbstractCvMat) = $op(x, MatExpr(y))
    end
end

transpose(x::MatExpr) = MatExpr(
    icxx"$(x.handle).t();")
transpose(x::AbstractCvMat) = transpose(MatExpr(x))
ctranspose(x::Union{MatExpr, AbstractCvMat}) = transpose(x)

### Linear algebra ###

import Base: inv, ^

inv(x::MatExpr, method=DECOMP_SVD) = MatExpr(
    icxx"$(x.handle).inv($method);")
inv(x::AbstractCvMat, method=DECOMP_SVD) = inv(MatExpr(x), method)

^(x::MatExpr, p::Integer) = (p == -1) ? inv(x) : error("not supported")
^(x::AbstractCvMat, p::Integer) = ^(MatExpr(x), p)


text:
using PolygonClipping
pc = PolygonClipping
using Test

@testset "iterators" begin
    poly = Polygon()
    vert1 = Vertex(1,2)
    push!(poly, vert1)
    for i = 1:10
        push!(poly, Vertex(rand(2)))
    end

    i = 0
    for vert in poly
        i += 1
    end
    i

    @test i == 11
    @test length(poly) == 11

    vert1 = poly.start
    vert2 = vert1.next
    for i = 1:10
        @test vert1 === vert2.prev
        vert1 = vert1.next
        vert2 = vert2.next
    end
end

@testset "vertex constructors" begin
    vert1 = Vertex(0,0)
    vert2 = Vertex(2,1)
    poly = Polygon()
    push!(poly, vert1)
    push!(poly, vert2)
    a = Vertex(poly.start, poly.start.next, 0.25)
    @test a.location == [0.5,0.25]
    @test poly.start === vert1
    @test a.prev === vert1
    @test a.next === vert2
    @test vert1.next === a
    @test vert2.prev === a
end

@testset "point in polygon" begin
    poly1 = Polygon()
    push!(poly1, Vertex(0,0))
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(1,1))
    push!(poly1, Vertex(0,1))
    @test isinside(Vertex(0.5,0.5), poly1) == true
    @test isinside(Vertex(-0.5,0.5), poly1) == false
    @test_throws VertexException isinside(Vertex(0,0), poly1)
    @test_throws VertexException isinside(Vertex(1,0), poly1)
    @test_throws EdgeException isinside(Vertex(0.5,0.0), poly1)
    for i = 1:100
        @test isinside(Vertex(rand(2)), poly1) == true
    end

    # add non-convexity
    push!(poly1, Vertex(0.9,0.5))
    @test isinside(Vertex(0.9,0.9), poly1) == true
    @test isinside(Vertex(0.9,0.9), poly1) == true
    @test isinside(Vertex(0.99,0.1), poly1) == true
    @test isinside(Vertex(0.99,0.5), poly1) == true
    @test isinside(Vertex(0.9,0.1), poly1) == true

    # add self-intersection
    push!(poly1, Vertex(1.1, 0.25))
    @test isinside(Vertex(1.09,0.25), poly1) == true
    @test isinside(Vertex(1.09,0.26), poly1) == true
    @test isinside(Vertex(0.9,0.25), poly1) == false
    @test isinside(Vertex(0.9,0.05), poly1) == true
    @test isinside(Vertex(0.9,0.9), poly1) == true
end

@testset "line intersection method" begin
    @test (intersection(Vertex(0,0), Vertex(1,1), Vertex(1,0), Vertex(0,1))
            == (true, 0.5, 0.5))
    @test (intersection(Vertex(0,0), Vertex(1,1), Vertex(1/2,0), Vertex(0,1/2))
            == (true, 0.25, 0.5))
    @test (intersection(Vertex(0,0), Vertex(2,2), Vertex(1/2,0), Vertex(0,1/2))
            == (true, 0.125, 0.5))
    @test (intersection(Vertex(0,0), Vertex(0,1), Vertex(1,0), Vertex(1,1))
            == (false, 0, 0))
    @test (intersection(Vertex(0.5,0), Vertex(0.5,1), Vertex(0,0.5), Vertex(1,0.5))
            == (true, 0.5, 0.5))
    @test (intersection(Vertex(0,0), Vertex(1,1), Vertex(0,0), Vertex(-1,-1))
            == (false, 0, 0))
    @test (intersection(Vertex(0,0), Vertex(0,0), Vertex(1,1), Vertex(1,1))
            == (false, 0, 0))
    @test (intersection(Vertex(0,0), Vertex(1,0), Vertex(0.25,0.5), Vertex(0.25,-0.5))
            == (true, 0.25, 0.5))
    @test (intersection(Vertex(1,0), Vertex(0,0), Vertex(0.25,0.5), Vertex(0.25,-0.5))
            == (true, 0.75, 0.5))
    @test (intersection(Vertex(0,0), Vertex(1,0), Vertex(0.25,-0.5), Vertex(0.25,0.25))
            == (true, 0.25, 2/3))
    @test (intersection(Vertex(0,0), Vertex(1,0), Vertex(0.25,0.25), Vertex(0.25,-0.5))
            == (true, 0.25, 1/3))
    @test (intersection(Vertex(0,0), Vertex(1,0), Vertex(1.5,0.5), Vertex(1.5,-0.5))
            == (false, 0, 0))
    @test_throws DegeneracyException intersection(Vertex(0,0), Vertex(1,0), Vertex(0.5,0), Vertex(0.5,1.0))
end

@testset "Vertex insertion" begin
    vert1 = Vertex(0,0)
    vert2 = Vertex(1,1)
    vert3 = Vertex(vert1, vert2, 0.25)
    @test vert3.next === vert2
    @test vert3.prev === vert1
    @test vert3.location == [0.25,0.25]
    vert3 = Vertex(vert2, vert1, 0.25)
    @test vert3.next === vert1
    @test vert3.prev === vert2
    @test vert3.location == [0.75,0.75]
end

@testset "unprocessed" begin
    poly1 = Polygon()
    push!(poly1, Vertex(0,0))
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(1,1))
    push!(poly1, Vertex(0,1))
    @test unprocessed(poly1) == false
    poly1.start.next.intersect = true
    @test unprocessed(poly1) == true
    poly1.start.next.visited = true
    @test unprocessed(poly1) == false
end

@testset "remove vertices" begin
    poly1 = Polygon()
    vert1 = Vertex(0,1)
    push!(poly1, vert1)
    @test length(poly1) == 1
    @test poly1.start === vert1
    remove(vert1, poly1)
    @test length(poly1) == 0
    @test poly1.start === nothing
    @test vert1.next === nothing
    @test vert1.prev === nothing
    push!(poly1, vert1)
    vert2 = Vertex(1,1)
    push!(poly1, vert2)
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(0,0))
    @test length(poly1) == 4
    remove(vert1, poly1)
    @test length(poly1) == 3
    @test poly1.start === vert2
end

@testset "clipping" begin
    poly1 = Polygon()
    push!(poly1, Vertex(0,1))
    push!(poly1, Vertex(1,1))
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(0,0))
    poly3 = Polygon()
    push!(poly3, Vertex(0.9,-0.1))
    push!(poly3, Vertex(0.9,0.1))
    push!(poly3, Vertex(0.5,-0.05))
    push!(poly3, Vertex(0.1,0.1))
    push!(poly3, Vertex(0.1,-0.1))

    #test phase1
    @test length(poly1) == 4
    @test length(poly3) == 5
    pc.phase1!(poly3, poly1)
    @test length(poly1) == 8
    @test length(poly3) == 9
    intersects = 0
    for p1 in poly1, p2 in poly3
        if p1.location == p2.location
            @test p1.neighbor === p2
            @test p2.neighbor === p1
            @test p1.intersect
            @test p2.intersect
            intersects += 1
        end
    end
    @test intersects == 4 # make sure we found 4 intersections

    pc.phase2!(poly3, poly1)
    pc.phase2!(poly1, poly3)
    entries = 0
    for vert in poly1
        if vert.entry
            entries += 1
        end
    end
    @test entries == 6
    entries = 0
    for vert in poly3
        if vert.entry
            entries += 1
        end
    end
    @test entries == 5

    # Test phase1 with holes
    poly1 = Polygon()
    push!(poly1, Vertex(0.3,1.5))
    push!(poly1, Vertex(0.6,1.5))
    push!(poly1, Vertex(0.6,-1.5))
    push!(poly1, Vertex(0.3,-1.5))
    poly2 = Polygon()
    push!(poly2, Vertex(0.25,0.25))
    push!(poly2, Vertex(0.25,0.75))
    push!(poly2, Vertex(0.75,0.75))
    push!(poly2, Vertex(0.75,0.25))
    poly3 = Polygon()
    push!(poly3, Vertex(0,1))
    push!(poly3, Vertex(1,1))
    push!(poly3, Vertex(1,0))
    push!(poly3, Vertex(0,0))
    pc.phase1!(poly1,poly2,poly3)
    @test length(poly1) == 12
    @test length(poly2) == 8
    @test length(poly3) == 8
    intersects = 0
    for p1 in poly1, p2 in poly2
        if p1.location == p2.location
            @test p1.neighbor === p2
            @test p2.neighbor === p1
            @test p1.intersect
            @test p2.intersect
            intersects += 1
        end
    end
    for p1 in poly1, p2 in poly3
        if p1.location == p2.location
            @test p1.neighbor === p2
            @test p2.neighbor === p1
            @test p1.intersect
            @test p2.intersect
            intersects += 1
        end
    end
    @test intersects == 8 # make sure we found 4 intersections

    poly1 = Polygon()
    push!(poly1, Vertex(0,1))
    push!(poly1, Vertex(1,1))
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(0,0))
    poly3 = Polygon()
    push!(poly3, Vertex(0.9,-0.1))
    push!(poly3, Vertex(0.9,0.1))
    push!(poly3, Vertex(0.5,-0.05))
    push!(poly3, Vertex(0.1,0.1))
    push!(poly3, Vertex(0.1,-0.1))

    results = intersection(poly3, poly1)

    @test length(results) == 2
    @test length(results[1]) == 4
    @test length(results[2]) == 4

    poly1 = Polygon()
    push!(poly1, Vertex(0,1))
    push!(poly1, Vertex(1,1))
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(0,0))

    poly2 = Polygon()
    push!(poly2, Vertex(0.5,0.5))
    push!(poly2, Vertex(1.5,0.5))
    push!(poly2, Vertex(1.5,-0.5))
    push!(poly2, Vertex(0.5,-0.5))
    results = intersection(poly1, poly2)
    @test length(results) == 1
    @test length(results[1]) == 5

    poly1 = Polygon()
    push!(poly1, Vertex(1,0))
    push!(poly1, Vertex(1,1))
    push!(poly1, Vertex(0,1))
    push!(poly1, Vertex(0,0))

    poly2 = Polygon()
    push!(poly2, Vertex(0.5,0.5))
    push!(poly2, Vertex(1.5,0.5))
    push!(poly2, Vertex(1.5,-0.5))
    push!(poly2, Vertex(0.5,-0.5))
    results = intersection(poly1, poly2)
    @test length(results) == 1
    @test length(results[1]) == 5

    # Issue 2

    polya = Polygon()
    push!(polya, Vertex(-0.5, 0))
    push!(polya, Vertex(1.5, 0.0))
    push!(polya, Vertex(1.5, 1))
    push!(polya, Vertex(-0.5, 1))

    polyb = Polygon()
    push!(polyb, Vertex(0, -0.5))
    push!(polyb, Vertex(0, 1.5))
    push!(polyb, Vertex(1, 1.5))
    push!(polyb, Vertex(1, -0.5))

    PolygonClipping.phase1!(polya, polyb)

    intersects = 0
    for p1 in polya, p2 in polyb
        if p1.location == p2.location
            @test p1.neighbor === p2
            @test p2.neighbor === p1
            @test p1.intersect
            @test p2.intersect
            intersects += 1
        end
    end
    @test intersects == 4 # make sure we found 4 intersections

    polya = Polygon()
    push!(polya, Vertex(-0.5, 0))
    push!(polya, Vertex(1.5, 0.0))
    push!(polya, Vertex(1.5, 1))
    push!(polya, Vertex(-0.5, 1))

    polyb = Polygon()
    push!(polyb, Vertex(0, -0.5))
    push!(polyb, Vertex(0, 1.5))
    push!(polyb, Vertex(1, 1.5))
    push!(polyb, Vertex(1, -0.5))

    PolygonClipping.phase1!(polyb, polya)

    intersects = 0
    for p1 in polya, p2 in polyb
        if p1.location == p2.location
            @test p1.neighbor === p2
            @test p2.neighbor === p1
            @test p1.intersect
            @test p2.intersect
            intersects += 1
        end
    end
    @test intersects == 4 # make sure we found 4 intersections

    PolygonClipping.phase2!(polyb, polya)
    PolygonClipping.phase2!(polya, polyb)
end

@testset "infill" begin
    polya = Polygon()
    push!(polya, Vertex(0.0, 1.0))
    push!(polya, Vertex(1.0, 1.0))
    push!(polya, Vertex(1.0, 0.0))
    push!(polya, Vertex(0.0, 0.0))

    polyb = Polygon()
    push!(polyb, Vertex(-0.15, -0.5))
    push!(polyb, Vertex(-0.15, 1.5))
    push!(polyb, Vertex(0.15, 1.5))
    push!(polyb, Vertex(0.15, -0.2))
    push!(polyb, Vertex(0.3, -0.2))
    push!(polyb, Vertex(0.3, 1.5))
    push!(polyb, Vertex(0.45, 1.5))
    push!(polyb, Vertex(0.45, -0.5))
    fill = infill(polya, polyb)
end


text:
using TerminalUI; using VT100

const thisdir = dirname(@__FILE__)

include(Pkg.dir("VT100","test","TerminalRegressionTests.jl"))
using TerminalRegressionTests: load_outputs, compare

function do_compare(em, i, s, output, decorator)
    TerminalUI.redraw(s, i.w)
    isa(output, String) && (output = Vector{UInt8}(output))
    isa(decorator, String) && (decorator = Vector{UInt8}(decorator))

    outbuf = IOBuffer()
    TerminalUI.render(outbuf, s)
    seekstart(outbuf); VT100.parseall!(em, outbuf)
    compare(em, output, decorator)
end

function test_against(em,file)
    output = open(read,joinpath(thisdir,file))
    compare(em, output)
end

function interact(w, inputs, file)
    outputs, decorators = load_outputs(joinpath(thisdir,file))

    outbuf = IOBuffer()
    i = TerminalUI.InlineDialog(w,
        TerminalUI.IOBufferTerminal(outbuf))
    s = TerminalUI.setup_screen(i)
    s.offset = (1,1)
    s.fullsize = Base.displaysize(i.tty)
    em = VT100.ScreenEmulator()
    TerminalUI.focus(i.w)

    (length(decorators) < length(outputs)) &&
        (decorators = Any[decorators; Vector{Void}(length(outputs))])
    @test do_compare(em, i, s, outputs[1], decorators[1])
    @assert length(outputs) >= length(inputs)+1
    # Run over inputs
    for (input, output, decorator) in zip(inputs, outputs[2:end], decorators[2:end])
        inbuf = IOBuffer(input)
        seekstart(inbuf)
        result = :ok
        while result != :done && !eof(inbuf)
            result = Base.LineEdit.match_input(i.w.ctx.focuss.keymap,nothing,inbuf)(s,nothing)
        end
        @test do_compare(em, i, s, output, decorator)
    end
end

# Test cellblock
let
    cells = Vector{Cell}[
      [Cell(TerminalUI.dummycell(:blue,:red); content = 'x') for i = 1:5]
    for i = 1:5];
    w = TerminalUI.CellDisplay(cells);
    em = TerminalUI.debug(w, false)
    @test test_against(em,"outputs/cellblock.out")
end

# Test simple form
let
    myform = form(Pair[
         :name => ask("Enter your name:", AbstractString),
         :address => ask("Enter your address:", AbstractString),
         :age => ask("Enter your age:", Int)
       ],style = :incremental);
    interact(myform.chain, ["The"," Wizard\n","Oz\n","10"], "outputs/form.multiout")
end

# Test simple list
let
    fruit = ["Apple","Pear","Kiwi"]
    w = ListWidget(fruit)
    interact(w, ["\e[B","\e[B"], "outputs/list.multiout")
end


text:
module UHDBindings

using Libdl
using Printf

# ----------------------------------------------------
# --- Architecture management
# ----------------------------------------------------
# We need to have a full word definition for time stamp format which is on 64 bits for 64 bits sytem and 32 bits for 32 bits system.
# For an unknown (at this stage) reason Clonglong is mapped to 64 bits even on ARM device so state a
# a special case if we use software on ARM device
# FIXME => Look for direct appropriate format on the e310
const ARCHI = Sys.CPU_NAME == "cortex-a9" ? "arm" : "pc";
if ARCHI == "arm"
    const FORMAT_LONG = Int32;
else
    const FORMAT_LONG = Clonglong;
end

# ----------------------------------------------------
# --- Artifact for LibUHD
# ----------------------------------------------------
# init globals and lib path
using Pkg.Artifacts;
const libUHD_rootpath = artifact"libUHD";
const libUHD = joinpath(libUHD_rootpath, "libuhd.so");

# ----------------------------------------------------
# --- Common configuration and structures
# ----------------------------------------------------
# --- Including the file
include("common.jl");
export Timestamp
export UHDBinding
# Exporting printing macros
export @infotx, @warntx;
export @inforx, @warnrx;

# ----------------------------------------------------
# --- Finding devices
# ----------------------------------------------------
include("Find.jl")
export uhd_find_devices

# ----------------------------------------------------
# --- Receiver Configuration
# ----------------------------------------------------
# All structures and functions for the Rx side
include("Rx.jl");
# Structures
export UHDRx
# Export functions
export recv,recv!;
export populateBuffer!
export getError, getTimestamp
export restartStreamer

# ----------------------------------------------------
# --- Transmitter Configuration
# ----------------------------------------------------
# All structures and functions for the Tx side
include("Tx.jl");
# Structures
export UHDTx
# Export functions
export send;
"""
Init the core parameter of the radio in Tx or in Rx mode and initiate RF parameters
# --- Syntax
openUHD(mode,sysImage,carrierFreq,samplingRate,txGain,antenna="RX2")
# --- Input parameters
- mode 			: String to open radio in "Tx" (transmitter) or in "Rx" (receive) mode
- carrierFreq	: Desired Carrier frequency [Union{Int,Float64}]
- samplingRate	: Desired bandwidth [Union{Int,Float64}]
- gain		: Desired Gain [Union{Int,Float64}]
- antenna		: Desired Antenna alias [String]
Keywords=
- args	  : String with the additionnal load parameters (for instance, path to the FPHGA image) [String]
# --- Output parameters
- uhd		  	: UHD object [UHDBinding]
"""
function openUHD(carrierFreq, samplingRate, gain, antenna = "RX2";args="")
	# ----------------------------------------------------
	# --- Handler
	# ----------------------------------------------------
	addressUSRP = Ref{Ptr{uhd_usrp}}();
	# --- Cal the init
	@assert_uhd ccall((:uhd_usrp_make, libUHD), uhd_error, (Ptr{Ptr{uhd_usrp}}, Cstring),addressUSRP,args);
	# --- Get the usable object
	pointerUSRP = addressUSRP[];
	# ----------------------------------------------------
	# --- Set Rx stage
	# ----------------------------------------------------
	rx = openUHDRx(pointerUSRP,carrierFreq, samplingRate, gain, antenna;args=args);
	# ----------------------------------------------------
	# --- Set Tx stage
	# ----------------------------------------------------
	tx = openUHDTx(pointerUSRP,carrierFreq, samplingRate, gain, antenna;args=args);
	# ----------------------------------------------------
	# --- Create radio
	# ----------------------------------------------------
	uhdBinding = UHDBinding(addressUSRP,rx,tx);
	return uhdBinding;
end
export openUHD;

"""
Close the USRP device (Rx or Tx mode) and release all associated objects

# --- Syntax

close(uhd)
# --- Input parameters
- uhd	: UHD object [UHDBinding]
# --- Output parameters
- []
"""
function Base.close(uhdBinding::UHDBinding)
	# --- Close Rx and Tx streams
	close(uhdBinding.rx);
	close(uhdBinding.tx);
	# --- Close the USRP main object
	@assert_uhd  ccall((:uhd_usrp_free, libUHD), uhd_error, (Ptr{Ptr{uhd_usrp}},),uhdBinding.addressUSRP);
	# --- Print a flag
	print("\n");
	@info "USRP device is now closed.";
end

"""
Print the radio configuration

# --- Syntax

print(radio)
# --- Input parameters
- radio		: UHD object [Union{UHDBinding,UHDTx,UHDRx}]
# --- Output parameters
- []
"""
function Base.print(radio::UHDBinding)
	# --- Print the configuration of Tx and Rx
	print(radio.rx);
	print(radio.tx);
end

# All functions are defined @UHDTx or @UHDRx level, we should define configuration functions @UHDBinding level
# In such case, we will both apply configuration at Rx and Tx sides.

"""
Update sampling rate of current radio device, and update radio object with the new obtained sampling frequency. If the input parameter is the UHDBinding object, the desired sampling frequency will be applied on both Rx and Tx sides.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx sampling frequency

# --- Syntax

updateSamplingRate!(radio,samplingRate)
# --- Input parameters
- radio	  : UHD device [Union{UHDBinding,UHDRx,UHDTx}]
- samplingRate	: New desired sampling rate
# --- Output parameters
-
"""
function updateSamplingRate!(radio::UHDBinding,samplingRate)
	@sync updateSamplingRate!(radio.rx,samplingRate);
	@sync updateSamplingRate!(radio.tx,samplingRate);
end

"""
Update carrier frequency of current radio device, and update radio object with the new obtained carrier frequency. If the input parameter is the UHDBinding object, the desired carrier frequency will be applied on both Rx and Tx sides.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx carrier frequency

# --- Syntax

updateCarrierFreq!(radio,carrierFreq)
# --- Input parameters
- radio	  : UHD device [Union{UHDBinding,UHDRx,UHDTx}]
- carrierFreq	: New desired carrier frequency
# --- Output parameters
-
"""
function updateCarrierFreq!(radio::UHDBinding,carrierFreq)
	@sync updateCarrierFreq!(radio.rx,carrierFreq);
	@sync updateCarrierFreq!(radio.tx,carrierFreq);
end

"""
Update gain of current radio device, and update radio object with the new obtained gain. If the input parameter is the UHDBinding object, the desired gain will be applied on both Rx and Tx sides.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx gain

# --- Syntax

updateGain!(radio,gain)
# --- Input parameters
- radio	  : UHD device [Union{UHDBinding,UHDRx,UHDTx}]
- gain	: New desired gain
# --- Output parameters
-
"""
function updateGain!(radio::UHDBinding,gain)
	@sync updateGain!(radio.rx,gain);
	@sync updateGain!(radio.tx,gain);
end

# When given to UHDBinding, recv and send will dispatch to the appropriate substructure
# Recv
recv(radio::UHDBinding,nbSamples)  = recv(radio.rx,nbSamples);
recv!(sig,radio::UHDBinding;kwargs...) = recv!(sig,radio.rx;kwargs...);
# Send
send(radio::UHDBinding,params...) = send(radio.tx,params...);

# ----------------------------------------------------
# --- Common functions and structures
# ----------------------------------------------------
export updateSamplingRate!
export updateGain!
export updateCarrierFreq!
export print;
export close;

end # module


text:
# Test MIP
#
#   maximize z
#
#   s.t.  3 p1 + 5 p2 + 4 p3 - z = 0
#         0.5 p1 + 2 p2 + 1 p3 <= 6
#
#         p1 is integer: 0 <= p1
#         p1 is integer: 0 <= p2
#         p1 is integer: 0 <= p3
#         z is binary

using Gurobi, Compat.Test

@testset "test_read" begin
# The function read() is used to produce a MIP start vector
# from a *.sol or *.mst file. These files must first be
# created using write_model().

    ## BUILD MODEL
    #-------------

    simple_model_env = Gurobi.Env()
    setparam!(simple_model_env, "OutputFlag", 0)

    simple_model = Gurobi.Model(simple_model_env, "simple_mip", :maximize)

    add_ivar!(simple_model, 0., 0, Inf)  # p1
    add_ivar!(simple_model, 0., 0, Inf)  # p2
    add_ivar!(simple_model, 0., 0, Inf)  # p3
    add_cvar!(simple_model, 1., 0., Inf) # z
    update_model!(simple_model)

    add_constr!(simple_model, [3., 5., 4., -1.], '<', 0.)
    add_constr!(simple_model, [3., 5., 4., -1.], '>', 0.)

    add_constr!(simple_model, [0.5, 2., 1., 0.], '<', 6.)

    setparam!(simple_model, "Heuristics", 0.0)
    setparam!(simple_model, "Presolve", 0)

    update_model!(simple_model)

    optimize(simple_model)

    first_solution = Gurobi.get_dblattrarray(simple_model, "X", 1, 4)

    ## WRITE SOLUTION FILES
    #----------------------

    write_model(simple_model, "simple_out.sol")
    write_model(simple_model, "simple_out.mst")

    @test isfile("simple_out.sol")
    @test isfile("simple_out.mst")

    ## READ SOLUTION FILE
    #--------------------

    start_1 = Gurobi.get_dblattrarray(simple_model, "Start", 1, 4)
    @test start_1 == [1.0e101, 1.0e101, 1.0e101, 1.0e101]

    Gurobi.read(simple_model, "simple_out.sol")
    update_model!(simple_model)

    start_2 = Gurobi.get_dblattrarray(simple_model, "Start", 1, 4)
    @test start_2 == [12.0, 0.0, 0.0, 36.0]
    @test start_2 == first_solution

    ## VERIFY MIP START
    #------------------

    # The parameters Presolve and Heuristics must be set to 0
    # so the MIP start is used.
    optimize(simple_model)

    ## DELETE FILES
    #--------------

    rm("simple_out.sol")
    rm("simple_out.mst")

end


text:
module TestGlobalScheduler

import Schedulers
using Schedulers.Internal: WorkStealingScheduler, PrioritizedScheduler
using Test

function test_workstealing_trivial()
    scheduler = Schedulers.global_workstealing() do
        Schedulers.current_scheduler()
    end
    @test scheduler isa WorkStealingScheduler
end

function test_prioritized_trivial()
    scheduler = Schedulers.global_prioritized() do
        Schedulers.current_scheduler()
    end
    @test scheduler isa PrioritizedScheduler
end

end  # module


text:
# import Pkg;
# Pkg.add("Dash")
# Pkg.add("DashHtmlComponents")
# Pkg.add("DashCoreComponents")
# Pkg.add("PlotlyJS")
# Pkg.add("JSON")

using Dash, DashHtmlComponents, DashCoreComponents
using PlotlyJS
using JSON
using Dates

#=
function plotPM()
    # NOTE: we must use the semicolon in beginning of plotly commands as everything works via kwargs

    pm1 = scatter(;
                  x = [1, 2, 3, 4, 5],
                  y = [5, 10, 15, 10, 8],
                  mode = "lines+markers",
                  name = "PM 1",
                  marker_color = :red,
                  marker_size = 5,
                  )

    layout = Layout(;
                    title = "PM Levels",
                    xaxis_title = "time",
                    yaxis_title = "PM [μg/m³]",
                    plot_bgcolor = :transparent,
                    paper_bgcolor = :transparent,
                    )
    plot([pm1], layout)
end
=#

#=

function plotContour()

    data = contour(;
                   x = sensors[:OPCN3][:dt],
                   y = binCenters,
                   z = log10.(hcat(sensors[:OPCN3][:bins]...) .+ 1.0 ),  # add 1 to each count so we don't have issues with log(0)
                   colorscale = "Jet",
                   ncontours = 100,
                   contours_showlines = false,
                   )

    layout = Layout(;
                    title = "Particle Size Distribution",
                    xaxis_title = "time",
                    yaxis_title = "Particle Radius [μm]",
                    yaxis_type = "log",
                    plot_bgcolor = :transparent,
                    paper_bgcolor = :transparent,
                    )
    plot(data, layout)
end

=#

#=

function update_pm()
    # load in some test data to play with
    APDS9002()
    BME280()
    GUV001()
    MGS001()
    OPCN3()
    SI114X()

    plotPM()
end

function update_contour()
    plotContour()
end

=#

pm1 = scatter(;
x = [1, 2, 3, 4, 5],
y = [5, 10, 15, 10, 8],
mode = "lines+markers",
name = "PM 1",
marker_color = :red,
marker_size = 5,
);

layout = Layout(;
title = "Particle Size Distribution",
xaxis_title = "time",
yaxis_title = "Particle Radius [μm]",
yaxis_type = "log",
plot_bgcolor = :transparent,
paper_bgcolor = :transparent,
)

app = dash()

app.layout = html_div() do
    html_h1(
        "MINTS Dashboards: UTD Nodes",
        style=Dict("color"=>:black, "textAlign"=>"center"),
    ),
    html_div(
        children = [
            dcc_graph(
                id="pm-graph",
                #figure = plotPM()

                figure = plot([pm1], layout)
            ),
            #dcc_graph(
             #   id="size-contour",
              #  figure = plotContour()
            #),

            dcc_interval(
                id="interval-component",
                interval = 250, # 250 miliseconds
                n_intervals=0
            ),
        ],
#        className = "two-thirds column wind__speed__container",
    )
end

#=
callback!(app,
          Output("pm-graph", "figure"),
          Input("interval-component", "n_intervals")) do n
    update_pm()
end

callback!(app, Output("size-contour", "figure"), Input("interval-component", "n_intervals")) do n
    update_contour()
end
=#

run_server(app, "0.0.0.0", debug=true)


text:
"""
    IdentityMap

An identity map,

```math
    x ↦ x.
```

### Fields

- `dim` -- dimension
"""
struct IdentityMap <: AbstractMap
    dim::Int
end
statedim(m::IdentityMap) = m.dim
outputdim(m::IdentityMap) = m.dim
inputdim(::IdentityMap) = 0
islinear(::IdentityMap) = true
isaffine(::IdentityMap) = true
apply(m::IdentityMap, x) = x

"""
    ConstrainedIdentityMap

An identity map with state constraints of the form:

```math
    x ↦ x, x(t) ∈ \\mathcal{X}.
```

### Fields

- `dim` -- dimension
- `X`   -- state constraints
"""
struct ConstrainedIdentityMap{ST} <: AbstractMap
    dim::Int
    X::ST
end
statedim(m::ConstrainedIdentityMap) = m.dim
stateset(m::ConstrainedIdentityMap) = m.X
outputdim(m::ConstrainedIdentityMap) = m.dim
inputdim(::ConstrainedIdentityMap) = 0
islinear(::ConstrainedIdentityMap) = true
isaffine(::ConstrainedIdentityMap) = true
apply(::ConstrainedIdentityMap, x) = x

"""
    LinearMap

A linear map,

```math
    x ↦ Ax
```

### Fields

- `A` -- matrix
"""
struct LinearMap{T, MT<:AbstractMatrix{T}} <: AbstractMap
    A::MT
end
statedim(m::LinearMap) = size(m.A, 2)
outputdim(m::LinearMap) = size(m.A, 1)
inputdim(::LinearMap) = 0
islinear(::LinearMap) = true
isaffine(::LinearMap) = true
apply(m::LinearMap, x) = m.A * x

"""
    ConstrainedLinearMap

A linear map with state constraints of the form:

```math
    x ↦ Ax, x(t) ∈ \\mathcal{X}.
```

### Fields

- `A` -- matrix
- `X` -- state constraints
"""
struct ConstrainedLinearMap{T, MT<:AbstractMatrix{T}, ST} <: AbstractMap
    A::MT
    X::ST
end
statedim(m::ConstrainedLinearMap) = size(m.A, 2)
stateset(m::ConstrainedLinearMap) = m.X
outputdim(m::ConstrainedLinearMap) = size(m.A, 1)
inputdim(::ConstrainedLinearMap) = 0
islinear(::ConstrainedLinearMap) = true
isaffine(::ConstrainedLinearMap) = true
apply(m::ConstrainedLinearMap, x) = m.A * x

"""
    AffineMap

An affine map,

```math
    x ↦ Ax + b.
```

### Fields

- `A` -- matrix
- `b` -- vector
"""
struct AffineMap{T, MT<:AbstractMatrix{T}, VT<:AbstractVector{T}} <: AbstractMap
    A::MT
    b::VT
    function AffineMap(A::MT, b::VT) where {T, MT<:AbstractMatrix{T}, VT<:AbstractVector{T}}
        @assert size(A, 1) == length(b)
        return new{T, MT, VT}(A, b)
    end
end
statedim(m::AffineMap) = size(m.A, 2)
outputdim(m::AffineMap) = length(m.b)
inputdim(::AffineMap) = 0
islinear(::AffineMap) = false
isaffine(::AffineMap) = true
apply(m::AffineMap, x) = m.A * x + m.b

"""
    ConstrainedAffineMap

An affine map with state constraints of the form:

```math
    x ↦ Ax + b, x(t) ∈ \\mathcal{X}.
```

### Fields

- `A` -- matrix
- `b` -- vector
- `X` -- state constraints
"""
struct ConstrainedAffineMap{T, MT<:AbstractMatrix{T}, VT<:AbstractVector{T}, ST} <: AbstractMap
    A::MT
    b::VT
    X::ST
    function ConstrainedAffineMap(A::MT, b::VT, X::ST) where {T, MT<:AbstractMatrix{T}, VT<:AbstractVector{T}, ST}
        @assert size(A, 1) == length(b)
        return new{T, MT, VT, ST}(A, b, X)
    end
end
statedim(m::ConstrainedAffineMap) = size(m.A, 2)
stateset(m::ConstrainedAffineMap) = m.X
outputdim(m::ConstrainedAffineMap) = length(m.b)
inputdim(::ConstrainedAffineMap) = 0
islinear(::ConstrainedAffineMap) = false
isaffine(::ConstrainedAffineMap) = true
apply(m::ConstrainedAffineMap, x) = m.A * x + m.b

"""
    LinearControlMap

A linear control map,

```math
    (x, u) ↦ Ax + Bu.
```

### Fields

- `A` -- matrix
- `B` -- matrix
"""
struct LinearControlMap{T, MTA<:AbstractMatrix{T}, MTB<:AbstractMatrix{T}} <: AbstractMap
    A::MTA
    B::MTB
    function LinearControlMap(A::MTA, B::MTB) where {T, MTA<:AbstractMatrix{T}, MTB<:AbstractMatrix{T}}
        @assert size(A, 1) == size(B, 1)
        return new{T, MTA, MTB}(A, B)
    end
end
statedim(m::LinearControlMap) = size(m.A, 2)
inputdim(m::LinearControlMap) = size(m.B, 2)
outputdim(m::LinearControlMap) = size(m.A, 1)
islinear(::LinearControlMap) = true
isaffine(::LinearControlMap) = true
apply(m::LinearControlMap, x, u) = m.A * x + m.B * u

"""
    ConstrainedLinearControlMap

A linear control map with state and input constraints,

```math
    (x, u) ↦ Ax + Bu, x ∈ \\mathcal{X}, u ∈ \\mathcal{U}.
```

### Fields

- `A` -- matrix
- `B` -- matrix
- `X` -- state constraints
- `U` -- input constraints
"""
struct ConstrainedLinearControlMap{T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST, UT} <: AbstractMap
    A::MTA
    B::MTB
    X::ST
    U::UT
    function ConstrainedLinearControlMap(A::MTA, B::MTB, X::ST, U::UT) where {T, MTA<:AbstractMatrix{T}, MTB<:AbstractMatrix{T}, ST, UT}
        @assert size(A, 1) == size(B, 1)
        return new{T, MTA, MTB, ST, UT}(A, B, X, U)
    end
end
statedim(m::ConstrainedLinearControlMap) = size(m.A, 2)
stateset(m::ConstrainedLinearControlMap) = m.X
outputdim(m::ConstrainedLinearControlMap) = size(m.A, 1)
inputdim(m::ConstrainedLinearControlMap) = size(m.B, 2)
inputset(m::ConstrainedLinearControlMap) = m.U
islinear(::ConstrainedLinearControlMap) = true
isaffine(::ConstrainedLinearControlMap) = true
apply(m::ConstrainedLinearControlMap, x, u) = m.A * x + m.B * u

"""
    AffineControlMap

An affine control map,

```math
    (x, u) ↦ Ax + Bu + c.
```

### Fields

- `A` -- matrix
- `B` -- matrix
- `c` -- vector
"""
struct AffineControlMap{T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT<:AbstractVector{T}} <: AbstractMap
    A::MTA
    B::MTB
    c::VT
    function AffineControlMap(A::MTA, B::MTB, c::VT) where {T, MTA<:AbstractMatrix{T}, MTB<:AbstractMatrix{T}, VT<:AbstractVector{T}}
        @assert size(A, 1) == size(B, 1) == length(c)
        return new{T, MTA, MTB, VT}(A, B, c)
    end
end
statedim(m::AffineControlMap) = size(m.A, 2)
outputdim(m::AffineControlMap) = size(m.A, 1)
inputdim(m::AffineControlMap) = size(m.B, 1)
islinear(::AffineControlMap) = false
isaffine(::AffineControlMap) = true
apply(m::AffineControlMap, x, u) = m.A * x + m.B * u + m.c

"""
    ConstrainedAffineControlMap

An affine control map with state and input constraints,

```math
    (x, u) ↦ Ax + Bu + c, x ∈ \\mathcal{X}, u ∈ \\mathcal{U}.
```

### Fields

- `A` -- matrix
- `B` -- matrix
- `c` -- vector
- `X` -- state constraints
- `U` -- input constraints
"""
struct ConstrainedAffineControlMap{T, MTA<:AbstractMatrix{T}, MTB<:AbstractMatrix{T}, VT<:AbstractVector{T}, ST, UT} <: AbstractMap
    A::MTA
    B::MTB
    c::VT
    X::ST
    U::UT
    function ConstrainedAffineControlMap(A::MTA, B::MTB, c::VT, X::ST, U::UT) where {T, MTA<:AbstractMatrix{T}, MTB<:AbstractMatrix{T}, VT<:AbstractVector{T}, ST, UT}
        @assert size(A, 1) == size(B, 1) == length(c)
        return new{T, MTA, MTB, VT, ST, UT}(A, B, c, X, U)
    end
end
statedim(m::ConstrainedAffineControlMap) = size(m.A, 2)
stateset(m::ConstrainedAffineControlMap) = m.X
inputdim(m::ConstrainedAffineControlMap) = size(m.B, 2)
inputset(m::ConstrainedAffineControlMap) = m.U
outputdim(m::ConstrainedAffineControlMap) = size(m.A, 1)
islinear(::ConstrainedAffineControlMap) = false
isaffine(::ConstrainedAffineControlMap) = true
apply(m::ConstrainedAffineControlMap, x, u) = m.A * x + m.B * u + m.c

"""
    ResetMap

A reset map,

```math
    x ↦ R(x),
```
such that a subset of the variables is given a specified value, and the rest
are unchanged.

### Fields

- `dim`  -- dimension
- `dict` -- dictionary whose keys are the indices of the variables that are reset,
            and whose values are the new values
"""
struct ResetMap{N} <: AbstractMap
    dim::Int
    dict::Dict{Int, N}
end
statedim(m::ResetMap) = m.dim
inputdim(::ResetMap) = 0
outputdim(m::ResetMap) = m.dim
islinear(::ResetMap) = false
isaffine(::ResetMap) = true

# convenience constructor for a list of pairs instead of a dictionary
ResetMap(dim::Int, args::Pair{Int, <:N}...) where {N} = ResetMap(dim, Dict{Int, N}(args))

"""
    ConstrainedResetMap

A reset map with state constraints of the form:

```math
    x ↦ R(x), x ∈ \\mathcal{X},
```
such that the specified variables are assigned a given value, and the remaining
variables are unchanged.

### Fields

- `dim`  -- dimension
- `X`    -- state constraints
- `dict` -- dictionary whose keys are the indices of the variables that are
            reset, and whose values are the new values
"""
struct ConstrainedResetMap{N, ST} <: AbstractMap
    dim::Int
    X::ST
    dict::Dict{Int, N}
end
statedim(m::ConstrainedResetMap) = m.dim
stateset(m::ConstrainedResetMap) = m.X
inputdim(::ConstrainedResetMap) = 0
outputdim(m::ConstrainedResetMap) = m.dim
islinear(::ConstrainedResetMap) = false
isaffine(::ConstrainedResetMap) = true

# convenience constructor for a list of pairs instead of a dictionary
ConstrainedResetMap(dim::Int, X::ST, args::Pair{Int, <:N}...) where {N, ST} =
    ConstrainedResetMap(dim, X, Dict{Int, N}(args))

function apply(m::Union{ResetMap, ConstrainedResetMap}, x)
    y = copy(x)
    for (index, value) in pairs(m.dict)
        y[index] = value
    end
    return y
end


text:
"""
    AdiabaticASEModel{A} <: AdiabaticModel

Wrapper for an `ase.Atoms` object that has a calculator attached.
This will synchronise the positions with the `ase` object and handle the unit conversions.

Implements both `potential` and `derivative!`.
"""
struct AdiabaticASEModel{A} <: AdiabaticModel
    atoms::A
end

NQCModels.ndofs(::AdiabaticASEModel) = 3

function NQCModels.potential(model::AdiabaticASEModel, R::AbstractMatrix)
    set_coordinates!(model, R)
    V = model.atoms.get_potential_energy()
    return austrip(V * u"eV")
end

function NQCModels.derivative!(model::AdiabaticASEModel, D::AbstractMatrix, R::AbstractMatrix)
    set_coordinates!(model, R)
    D .= -model.atoms.get_forces()'
    @. D = austrip(D * u"eV/Å")
    return D
end

function set_coordinates!(model::AdiabaticASEModel, R)
    model.atoms.set_positions(ustrip.(auconvert.(u"Å", R')))
end


text:
#Get displacement at a single point over the time of the simulation
function get_field_point(fname::String, x_read::Float64 , y_read::Float64)

    #Get Variables
    fid         = h5open( fname , "r")

    max_iter = read_attribute(fid, "max_iter")

    xmin = read_attribute(fid, "xmin")
    xmax = read_attribute(fid, "xmax")
    Nx   = read_attribute(fid, "xnodes")

    ymin = read_attribute(fid, "ymin")
    ymax = read_attribute(fid, "ymax")
    Ny   = read_attribute(fid, "ynodes")

    #Get associated index
    xi_read =floor(Int , Nx * (x_read - xmin) / (xmax - xmin))
    yi_read =floor(Int , Ny * (y_read - ymin) / (ymax - ymin))

    close(fid)

    #Check bounds
    if(xi_read > Nx || yi_read > Ny)
        println("ERROR: Index out of bounds. Max index is " , Nx)
        return
    end

    #Max iteration time
    _ , _ , max_sim_iter = get_time_variables(fname)

    #Storage vectors
    pos_vector = Float64[]
    dpos_vector = Float64[]
    time_vector = Float64[]

    #ReadOnce
    t, x , y , ψ , dψ = get_fields(fname,1)

    #Cycle
    for i in 1:max_sim_iter

        #Time and fields
        t, _ , _ , ψ , dψ = get_fields(fname,i)

        push!(time_vector, t)
        push!(pos_vector, ψ[xi_read , yi_read] )
        push!(dpos_vector, dψ[xi_read , yi_read] )

    end

    #return vectors
    return time_vector , pos_vector , dpos_vector

end


text:
using Trilinos
using Test
using MPI

pl = Teuchos.ParameterList("test1")
@test Teuchos.name(pl) == "test1"
@test isempty(pl)

Teuchos.setName(pl, "test2")
@test Teuchos.name(pl) == "test2"

@test Teuchos.numParams(pl) == 0
@test !haskey(pl, "intparam")

Teuchos.set(pl, "intparam", 1)
@test haskey(pl, "intparam")
@test Teuchos.isParameter(pl, "intparam")
@test !Teuchos.isParameter(pl, "dummy")
@test Teuchos.isType(Int, pl, "intparam")
@test !Teuchos.isType(Int, pl, "dummy")
@test_throws ErrorException Teuchos.get(Int, pl, "dummy")
@test Teuchos.get(Int, pl, "intparam") == 1
@test Teuchos.get_type(pl, "intparam") == Int
@test pl["intparam"] == 1

pl["intparam"] = 2
@test pl["intparam"] == 2

pl["dblparam"] = 2.0
@test length(pl) == 2
@test pl["dblparam"] == 2.0

@test_throws KeyError pl["dummy"]

pl2 = Teuchos.sublist(pl, "child")
pl2["param1"] = "test"

pl3 = Teuchos.ParameterList()
pl2["Manual Sublist"] = pl3
@test Teuchos.name(pl2["Manual Sublist"]) == "ANONYMOUS"

@test haskey(pl, "child")
@test Teuchos.isSublist(pl, "child")
@test pl["child"]["param1"] == "test"

pl2_ref = pl["child"]
pl2_ref["param1"] = "test2"

@test pl["child"]["param1"] == "test2"

@test get(pl, "intparam", 0) == 2
@test get(pl, "thisdoesnotexist", "somedefault") == "somedefault"

pl_copy = Teuchos.ParameterList(pl)
pl_copy["intparam"] = 3
@test pl["intparam"] == 2
@test pl_copy["intparam"] == 3
pl_copy["child"]["param1"] = "copy_modified"
@test pl["child"]["param1"] == "test2"
@test pl_copy["child"]["param1"] == "copy_modified"

pl["enum"] = Teuchos.VERB_DEFAULT
@test pl["enum"] == -1

display(pl[])
println()


text:
using PlotlyJS

function get_plotly_data(E::Vector{Tuple},
                   xcoord::Vector{T}, ycoord::Vector{T}) where T <: Real

    # E: Vector of tuples representing the graph edges
    # xcoords, ycoords: Vector(s) of node coordinates returned by a graph layout

    Xedges = Float64[]
    Yedges = Float64[]
    for e in E
        append!(Xedges, [xcoord[e[1]], xcoord[e[2]], NaN])
        append!(Yedges, [ycoord[e[1]], ycoord[e[2]], NaN])
    end
    return Xedges, Yedges
end;

function get_node_trace(x::Vector{T}, y::Vector{T}, labels::Vector{String};
                   marker_size=5, marker_color="#6959CD",
                   linecolor="rgb(50,50,50)", linewidth=0.5) where T <: Real

    return PlotlyJS.scatter(
                x=x,
                y=y,
                mode="markers",
                marker=attr(
                            size=marker_size,
                            color=marker_color,
                            line=attr(color=linecolor, width=linewidth)
                             ),
                text=labels,
                hoverinfo="text")
end;

function get_edge_trace(xedge::Vector{T}, yedge::Vector{T};
                 linecolor="rgb(210,210,210)", linewidth=1) where T <: Real
    return PlotlyJS.scatter(
                x=xedge,
                y=yedge,
                mode="lines",
                line_color=linecolor,
                line_width=linewidth,
                hoverinfo="none"
               )
end;

function perfect_balanced_tree_edges(;bn=2, h=6)

    # bn:  branching number; each node has bn children
    # h: the tree height, i.e. each leave is at the distance h of the tree root
    # returns the edges of such a tree

    if bn == 1
        n = h+1
    else
       n = (bn^(h + 1)-1) ÷ (bn -1)
    end

    nodes = 1:n
    k = 1
    parents = [nodes[k]]
    edges = Tuple[]

    while !isempty(parents)
        s = popfirst!(parents)
        for _ in 1:bn
            k += 1
            if k <= n
                push!(parents, nodes[k])
                push!(edges, (s, k))
            else
                break
            end
        end
    end
    return edges
end;

function  Bezier_edges(xc::Vector{T}, yc::Vector{T},
                edges::Vector{Tuple{Int, Int}};
                linewidth=1.5, linecolor="rgb(210, 210, 210)") where T <: Real

    shapes =  []
    for e in edges[1:end]
        #set the coordinates of the Bezier control points d_i(xi, yi),
        #for the cubic Bezier representing a tree edge
        x1 = xc[e[1]]
        y1 = yc[e[1]]
        d1 = sqrt(x1^2 + y1^2)
        x4 = xc[e[2]]
        y4 = yc[e[2]]
        d4 = sqrt(x4^2 + y4^2)
        d = 0.5*(d1+d4)
        if d1 == 0 || d4 == 0 # if the first Bezier  control point is  the  root, def edge as a segment of line
            x2 = x1
            y2 = y1
            x3 = x4
            y3 = 0.3*y1 + 0.7*y4
        else
            x2 = d*x1/d1
            y2 = d*y1/d1
            x3 = d*x4/d4
            y3 = d*y4/d4
        end
        push!(shapes, attr(type="path",
                           layer="below",
                           path= "M $x1 $y1, C $x2 $y2, $x3 $y3, $x4  $y4",
                           line=attr(color=linecolor, width=linewidth)
                                ))
    end
    return shapes
end;

function proj2circle(xpos::Vector{T}, ypos::Vector{T};
                dirfactor=-1,  θ₀=-π/2) where T <: Real

    xmin, xmax = extrema(xpos)
    ymin, ymax = extrema(ypos)
    θ = 2π*(xpos .- xmin)/(xmax-xmin) .+ θ₀
    r = ymax .-  ypos
    xc = r .* cos.(θ)
    yc = dirfactor * r .* sin.(θ)
    return xc, yc
end;

max_stars_repo_name: empet/Julia-Networks

max_issues_repo_name: empet/Julia-Networks

max_forks_repo_name: empet/Julia-Networks


text:
export DQNLearner

using Setfield: @set
using Random: AbstractRNG, GLOBAL_RNG
import Functors

Base.@kwdef mutable struct DQNLearner{A<:Approximator{<:TwinNetwork}} <: AbstractLearner
    approximator::A
    loss_func::Any
    n::Int = 1
    γ::Float32 = 0.99f0
    is_enable_double_DQN::Bool = true
    rng::AbstractRNG = GLOBAL_RNG
    # for logging
    loss::Float32 = 0.0f0
end

(L::DQNLearner)(s::AbstractArray) = L.approximator(s)

Functors.functor(x::DQNLearner) = (; approximator=x.approximator), y -> @set x.approximator = y.approximator

function RLBase.optimise!(learner::DQNLearner, batch::Union{NamedTuple{SS′ART},NamedTuple{SS′L′ART}})
    A = learner.approximator
    Q = A.model.source
    Qₜ = A.model.target
    γ = learner.γ
    loss_func = learner.loss_func
    n = learner.n

    s, s′, a, r, t = map(x -> batch[x], SS′ART)
    a = CartesianIndex.(a, 1:length(a))

    q′ = learner.is_enable_double_DQN ? Q(s′) : Qₜ(s′)

    if haskey(batch, :next_legal_actions_mask)
        q′ .+= ifelse.(batch[:next_legal_actions_mask], 0.0f0, typemin(Float32))
    end

    q′ₐ = learner.is_enable_double_DQN ? Qₜ(s′)[dropdims(argmax(q′, dims=1), dims=1)] : dropdims(maximum(q′; dims=1), dims=1)

    G = r .+ γ^n .* (1 .- t) .* q′ₐ

    gs = gradient(params(A)) do
        qₐ = Q(s)[a]
        loss = loss_func(G, qₐ)
        ignore() do
            learner.loss = loss
        end
        loss
    end

    optimise!(A, gs)
end