text:
###############################################################################
#
#   fq.jl : Flint finite fields
#
###############################################################################

export FlintFiniteField, characteristic, order, fq, FqFiniteField, frobenius,
       pth_root, trace, norm

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent_type(::Type{fq}) = FqFiniteField

elem_type(::Type{FqFiniteField}) = fq

doc"""
    base_ring(a::FqFiniteField)
> Returns `Union{}` as this field is not dependent on another field.
"""
base_ring(a::FqFiniteField) = Union{}

doc"""
    base_ring(a::fq)
> Returns `Union{}` as this field is not dependent on another field.
"""
base_ring(a::fq) = Union{}

doc"""
    parent(a::fq)
> Returns the parent of the given finite field element.
"""
parent(a::fq) = a.parent

isdomain_type(::Type{fq}) = true

function check_parent(a::fq, b::fq)
   a.parent != b.parent && error("Operations on distinct finite fields not supported")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function Base.hash(a::fq, h::UInt)
   b = 0xb310fb6ea97e1f1a%UInt
   for i in 1:degree(parent(a)) + 1
      b = xor(b, xor(hash(coeff(a, i), h), h))
      b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
   end
   return b
end

doc"""
    coeff(x::fq, n::Int)
> Return the degree $n$ coefficient of the polynomial representing the given
> finite field element.
"""
function coeff(x::fq, n::Int)
   n < 0 && throw(DomainError())
   z = fmpz()
   ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void,
               (Ref{fmpz}, Ref{fq}, Int), z, x, n)
   return z
end

doc"""
    zero(a::FqFiniteField)
> Return the additive identity, zero, in the given finite field.
"""
function zero(a::FqFiniteField)
   d = a()
   ccall((:fq_zero, :libflint), Void, (Ref{fq}, Ref{FqFiniteField}), d, a)
   return d
end

doc"""
    one(a::FqFiniteField)
> Return the multiplicative identity, one, in the given finite field.
"""
function one(a::FqFiniteField)
   d = a()
   ccall((:fq_one, :libflint), Void, (Ref{fq}, Ref{FqFiniteField}), d, a)
   return d
end

doc"""
    gen(a::FqFiniteField)
> Return the generator of the finite field. Note that this is only guaranteed
> to be a multiplicative generator if the finite field is generated by a
> Conway polynomial automatically.
"""
function gen(a::FqFiniteField)
   d = a()
   ccall((:fq_gen, :libflint), Void, (Ref{fq}, Ref{FqFiniteField}), d, a)
   return d
end

doc"""
    iszero(a::fq)
> Return `true` if the given finite field element is zero, otherwise return
> `false`.
"""
iszero(a::fq) = ccall((:fq_is_zero, :libflint), Bool,
                     (Ref{fq}, Ref{FqFiniteField}), a, a.parent)

doc"""
    isone(a::fq)
> Return `true` if the given finite field element is one, otherwise return
> `false`.
"""
isone(a::fq) = ccall((:fq_is_one, :libflint), Bool,
                    (Ref{fq}, Ref{FqFiniteField}), a, a.parent)

doc"""
    isgen(a::fq)
> Return `true` if the given finite field element is the generator of the
> finite field, otherwise return `false`.
"""
isgen(a::fq) = a == gen(parent(a))

doc"""
    isunit(a::fq)
> Return `true` if the given finite field element is invertible, i.e. nonzero,
> otherwise return `false`.
"""
isunit(a::fq) = ccall((:fq_is_invertible, :libflint), Bool,
                     (Ref{fq}, Ref{FqFiniteField}), a, a.parent)

doc"""
    characteristic(a::FqFiniteField)
> Return the characteristic of the given finite field.
"""
function characteristic(a::FqFiniteField)
   d = fmpz()
   ccall((:__fq_ctx_prime, :libflint), Void,
         (Ref{fmpz}, Ref{FqFiniteField}), d, a)
   return d
end

doc"""
    order(a::FqFiniteField)
> Return the order, i.e. the number of elements in, the given finite field.
"""
function order(a::FqFiniteField)
   d = fmpz()
   ccall((:fq_ctx_order, :libflint), Void,
         (Ref{fmpz}, Ref{FqFiniteField}), d, a)
   return d
end

doc"""
    degree(a::FqFiniteField)
> Return the degree of the given finite field.
"""
function degree(a::FqFiniteField)
   return ccall((:fq_ctx_degree, :libflint), Int, (Ref{FqFiniteField},), a)
end

function deepcopy_internal(d::fq, dict::ObjectIdDict)
   z = fq(parent(d), d)
   return z
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(x::fq) = x

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::fq)
   cstr = ccall((:fq_get_str_pretty, :libflint), Ref{UInt8},
                (Ref{fq}, Ref{FqFiniteField}), x, x.parent)

   print(io, unsafe_string(cstr))

   ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
end

function show(io::IO, a::FqFiniteField)
   print(io, "Finite field of degree ", degree(a))
   print(io, " over F_", characteristic(a))
end

needs_parentheses(x::fq) = x.length > 1

isnegative(x::fq) = false

show_minus_one(::Type{fq}) = true

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

function -(x::fq)
   z = parent(x)()
   ccall((:fq_neg, :libflint), Void,
         (Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, x.parent)
   return z
end

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

function +(x::fq, y::fq)
   check_parent(x, y)
   z = parent(y)()
   ccall((:fq_add, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, y, y.parent)
   return z
end

function -(x::fq, y::fq)
   check_parent(x, y)
   z = parent(y)()
   ccall((:fq_sub, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, y, y.parent)
   return z
end

function *(x::fq, y::fq)
   check_parent(x, y)
   z = parent(y)()
   ccall((:fq_mul, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, y, y.parent)
   return z
end

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

function *(x::Int, y::fq)
   z = parent(y)()
   ccall((:fq_mul_si, :libflint), Void,
         (Ref{fq}, Ref{fq}, Int, Ref{FqFiniteField}), z, y, x, y.parent)
   return z
end

*(x::Integer, y::fq) = fmpz(x)*y

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

function *(x::fmpz, y::fq)
   z = parent(y)()
   ccall((:fq_mul_fmpz, :libflint), Void,
         (Ref{fq}, Ref{fq}, Ref{fmpz}, Ref{FqFiniteField}),
                                            z, y, x, y.parent)
   return z
end

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

+(x::fq, y::Integer) = x + parent(x)(y)

+(x::Integer, y::fq) = y + x

+(x::fq, y::fmpz) = x + parent(x)(y)

+(x::fmpz, y::fq) = y + x

-(x::fq, y::Integer) = x - parent(x)(y)

-(x::Integer, y::fq) = parent(y)(x) - y

-(x::fq, y::fmpz) = x - parent(x)(y)

-(x::fmpz, y::fq) = parent(y)(x) - y

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

function ^(x::fq, y::Int)
   if y < 0
      x = inv(x)
      y = -y
   end
   z = parent(x)()
   ccall((:fq_pow_ui, :libflint), Void,
         (Ref{fq}, Ref{fq}, Int, Ref{FqFiniteField}), z, x, y, x.parent)
   return z
end

function ^(x::fq, y::fmpz)
   if y < 0
      x = inv(x)
      y = -y
   end
   z = parent(x)()
   ccall((:fq_pow, :libflint), Void,
         (Ref{fq}, Ref{fq}, Ref{fmpz}, Ref{FqFiniteField}),
                                            z, x, y, x.parent)
   return z
end

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

function ==(x::fq, y::fq)
   check_parent(x, y)
   ccall((:fq_equal, :libflint), Bool,
         (Ref{fq}, Ref{fq}, Ref{FqFiniteField}), x, y, y.parent)
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

==(x::fq, y::Integer) = x == parent(x)(y)

==(x::fq, y::fmpz) = x == parent(x)(y)

==(x::Integer, y::fq) = parent(y)(x) == y

==(x::fmpz, y::fq) = parent(y)(x) == y

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fq, y::fq)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   z = parent(y)()
   ccall((:fq_div, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, y, y.parent)
   return z
end

function divides(a::fq, b::fq)
   iszero(b) && error("Division by zero in divides")
   return true, divexact(a, b)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

divexact(x::fq, y::Integer) = divexact(x, parent(x)(y))

divexact(x::fq, y::fmpz) = divexact(x, parent(x)(y))

divexact(x::Integer, y::fq) = divexact(parent(y)(x), y)

divexact(x::fmpz, y::fq) = divexact(parent(y)(x), y)

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

doc"""
    inv(x::fq)
> Return $x^{-1}$.
"""
function inv(x::fq)
   iszero(x) && throw(DivideError())
   z = parent(x)()
   ccall((:fq_inv, :libflint), Void,
         (Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, x.parent)
   return z
end

###############################################################################
#
#   Special functions
#
###############################################################################

doc"""
    pth_root(x::fq)
> Return the $p$-th root of $a$ in the finite field of characteristic $p$. This
> is the inverse operation to the Frobenius map $\sigma_p$.
"""
function pth_root(x::fq)
   z = parent(x)()
   ccall((:fq_pth_root, :libflint), Void,
         (Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, x.parent)
   return z
end

doc"""
    trace(x::fq)
> Return the trace of $a$. This is an element of $\F_p$, but the value returned
> is this value embedded in the original finite field.
"""
function trace(x::fq)
   z = fmpz()
   ccall((:fq_trace, :libflint), Void,
         (Ref{fmpz}, Ref{fq}, Ref{FqFiniteField}), z, x, x.parent)
   return parent(x)(z)
end

doc"""
    norm(x::fq)
> Return the norm of $a$. This is an element of $\F_p$, but the value returned
> is this value embedded in the original finite field.
"""
function norm(x::fq)
   z = fmpz()
   ccall((:fq_norm, :libflint), Void,
         (Ref{fmpz}, Ref{fq}, Ref{FqFiniteField}), z, x, x.parent)
   return parent(x)(z)
end

doc"""
    frobenius(x::fq, n = 1)
> Return the iterated Frobenius $\sigma_p^n(a)$ where $\sigma_p$ is the
> Frobenius map sending the element $a$ to $a^p$ in the finite field of
> characteristic $p$. By default the Frobenius map is applied $n = 1$ times if
> $n$ is not specified.
"""
function frobenius(x::fq, n = 1)
   z = parent(x)()
   ccall((:fq_frobenius, :libflint), Void,
         (Ref{fq}, Ref{fq}, Int, Ref{FqFiniteField}), z, x, n, x.parent)
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::fq)
   ccall((:fq_zero, :libflint), Void,
        (Ref{fq}, Ref{FqFiniteField}), z, z.parent)
   return z
end

function mul!(z::fq, x::fq, y::fq)
   ccall((:fq_mul, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, y, y.parent)
   return z
end

function addeq!(z::fq, x::fq)
   ccall((:fq_add, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, z, x, x.parent)
   return z
end

function add!(z::fq, x::fq, y::fq)
   ccall((:fq_add, :libflint), Void,
        (Ref{fq}, Ref{fq}, Ref{fq}, Ref{FqFiniteField}), z, x, y, x.parent)
   return z
end

###############################################################################
#
#   Random functions
#
###############################################################################

function rand(K::FinField)
	p = characteristic(K)
	r = degree(K)
	alpha = gen(K)
	res = zero(K)
	for i = 0 : (r-1)
		c = rand(BigInt(0) : BigInt(p - 1))
		res += c * alpha^i
	end
	return res
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{fq}, ::Type{T}) where {T <: Integer} = fq

promote_rule(::Type{fq}, ::Type{fmpz}) = fq

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FqFiniteField)()
   z = fq(a)
   return z
end

(a::FqFiniteField)(b::Integer) = a(fmpz(b))

function (a::FqFiniteField)(b::Int)
   z = fq(a, b)
   return z
end

function (a::FqFiniteField)(b::fmpz)
   z = fq(a, b)
   return z
end

function (a::FqFiniteField)(b::fq)
   parent(b) != a && error("Coercion between finite fields not implemented")
   return b
end

###############################################################################
#
#   FlintFiniteField constructor
#
###############################################################################

doc"""
    FlintFiniteField(char::fmpz, deg::Int, s::AbstractString)
> Returns a tuple $S, x$ consisting of a finite field parent object $S$ and
> generator $x$ for the finite field of the given characteristic and degree.
> The string $s$ is used to designate how the finite field generator will be
> printed. The characteristic must be prime. When a Conway polynomial is known,
> the field is generated using the Conway polynomial. Otherwise a random
> sparse, irreducible polynomial is used. The generator of the field is
> guaranteed to be a multiplicative generator only if the field is generated by
> a Conway polynomial. We require the degree to be positive.
"""
function FlintFiniteField(char::fmpz, deg::Int, s::AbstractString; cached = true)
   S = Symbol(s)
   parent_obj = FqFiniteField(char, deg, S, cached)
   return parent_obj, gen(parent_obj)
end

doc"""
    FlintFiniteField(char::Integer, deg::Int, s::AbstractString)
> Returns a tuple $S, x$ consisting of a finite field parent object $S$ and
> generator $x$ for the finite field of the given characteristic and degree.
> The string $s$ is used to designate how the finite field generator will be
> printed. The characteristic must be prime. When a Conway polynomial is known,
> the field is generated using the Conway polynomial. Otherwise a random
> sparse, irreducible polynomial is used. The generator of the field is
> guaranteed to be a multiplicative generator only if the field is generated by
> a Conway polynomial. We require the degree to be positive.
"""
function FlintFiniteField(char::Integer, deg::Int, s::AbstractString; cached = true)
   return FlintFiniteField(fmpz(char), deg, s; cached = cached)
end

doc"""
    FlintFiniteField(pol::fmpz_mod_poly, s::AbstractString)
> Returns a tuple $S, x$ consisting of a finite field parent object $S$ and
> generator $x$ for the finite field over $F_p$ defined by the given
> polynomial, i.e. $\mathbb{F}_p[t]/(pol)$. The characteristic is specified by
> the modulus of `pol`. The polynomial is required to be irreducible, but this
> is not checked. The string $s$ is used to designate how the finite field
> generator will be printed. The generator will not be multiplicative in
> general.
"""
function FlintFiniteField(pol::fmpz_mod_poly, s::AbstractString; cached = true)
   S = Symbol(s)
   parent_obj = FqFiniteField(pol, S, cached)

   return parent_obj, gen(parent_obj)
end


text:
using CUDArt,CUDNN,KUDNN,Base.Test
include("isapprox.jl")

cd4 = ConvolutionDescriptor(padding=(0,0), stride=(1,1), upscale=(1,1), mode=CUDNN_CROSS_CORRELATION)
x1 = CudaArray(rand(10,20,3,4)) # W,H,C,N
w1 = CudaArray(rand(5,4,3,2))   # X1,Y,C,K
y1 = kudnnConvolutionForward(x1, w1; convDesc=cd4)
y2 = cudnnConvolutionForward(x1, w1; convDesc=cd4)
@show map(size, (x1, w1, y1,y2))
@show to_host(y1)==to_host(y2)
@test @show isapprox(y1,y2)
@test @show size(y1) == kudnnGetConvolutionNdForwardOutputDim(x1, w1; convDesc=cd4)

dy1 = CudaArray(rand(size(y1)...))
dw1 = kudnnConvolutionBackwardFilter(x1, dy1, copy(w1); convDesc=cd4)
dw2 = cudnnConvolutionBackwardFilter(x1, dy1, copy(w1); convDesc=cd4)
@show map(size, (x1, w1, dy1, dw1, dw2))
@show to_host(dw1)==to_host(dw2)
@test @show isapprox(dw1,dw2)

dx1 = kudnnConvolutionBackwardData(w1, dy1, copy(x1); convDesc=cd4)
dx2 = cudnnConvolutionBackwardData(w1, dy1, copy(x1); convDesc=cd4)
@show map(size, (x1, w1, dy1, dx1, dx2))
@show to_host(dx1)==to_host(dx2)
@test @show isapprox(dx1,dx2)

cd5 = ConvolutionDescriptor(padding=(0,0,0), stride=(1,1,1), upscale=(1,1,1), mode=CUDNN_CROSS_CORRELATION)
x2 = CudaArray(rand(10,20,30,3,4)) # 1,2,3,C,N
w2 = CudaArray(rand(6,5,4,3,2))   # 1,2,3,C,K
y3 = kudnnConvolutionForward(x2, w2; convDesc=cd5)
@show map(size, (x2, w2, y3))
@test @show size(y3) == kudnnGetConvolutionNdForwardOutputDim(x2, w2; convDesc=cd5)
# y4 = cudnnConvolutionForward(x2, w2; convDesc=cd5)
# @show isapprox(y3,y4)
# @show to_host(y3)==to_host(y4)

dy2 = CudaArray(rand(size(y3)...))
dw3 = kudnnConvolutionBackwardFilter(x2, dy2, copy(w2); convDesc=cd5)
@show map(size, (x2, w2, dy2, dw3))
# dw4 = cudnnConvolutionBackwardFilter(x2, dy2, copy(w2); convDesc=cd5)
# @show isapprox(dw3,dw4)
# @show to_host(dw3)==to_host(dw4)

dx3 = kudnnConvolutionBackwardData(w2, dy2, copy(x2); convDesc=cd5)
@show map(size, (x2, w2, dy2, dx3))
# dx4 = cudnnConvolutionBackwardData(w2, dy2, copy(x2); convDesc=cd5)
# @show isapprox(dx3,dx4)
# @show to_host(dx3)==to_host(dx4)

pd4 = PoolingDescriptor((3,2))
x10 = CudaArray(rand(5,4,3,2)) # W,H,C,N
y10 = similar(x10, kudnnGetPoolingNdForwardOutputDim(pd4,x10))
y11 = kudnnPoolingForward(pd4, x10, copy(y10))
y12 = cudnnPoolingForward(pd4, x10, copy(y10))
@show map(size, (x10,y10))
@show to_host(y11)==to_host(y12)
@test @show isapprox(y11,y12)

dy10 = CudaArray(rand(size(y10)...))
dx10 = similar(x10)
dx11 = kudnnPoolingBackward(pd4, y11, dy10, x10, copy(dx10))
dx12 = cudnnPoolingBackward(pd4, y11, dy10, x10, copy(dx10))
@show to_host(dx11) == to_host(dx12)
@test @show isapprox(dx11,dx12)

pd5 = PoolingDescriptor((3,4,5))
x13 = CudaArray(rand(6,5,4,3,2)) # W,H,C,N
y13 = similar(x13, kudnnGetPoolingNdForwardOutputDim(pd5,x13))
y14 = kudnnPoolingForward(pd5, x13, copy(y13))
# y15 = cudnnPoolingForward(pd4, x10, copy(y10))
@show map(size, (x13,y14))
# @show to_host(y11)==to_host(y12)
# @test @show isapprox(y11,y12)

dy13 = CudaArray(rand(size(y13)...))
dx13 = similar(x13)
dx16 = kudnnPoolingBackward(pd5, y14, dy13, x13, copy(dx13))
# dx17 = cudnnPoolingBackward(pd5, y15, dy13, x13, copy(dx13))
# @show to_host(dx11) == to_host(dx12)
# @test @show isapprox(dx11,dx12)

:ok

ext: jl

max_stars_repo_name: denizyuret/KUDNN

max_issues_repo_name: denizyuret/KUDNN

max_forks_repo_name: denizyuret/KUDNN


text:
# Issue #14

using JuliaInterpreter: JuliaInterpreter, pc_expr, evaluate_call!, finish_and_return!, @lookup, enter_call_expr
runframe(frame::Frame, pc=frame.pc[]) = Some{Any}(finish_and_return!(Compiled(), frame))

frame = @make_frame map(x->2x, 1:10)
state = dummy_state(frame)
execute_command(state, Val{:so}(), "so")
@test isnothing(state.frame)
@test state.overall_result == 2 .* [1:10...]

# Issue #12
function complicated_keyword_stuff(args...; kw...)
    args[1] == args[1]
    (args..., kw...)
end
frame = @make_frame complicated_keyword_stuff(1)
state = dummy_state(frame)
execute_command(state, Val{:n}(), "n")
execute_command(state, Val{:so}(), "so")
@test isnothing(state.frame)

@test runframe(JuliaInterpreter.enter_call(complicated_keyword_stuff, 1, 2)) ==
      runframe(@make_frame(complicated_keyword_stuff(1, 2)))
@test runframe(JuliaInterpreter.enter_call(complicated_keyword_stuff, 1, 2; x=7, y=33)) ==
      runframe(@make_frame(complicated_keyword_stuff(1, 2; x=7, y=33)))

# Issue #22
f22() = string(:(a+b))
@test step_through(enter_call_expr(:($f22()))) == "a + b"
f22() = string(QuoteNode(:a))
@test step_through(enter_call_expr(:($f22()))) == ":a"

# Breakpoints
f_outer_break(x) = f_inner_break(x)
f_inner_break(x) = x
JuliaInterpreter.breakpoint(f_inner_break)
frame = @make_frame f_outer_break(2)
state = dummy_state(frame)
execute_command(state, Val{:so}(), "c")
@test state.frame.framecode.scope.name === :f_inner_break
execute_command(state, Val{:so}(), "c")
@test state.frame === nothing
@test state.overall_result == 2

@inline fnothing(x) = 1
frame = @make_frame fnothing(0)
io = IOBuffer()
Debugger.print_next_expr(io, frame)
@test chomp(String(take!(io))) == "About to run: return 1"


text:
module NetworkPowerSwitch
using Sockets
using Libaudio

# Network Controlled Relay IP: 192.168.1.199:12345 (client)
# Host IP: 192.168.1.190:6000 (server)

"""
    server()

deprecated. use client function for better responsiveness.
"""
function server()
    op = Array{String,1}()
    #@async begin
        server = listen(IPv4(0), 6000)
        #while true
            sock = accept(server)
            write(sock, [0x41 0x54 0x0D 0x0A])
            status = readline(sock)
            push!(op, "AT: $status")

            write(sock, [0x41 0x54 0x2B 0x4C 0x49 0x4E 0x4B 0x53 0x54 0x41 0x54 0x3D 0x3F 0x0D 0x0A])
            status = readline(sock)
            push!(op, "AT+LINKSTAT: $status")
            write(sock, [0x41 0x54 0x2B 0x4D 0x4F 0x44 0x45 0x4C 0x3D 0x3F 0x0D 0x0A])
            status = readline(sock)
            push!(op, "AT+MODEL: $status")
            strtemp = "Init Status: "
            write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x30 0x3D 0x3F 0x0D 0x0A])
            for i = 1:4
                status = readline(sock)
                strtemp = strtemp * status
            end
            push!(op, strtemp)
            write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x30 0x3D 0x30 0x0D 0x0A])
            status = readline(sock)
            push!(op, "Turn off all switches: $status")
            sleep(5)
            write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x31 0x3D 0x31 0x0D 0x0A])
            status = readline(sock)
            push!(op, "Turn on switch 1: $status")
            sleep(20)
            write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x32 0x3D 0x31 0x0D 0x0A])
            status = readline(sock)
            push!(op, "Turn on switch 2: $status")
            strtemp = "Status: "
            write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x30 0x3D 0x3F 0x0D 0x0A])
            for i = 1:4
                status = readline(sock)
                strtemp = strtemp * status
            end
            push!(op, strtemp)
            close(sock)
        #end
    #end
    return op
end

"""
    npsopen(ip::IPv4=ip"192.168.1.199", port=12345)

open socket to nps

# Arguments
    - 'ip': IPv4 address of the remote switch
    - 'port': port number of the remote switch
"""
function npsopen(ip::IPv4=ip"192.168.1.199", port=12345)
    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    sock = connect(ip, port)

    write(sock, [0x41 0x54 0x0D 0x0A])
    status = readline(sock)
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsopen: AT $(status)")

    write(sock, [0x41 0x54 0x2B 0x4C 0x49 0x4E 0x4B 0x53 0x54 0x41 0x54 0x3D 0x3F 0x0D 0x0A])
    status = readline(sock)
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsopen: AT+LINKSTAT $(status)")

    write(sock, [0x41 0x54 0x2B 0x4D 0x4F 0x44 0x45 0x4C 0x3D 0x3F 0x0D 0x0A])
    status = readline(sock)
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsopen: AT+MODEL: $(status)")

    strtemp = "Init Status: "
    write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x30 0x3D 0x3F 0x0D 0x0A])
    for i = 1:4
        status = readline(sock)
        strtemp = strtemp * status
    end
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsopen: $(strtemp)")
    return sock
end

"""
    npsoff(sock)

turn off all switches
"""
function npsoff(sock)
    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x30 0x3D 0x30 0x0D 0x0A])
    status = readline(sock)
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsoff: turn off all switches $(status)")
    return nothing
end

"""
    npson(sock, i)

turn on i-th switch
"""
function npson(sock, i=1)
    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    pa = [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48]
    pb = [0x3D 0x31 0x0D 0x0A]
    write(sock, [pa convert(UInt8,0x30+i) pb])
    status = readline(sock)
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npson: turn on switch $(i) $(status)")
    return nothing
end

"""
    npsclose(sock)

get the updated register states, close the socket and return true
"""
function npsclose(sock)
    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    strtemp = "Updated Status: "
    write(sock, [0x41 0x54 0x2B 0x53 0x54 0x41 0x43 0x48 0x30 0x3D 0x3F 0x0D 0x0A])
    for i = 1:4
        status = readline(sock)
        strtemp = strtemp * status
    end
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsclose: $(strtemp)")
    close(sock)
    Libaudio.printl(root, :green, Libaudio.nows() * " | NetworkPowerSwitch.npsclose: socket closed")
    return true
end

"""
    npserr()

if something is wrong return false
"""
function npserr()
    root = joinpath(Libaudio.folder(), Libaudio.logfile())
    Libaudio.printl(root, :light_red, Libaudio.nows() * " | NetworkPowerSwitch.npserr: unknown reason")
    return false
end

"""
    reboot(t, ip, port) where T <: Real

Turn off all switches and turn on each switch with time sequence specified.
# Arguments
    - 't': time between off and on
           or visually, (all switches off)-->t[1]-->(switch 1 on)-->t[2]-->(switch 2 on)-->t[3]...
"""
function reboot(t::Array{T,1}) where T <: Real
    try
        nps = npsopen()
        npsoff(nps)
        for (j,k) in enumerate(t)
            sleep(k)
            npson(nps, j)
        end
        npsclose(nps)
    catch
        npserr()
    end
end

"""
    alloff()

turn all switches off
"""
function alloff()
    try
        nps = npsopen()
        npsoff(nps)
        npsclose(nps)
    catch
        npserr()
    end
end

end # module


text:
## lattice.jl : unit testing for lattice.jl
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

@testset "LatticeRule32                " begin

lattice_rule = LatticeRule32([0x00000001,0x00000022], 2, 0x00000037)
point = get_point(lattice_rule, 2)
@test lattice_rule isa LatticeRule32
@test length(point) == 2
@test all(0 .≤ point .< 1)

z_file_1 = joinpath(dirname(pathof(MultilevelEstimators)), "generating_vectors", "K_3600_32.txt")
z_file_2 = joinpath(dirname(pathof(MultilevelEstimators)), "generating_vectors", "CKN_250_20.txt")
lattice_rule = LatticeRule32(z_file_1, 101)
point = get_point(lattice_rule, 100_000)
@test lattice_rule isa LatticeRule32
@test length(point) == 101
@test all(0 .≤ point .< 1)

function approx_π(n::Integer)
    lattice_rule = LatticeRule32(z_file_2, 2)
    count = 0
    for i in 0:n-1
        x = get_point(lattice_rule, i)
        if x[1]*x[1] + x[2]*x[2] < 1
            count += 1
        end
    end
    return 4*count/n
end

for (n, ϵ) in zip(2 .^(0:20), exp10.([0 0 0 0 0 0 1 1 2 2 2 3 3 3 3 3 3 4 4 5]))
    @test approx_π(n) ≈ π atol=ϵ
end

@test LatticeRule32(3600) isa LatticeRule32
@test_throws ArgumentError LatticeRule32(3601)
@test LatticeRule32(z_file_2, 250) isa LatticeRule32
@test_throws ArgumentError LatticeRule32(z_file_2, 251)
@test_throws ArgumentError LatticeRule32(-2)
@test LatticeRule32(z_file_2, 1) isa LatticeRule32
@test_throws ArgumentError LatticeRule32(0)
@test_throws ArgumentError LatticeRule32(z_file_1, 24, 4294967297)

end


text:
using GLFW
using Colors
using Glutils
using CSyntax
using ModernGL
using StaticArrays

macro ignore(args...) end

@ignore begin
   include("basic.jl")
end

if isempty(ARGS)
   error("Missing obrigatory argument <example>.

      Usage: julia renderer.jl <example>
   ")
end

example = first(ARGS)

title  = "Example"
width  = 800
height = 800
background = "#4e26ff"

const window = Ref{GLFW.Window}()

var"@size"(args...)   = :( collect(GLFW.GetFramebufferSize(window[])) )
var"@width"(args...)  = :( first(@size) )
var"@height"(args...) = :( last(@size) )
var"@window"(args...) = :( window[] )

GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 4)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 6)
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, true)
GLFW.WindowHint(GLFW.OPENGL_DEBUG_CONTEXT , true)
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)

window[] = GLFW.CreateWindow(width, height, title)
GLFW.MakeContextCurrent(window[])

ondebug_msg() do info::GLDebugInfo
   type, source, message, severity = info
   @error "$(message)" type source severity
end

include("$example.jl")

if @isdefined(vsync) && vsync
   GLFW.SwapInterval(1)
end

GLFW.SetWindowTitle(@window, title)
GLFW.SetWindowSize(@window, width, height)

then = time()
t = 0.0

try
   while !GLFW.WindowShouldClose(@window)
      GLFW.PollEvents()

      now = time()
      dt = now - then
      global then = now
      global t += dt

      bg = parse(RGBA, background)

      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
      glClearColor(bg.r, bg.g, bg.b, bg.alpha)
      glViewport(0, 0, (@size)...)

      onupdate(dt)

      GLFW.SwapBuffers(@window)
   end
catch err
   @error "Error on update callback" exception=err
   Base.show_backtrace(stderr, catch_backtrace())
finally
   if @isdefined(onclose)
      onclose()
   end
   GLFW.DestroyWindow(@window)
end

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
using GuillotineModels.InstanceReader

# What do we want to test:
# 0) If two instances have the same number of pieces, L, and W.
# 1) If the two instances have the same plate types.
# 1.1) If they also share profit.
# 1.2) If they also share demand.
# 2) If reversing L and W (l and w) we then get the same as (1).

function b2s(b)
	return b ? "EQUAL" : "different"
end

macro tell_if_equal(a, b)
	v_a, v_b = esc(a), esc(b)
	s_a, s_b = string(a), string(b)
	return quote
		equal = $v_a == $v_b
		println($s_a * " == " * $s_b * ": " * b2s(equal))
		equal
	end
end

# NOTE: assume the pieces are sorted by their fields, in order.
function check_ptype_uniqueness(l, w, p, d)
	@assert isone(length(unique(length.([l, w, p, d]))))
	@assert length(l) > 0
	is_first = true
	llv, lwv, lpv, ldv = first.((l, w, p, d))
	for (lv, wv, pv, dv) in zip(l, w, p, d)
		if is_first
			is_first = false
			continue
		end
		if lv == llv && wv == lwv
			if pv == lpv
				if dv == ldv
					println("WARNING: two or more pieces have the same length, width," *
					" profit, and demand.")
				else
					println("WARNING: two or more pieces have the same length, width," *
					" and profit.")
				end
			else
				println("WARNING: two or more pieces have the same length and width.")
			end
			return
		end
		llv, lwv, lpv, ldv = lv, wv, pv, dv
	end
	return
end

# NOTE:
# We disregard:
# 1. The order of the pieces inside the instance.
# 2. The fact what one considers to be length the other considers to be width.
# We should also:
# 1. If the same piece type appears multiple times the lines should be combined.
function compare(path_1, path_2)
	N1, L1, W1, l1, w1, p1, q1 = read_from_file(path_1)
	N2, L2, W2, l2, w2, p2, q2 = read_from_file(path_2)
	# Sort the pieces in both instances.
	l1, w1, p1, q1 = map(i -> getindex.(sort(collect(zip(l1, w1, p1, q1))), i), 1:4)
	l2, w2, p2, q2 = map(i -> getindex.(sort(collect(zip(l2, w2, p2, q2))), i), 1:4)

	println("Comparing $(basename(path_1)) with $(basename(path_2))")
	check_ptype_uniqueness(l1, w1, p1, q1)
	check_ptype_uniqueness(l2, w2, p2, q2)
	if @tell_if_equal N1 N2
		if @tell_if_equal (L1, W1) (L2, W2)
			if @tell_if_equal (l1, w1) (l2, w2)
				@tell_if_equal p1 p2
				@tell_if_equal q1 q2
			end
		elseif @tell_if_equal (L1, W1) (W2, L2)
			if @tell_if_equal (l1, w1) (w2, l2)
				@tell_if_equal p1 p2
				@tell_if_equal q1 q2
			end
		end
	end
end

function compare_pairwise(folder_path)
	list = readdir(folder_path; join = true, sort = true)
	n = length(list)
	for i = 1:n, j = (i+1):n
		compare(list[i], list[j])
	end
	return
end

#compare(ARGS[1], ARGS[2])
compare_pairwise(ARGS[1])


text:
module KalmanFilter

using LinearAlgebra, Distributions

include("src/types.jl")
include("src/filter.jl")
include("src/smoother.jl")
include("src/simulate.jl")
include("example/local_level.jl")
include("example/kca.jl")
include("example/multivariate_local_level.jl")

# Export types
export Model, Output

# Export functions
export kalmanfilter, kalmansmoother

end # module


text:
@testset "Retirement" begin

@testset "Constructor" begin
    retirement = Retirement()
    @test all( [retirement.maxCareerLength == 0, retirement.retirementAge == 0,
        retirement.freq == 1, retirement.offset == 0, retirement.isEither] )
end  # @testset "Constructor"

retirement = Retirement()

@testset "function setRetirementCareerLength!" begin
    setRetirementCareerLength!( retirement, 120 )
    @test retirement.maxCareerLength == 120
    @test !setRetirementCareerLength!( retirement, -60 )
end  # @testset "function setRetirementCareerLength!"

@testset "function setRetirementAge!" begin
    setRetirementAge!( retirement, 600 )
    @test retirement.retirementAge == 600
    @test !setRetirementAge!( retirement, -60 )
end  # @testset "function setRetirementAge!"

@testset "function setRetirementSchedule!" begin
    setRetirementSchedule!( retirement, 12, 18 )
    @test ( retirement.freq == 12.0 ) &&
        ( retirement.offset == 6.0 )
    setRetirementSchedule!( retirement, 6, -4 )
    @test ( retirement.freq == 6.0 ) &&
        ( retirement.offset == 2.0 )
    @test !setRetirementSchedule!( retirement, -24 )
end  # @testset "function setRetirementSchedule!"

@testset "function setRetirementIsEither!" begin
    setRetirementIsEither!( retirement, false )
    @test !retirement.isEither
end  # @testset "function setRetirementIsEither!"

end  # @testset "Retirement"


text:
### 2-D Recipes

## Scatter plots
# real vectors
"""
    scatter(x, y [, axes] [, curve] ; args...) -> Gaston.Figure

Create a scatter plot with vectors `x` and `y`. `axes` and `curve` specify the
axes and curve configurations, respectively.

    scatter(c, args...) -> Gaston.Figure

If `c` is a complex vector, Create a scatter plot of the real and imaginary
parts of `c`.
"""
scatter(x, y ; args...) = plot(x, y, w=:points ; args...)
scatter(x, y, a::Axes ; args...) = plot(x, y, nothing, a, w=:points ; args...)

# complex vectors
scatter(y::ComplexCoord ; args...) = scatter(real(y), imag(y) ; args...)
scatter(y::ComplexCoord, a::Axes ; args...) = scatter(real(y), imag(y), a ; args...)

# add curves
"""
    scatter!(...) -> Gaston.Figure

Add a new scatter plot to an existing plot.
"""
scatter!(x, y ; args...) = plot!(x, y, ps=:points ; args...)
scatter!(y::ComplexCoord ; args...) = scatter!(real(y), imag(y); args...)

## Stem plots
"""
    stem(x, y, onlyimpulses=false, [,axes] [, curve], args...) -> Gaston.Figure

Create a stem plot using `x` and `y`. If `onlyimpulses` is true, then the
conventional circles indicating each sample are not drawn.

    stem(x, f::Function, onlyimpulses=false, [,axes] [, curve], args...) -> Gaston.Figure

Use `y = f.(x)`
"""
function stem(x, y, a::Axes = Axes() ; onlyimpulses=false, args...)
    p = plot(x, y, nothing, a ; w=:impulses, lc=:blue, lw=1.25, args...)
    onlyimpulses || (p = plot!(x, y ; w=:points, lc=:blue, pt="ecircle", ps=1.5, args...))
    return p
end
stem(y ; args...) = stem(1:length(y), y ; args...)
stem(y, a::Axes ; args...) = stem(1:length(y), y, a ; args...)
stem(x, f::Function ; args...) = stem( x, f.(x) ; args...)
stem(x, f::Function, a::Axes ; args...) = stem( x, f.(x), a ; args...)

## Bar plots
"""
    bar([x,] y, [axes,] args...) -> Gaston.Figure

Generate a bar plot with axes configuration `axes`.
"""
function bar(x, y, a::Axes = Axes() ; args...)
    plot(x, y, nothing,
         merge(Axes(boxwidth="0.8 relative", style="fill solid 0.5"), a),
         w=:boxes ; args...)
end
bar(y ; args...) = bar(1:length(y), y ; args...)
bar(y, a::Axes ; args...) = bar(1:length(y), y, a ; args...)

## Histograms
"""
    histogram(data, [axes,] bins=10, norm=1.0, args...) -> Gaston.Figure

Plot a histogram of `data`. `bins` specifies the number of bins (default 10),
and the histogram area is normalized to `norm` (default 1.0).
"""
function histogram(data, a::Axes = Axes() ; bins::Int=10, norm::Real=1.0, args...)
    # validation
    bins < 1 && throw(DomainError(bins, "at least one bin is required"))
    norm < 0 && throw(DomainError(norm, "norm must be a positive number."))
    # calculate data
    x, y = hist(data,bins)
    # make area under histogram equal to norm
    y = norm*y/(step(x)*sum(y))
    # plot
    bar(x, y, a; args...)
end

## Images
"""
    imagesc([x,] [y,] z, [axes,] args...) -> Gaston.Figure

Plot an image given by array `z`. If the array is a matrix, a grayscale image
is assumed. If the array is three-dimensional, an rgbimage is assumed, with
`z[1,:,:]` the red channel, `z[2,:,:]` the blue channel, and `z[3,:,:]` the
blue channel.
"""
function imagesc(x, y, z, a::Axes = Axes() ; args...)
    ps = "image"
    ndims(z) == 3 && (ps = "rgbimage")
    plot(x, y, z, a, w=ps ; args...)
end
# grayscale
function imagesc(z::Matrix, a::Axes = Axes() ; args...)
    imagesc(1:size(z)[2], 1:size(z)[1], z, a ; args...)
end
# rgb
function imagesc(z::AbstractArray{<:Real,3}, a::Axes = Axes() ; args...)
    imagesc(1:size(z)[3], 1:size(z)[2], z, a ; args...)
end

### 3-D recipes

## Surfaces
"""
    surf([x,] [y,] z, [axes,] args...)-> Gaston.Figure

Plot the surface specified by `z` with axes specification `axes`.

    surf(x, y f::Function, args...)-> Gaston.Figure

Plot the surface `f.(x,y)`.
"""
surf(x, y, z, a::Axes = Axes() ; args...) = plot(x, y, z, a, dims=3 ; args...)
surf!(x, y, z ; args...) = plot!(x, y, z, dims=3 ; args...)
# function
function surf(x, y, f::Function, a::Axes = Axes() ; args...)
    plot(x, y, meshgrid(x, y, f), a, dims=3 ; args...)
end

"""
    surf!(...)

Add a surface to an existing figure.
"""
surf!(x, y, f::Function ; args...) = plot!(x, y, meshgrid(x, y, f), dims=3 ; args...)
# matrix
function surf(z::Matrix, a::Axes = Axes() ; args...)
    surf(1:size(z)[2], 1:size(z)[1], z, a, dims=3 ; args...)
end
surf!(z::Matrix ; args...) = surf!(1:size(z)[2], 1:size(z)[1], z, dims=3 ; args...)

## 3-D scatter
"""
    scatter3(x, y, z, [axes,] [curve,], args...) -> Gaston.Figure

Create a 3-D scatter plot with vectors `x`, `y` and `z`. `axes` and `curve`
specify the axes and curve configurations, respectively.

    scatter3(x, f1::Function, f2::Function, args...) -> Gaston.Figure

Create the scatter plot of `x`, `f1.(x)`, and `f2.(x)`.
"""
function scatter3(x, y, z, a::Axes = Axes() ; args...)
    surf(x, y, z, a, w=:points ; args...)
end
scatter3!(x, y, z ; args...) = surf!(x, y, z, w=:points ; args...)
# functions
function scatter3(x, f1::Function, f2::Function, a::Axes = Axes() ; args...)
    scatter3(x, f1.(x), f2.(x), a, dims=3 ; args...)
end
""""
    scatter3!(...) -> Gaston.Figure

Add a new scatter plot to an existing plot.
"""
function scatter3!(x, f1::Function, f2::Function ; args...)
    scatter3!(x, f1.(x), f2.(x), dims=3 ; args...)
end

# 3-D contour plots
"""
    contour([x,] [y,] z::Matrix, labels = true, [axes,] args...) -> Gaston.Figure

Plot the contour lines given by the surface `z`. If `labels = false`, no labels
are included in the plot.

    contour(x, y, f::Function, args...) -> Gaston.Figure

Plot the contours of the surface `f.(x,y)`.
"""
function contour(x, y, z, a::Axes = Axes() ; labels=true, args...)
    conaxes = Axes(key = false,
                   view = "map",
                   contour = "base",
                   surface = false,
                   cntrlabel = "font '7'",
                   cntrparam = "levels 10"
                  )
    p = surf(x, y, z, merge(conaxes, a) ; args...)
    labels && (p = surf!(x, y, z ; plotstyle=:labels))
    return p
end
contour(x, y, f::Function ; args...) = contour(x, y, meshgrid(x, y, f) ; args...)
contour(z::Matrix ; args...) = contour(1:size(z)[2], 1:size(z)[1], z ; args...)

# 3-D heatmaps
"""
    heatmap(x, y, z::Matrix, [axes,] args...)

Plot the heatmap of the surface specified by `z`.
"""
function heatmap(x, y, z, a::Axes = Axes(); args...)
    surf(x, y, z, merge(Axes(view="map"), a), w=:pm3d ; args...)
end


text:
export chamfer_distance

"""
    chamfer_distance(A::PointCloud, B::PointCloud; w1=1.0, w2=1.0)

Computes the chamfer distance between PointCloud `A` and `B`.

`w1` and `w2` are optional arguments for specifying the weighted mean of between two
pcloud for computing metrics.
"""
chamfer_distance(A::PointCloud, B::PointCloud; w1::Number = 1.0, w2::Number = 1.0) =
    _chamfer_distance(A.points, B.points, Float32(w1), Float32(w2))

chamfer_distance(
    A::AbstractArray{<:Number},
    B::AbstractArray{<:Number};
    w1::Number = 1.0,
    w2::Number = 1.0,
) = _chamfer_distance(Float32.(A), Float32.(B), Float32(w1), Float32(w2))

chamfer_distance(
    A::AbstractArray{Float32},
    B::AbstractArray{Float32};
    w1::Number = 1.0,
    w2::Number = 1.0,
) = _chamfer_distance(A, B, Float32(w1), Float32(w2))

function _chamfer_distance(
    A::AbstractArray{Float32,2},
    B::AbstractArray{Float32,2},
    w1::Float32 = 1.0f0,
    w2::Float32 = 1.0f0,
)
    A = reshape(A, size(A)..., 1)
    B = reshape(B, size(B)..., 1)
    return _chamfer_distance(A, B, w1, w2)
end

function _chamfer_distance(
    A::AbstractArray{Float32,3},
    B::AbstractArray{Float32,3},
    w1::Float32 = 1.0f0,
    w2::Float32 = 1.0f0,
)
    nn_for_A, nn_for_B = @ignore _nearest_neighbors(A, B)

    dist_A_to_B = mean((A .- B[:, nn_for_A]) .^ 2) * 3.0f0
    dist_B_to_A = mean((B .- A[:, nn_for_B]) .^ 2) * 3.0f0

    distance = (w1 * dist_A_to_B) + (w2 * dist_B_to_A)
    return distance
end

function _nearest_neighbors(x::Array{Float32,3}, y::Array{Float32,3})
    nn_for_x = cat(
        [
            CartesianIndex.(reduce(vcat, knn(KDTree(y[:, :, i]), x[:, :, i], 1)[1]), i)
            for i = 1:size(x, 3)
        ]...,
        dims = 2,
    )
    nn_for_y = cat(
        [
            CartesianIndex.(reduce(vcat, knn(KDTree(x[:, :, i]), y[:, :, i], 1)[1]), i)
            for i = 1:size(x, 3)
        ]...,
        dims = 2,
    )
    return nn_for_x, nn_for_y
end

function _nearest_neighbors(x::CuArray{Float32,3}, y::CuArray{Float32,3})
    xx = sum(x .^ 2, dims = 1)
    yy = sum(y .^ 2, dims = 1)
    zz = Flux.batched_mul(permutedims(x, (2, 1, 3)), y)
    rx = reshape(xx, size(xx, 2), 1, :)
    ry = reshape(yy, 1, size(yy, 2), :)
    P = (rx .+ ry) .- (2 .* zz)

    nn_for_x = argmin(P; dims = 2) |> Array
    nn_for_y = argmin(P; dims = 1) |> Array

    nn_for_x = reshape(map(x -> CartesianIndex(x[2], x[3]), nn_for_x), :, size(x, 3))
    nn_for_y = reshape(map(y -> CartesianIndex(y[1], y[3]), nn_for_y), :, size(y, 3))
    return nn_for_x, nn_for_y
end


text:
# https://en.wikipedia.org/wiki/Tak_(function)
module BenchTarai

using BenchmarkTools

module SeqTarai

tarai(x, y, z) =
    if y < x
        tarai(tarai(x - 1, y, z), tarai(y - 1, z, x), tarai(z - 1, x, y))
    else
        y
    end

end  # module SeqTarai

module BaseTarai

tarai(x, y, z) =
    if y < x
        a = Threads.@spawn tarai(x - 1, y, z)
        b = Threads.@spawn tarai(y - 1, z, x)
        c = tarai(z - 1, x, y)
        tarai(fetch(a)::Int, fetch(b)::Int, c)
    else
        y
    end

end  # module BaseTarai

module TapirTarai

using Base.Experimental: Tapir

tarai(x, y, z) =
    if y < x
        Tapir.@output a b c
        Tapir.@sync begin
            Tapir.@spawn a = tarai(x - 1, y, z)
            Tapir.@spawn b = tarai(y - 1, z, x)
            c = tarai(z - 1, x, y)
        end
        tarai(a, b, c)
    else
        y
    end

end  # module TapirTarai

module WSTarai

using Base.Experimental: Tapir
using TapirSchedulers

tarai(x, y, z) =
    if y < x
        Tapir.@output a b c
        Tapir.@sync WorkStealingTaskGroup() begin
            Tapir.@spawn a = tarai(x - 1, y, z)
            Tapir.@spawn b = tarai(y - 1, z, x)
            c = tarai(z - 1, x, y)
        end
        tarai(a, b, c)
    else
        y
    end

end  # module WSTarai

function setup(xyz = [(3, 1, 10)])
    suite = BenchmarkGroup()
    for (x::Int, y::Int, z::Int) in xyz
        s0 = suite["x=$x, y=$y, z=$z"] = BenchmarkGroup()
        s0["seq"] = @benchmarkable SeqTarai.tarai($x, $y, $z)
        s0["base"] = @benchmarkable BaseTarai.tarai($x, $y, $z)
        s0["tapir"] = @benchmarkable TapirTarai.tarai($x, $y, $z)
        s0["ws"] = @benchmarkable WSTarai.tarai($x, $y, $z)
    end
    return suite
end

function clear() end

end  # module BenchTarai


text:
addedge!(srccont::C, dstcont, dst, op, conf) where C<:AbstractContainer =
  addedge!(root(srccont), dstcont, dst, op, conf)
deledge!(srccont::C, dst, op) where C<:AbstractContainer =
  deledge!(root(srccont), dst, op)
deledges!(srccont::C) where C<:AbstractContainer =
  deledges!(root(srccont))

nupdate!(cont::CPUContainer{C}) where C<:AbstractNode = nupdate!(transient(cont))
nupdate!(x::Any) = ()

sync!(cont::C) where C<:AbstractContainer = sync!(cont.root) # TODO: First pull down transient?
sync!(num::N) where N<:Number = ()
sync!(arr::A) where A<:AbstractArray = Mmap.sync!(arr)

mutate!(cont::CPUContainer{C}, args...) where C = mutate!(transient(cont), args...)

defaultvalue(obj, param1, params...) = defaultvalue(obj[param1], params...)
defaultvalue(obj, param::T) where T<:Val = defaultvalue(obj[param])
defaultvalue(x::T) where T<:Number = zero(T)

bounds(obj, param1, params...) = bounds(obj[param1], params...)
bounds(obj, param::T) where T<:Val = bounds(obj[param])
bounds(x::T) where T<:Number = (zero(T), one(T))

Base.size(cont::CPUContainer{C}) where C<:AbstractNode =
  size(root(cont))

Base.getindex(cont::CPUContainer{C}, idx...) where C<:AbstractNode =
  getindex(root(cont), idx...)


text:
#   This file is part of UnitSystems.jl
#   It is licensed under the MIT license
#   UnitSystems Copyright (C) 2022 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

export Constant, constant

# constant

struct Constant{D} <: Real
    @pure Constant{D}() where D = new{D}()
end

@pure param(::Constant{D}) where D = D
Base.float(x::Constant) = float(constant(x))
Base.convert(::Type{Float64},c::Constant) = float(c)
logdb(x::Constant{D}) where D = Constant{logdb(D)}()
Base.:+(a::Constant,b::Constant) = Constant(constant(a)+constant(b))
Base.:-(a::Constant,b::Constant) = Constant(constant(a)-constant(b))
Base.:+(a::Constant{D},::Constant{D}) where D = 𝟐*a
Base.:-(a::Constant{D},::Constant{D}) where D = Constant(0)
Base.:+(a::Number,b::Constant) = a+constant(b)
Base.:+(a::Constant,b::Number) = constant(a)+b
Base.:-(a::Number,b::Constant) = a-constant(b)
Base.:-(a::Constant,b::Number) = constant(a)-b
Base.:*(a::Real,b::Constant) = a*constant(b)
Base.:*(a::Constant,b::Real) = constant(a)*b
Base.:*(a::Constant{A},b::Constant{B}) where {A,B} = Constant{A*B}()
Base.:/(a::Constant{A},b::Constant{B}) where {A,B} = Constant{A/B}()
Base.:/(a::Number,b::Constant) = a*inv(b)
Base.:/(a::Constant,b::Number) = a*inv(b)
Base.inv(a::Constant{D}) where D = Constant{inv(D)}()
Base.sqrt(a::Constant{D}) where D = Constant{sqrt(D)}()
Base.cbrt(a::Constant{D}) where D = Constant{cbrt(D)}()
Base.log(x::Constant{D}) where D = Constant{log(D)}()
Base.log2(x::Constant{D}) where D = Constant{log2(D)}()
Base.log10(x::Constant{D}) where D = Constant{log10(D)}()
Base.log(b::Number,x::Constant{D}) where D = Constant{log(b,D)}()
Base.exp(x::Constant{D}) where D = Constant{exp(D)}()
Base.exp2(x::Constant{D}) where D = Constant{exp2(D)}()
Base.exp10(x::Constant{D}) where D = Constant{exp10(D)}()
Base.:^(a::Number,b::Constant{D}) where D = Constant{a^D}()
Base.:^(a::Constant{D},b::Number) where D = Constant{D^b}()
Base.:^(a::Constant{D},b::Integer) where D = Constant{D^b}()
Base.:^(a::Constant{D},b::Rational{Int}) where D = Constant{D^b}()

Base.isone(x::Constant{D}) where D = isone(D)
Base.iszero(x::Constant{D}) where D = iszero(D)

Base.:(==)(a::Real,b::Constant) = a == constant(b)
Base.:(==)(a::Constant,b::Real) = constant(a) == b
Base.:(==)(a::Constant,b::Constant) = constant(a) == constant(b)
Base.isapprox(a::Real,b::Constant) = isapprox(a,constant(b))
Base.isapprox(a::Constant,b::Real) = isapprox(constant(a),b)
Base.isapprox(a::Constant,b::Constant) = isapprox(constant(a),constant(b))


text:
export ft_purge
export FT_PURGE_RX, FT_PURGE_TX

#Purge RX and TX Buffers (see ft_purge)
const FT_PURGE_RX = 1
const FT_PURGE_TX = 2

function ft_purge(ft_handle::Culong, mask::Integer)
  ft_status = ccall((:FT_Purge, d2xx),
                     Cuint,
                     (Culong, Cuint),
                     ft_handle, mask)
  checkstatus(ft_status)
  return nothing
end


text:
type DiscretizeUniformCount <: DiscretizationAlgorithm
    nbins::Int
end

function binedges{N<:AbstractFloat}(alg::DiscretizeUniformCount, data::AbstractArray{N})

    nbins = alg.nbins

    n = length(data)
    n ≥ nbins || error("too many bins requested")

    p = sortperm(data)
    counts_per_bin, remainder = div(n,nbins), rem(n,nbins)
    retval = Array(N, nbins+1)
    retval[1] = data[p[1]]
    retval[end] = data[p[end]]

    ind = 0
    for i in 2 : nbins
        counts = counts_per_bin + (remainder > 0.0 ? 1 : 0)
        remainder -= 1.0
        ind += counts
        retval[i] = (data[p[ind]] + data[p[ind+1]])/2
        retval[i-1] != retval[i] || error("binedges non-unique") # TODO(tim): should make the algorithm handle this
    end

    retval
end
function binedges{N<:Integer}(alg::DiscretizeUniformCount, data::AbstractArray{N})

    nbins = alg.nbins

    n = length(data)
    n ≥ nbins || error("too many bins requested")

    p = sortperm(data)
    counts_per_bin, remainder = div(n,nbins), rem(n,nbins)
    retval = Array(N, nbins+1)
    retval[1] = data[p[1]]
    retval[end] = data[p[end]]

    ind = 0
    for i in 2 : nbins
        counts = counts_per_bin + (remainder > 0.0 ? 1 : 0)
        remainder -= 1.0
        ind += counts
        retval[i] = ceil(Int, (data[p[ind]] + data[p[ind+1]])/2) # value V will be placed in bin B if V ∈ [Bₗ Bᵣ)
        retval[i-1] != retval[i] || error("binedges non-unique") # TODO(tim): should make the algorithm handle this
    end

    retval
end


text:
function get_corrected(obj::KalmanFilter)
    return jcall(obj, "getCorrected", ProcessEstimate, ())
end

function get_predicted(obj::KalmanFilter)
    return jcall(obj, "getPredicted", ProcessEstimate, ())
end


text:
module ComplexMixtures

using Printf
using Parameters
using ProgressMeter
using Statistics
using FortranFiles
using PDBTools
using StructTypes
using JSON3
using ThreadPools
using StaticArrays
using DocStringExtensions
import Random

export Selection, Trajectory, mddf, save, load, write, Options, Result
export contrib, merge
export overview, gr, grid3D
export VMDSelect

# Input and Output data structures
include("./FileOperations.jl")
include("./OutputFiles.jl")
include("./Density.jl")
include("./Volume.jl")
include("./Options.jl")
include("./Units.jl")

# Function to rigid-body move molecules
include("./random.jl")
include("./centerofcoordinates.jl")
include("./eulermat.jl")
include("./random_move.jl")
include("./move.jl")
include("./center_to_origin.jl")
include("./writexyz.jl")

# Structures and functions to deal with the solute and solvent selections
include("./Selection.jl")
include("./itype.jl")

# Select solute and solvent using VMD
include("./VMDselect.jl")

# Structures and functions to read different types of trajectories
include("./Trajectory.jl")
include("./trajectory_formats/ChemFiles.jl")
include("./trajectory_formats/NamdDCD.jl")
include("./trajectory_formats/PDBTraj.jl")

# Structures used to store results
include("./isautocorrelation.jl")
include("./SolSummary.jl")
include("./Result.jl")
include("./Samples.jl")

# Function to print some data about the run
include("./title.jl")

# Structure used for periodic boundary conditions
include("./wrap.jl")
include("./Box.jl")

# Functions to compute distances
include("./distance.jl")
include("./minimumdistance.jl")

# Functions to construct histograms
include("./sphericalshellvolume.jl")
include("./shellradius.jl")
include("./sphereradiusfromshellvolume.jl")
include("./setbin.jl")
include("./viewmol.jl")
include("./inbulk.jl")

# Structures to report results
include("./finalresults.jl")
include("./merge.jl")
include("./save.jl")
include("./load.jl")
include("./write.jl")
include("./which_types.jl")
include("./contrib.jl")

# Tools
include("./tools/gr.jl")
include("./tools/overview.jl")
include("./tools/grid3D.jl")
include("./isapprox.jl")

# Structures and functions for the linked cell method
include("./LinkedCells.jl")
include("./DminMol.jl")
include("./CutoffDistances_Struct.jl")
include("./partialsort_cutoff.jl")
include("./updatecounters.jl")
include("./update_counters_frame.jl")
include("./icell1D.jl")
include("./wrap_cell.jl")
include("./icell3D.jl")
include("./initcells.jl")
include("./increase_size.jl")
include("./cutoffdcell.jl")
include("./cutoffdistances.jl")
include("./mddf_linkedcells.jl")
include("./cutoffdcell_self.jl")
include("./cutoffdistances_self.jl")

# for the linked-cell-parallel version
include("./FrameData.jl")
include("./sum.jl")
include("./mddf_frame.jl")
include("./mddf_frame_self.jl")
include("./mddf_linkedcells_parallel.jl")

# Main function
include("./mddf.jl")

end


text:
export  alphanumeric_names

# ----------
# types

# ----------
# functions
function learn{S<:DSL_LearningAlgorithm}( mat::Matrix{Int}, alg::Type{S} )
	dset = matrix_to_dataset(mat)
	net  = Network()
	worked = learn!( net, dset, alg )
	(net, worked)
end
function learn{S<:DSL_LearnParams}( mat::Matrix{Int}, params::S )
	dset = matrix_to_dataset(mat)
	net  = Network()
	worked = learn!( net, dset, params )
	(net, worked)
end
function learn!{S<:DSL_LearningAlgorithm}( net::Network, mat::Matrix{Int}, alg::Type{S} )
	dset = matrix_to_dataset(mat)
	learn!( net, dset, alg )
end
function learn!{S<:DSL_LearnParams}( net::Network, mat::Matrix{Int}, params::S )
	dset = matrix_to_dataset(mat)
	learn!( net, dset, params )
end

function matrix_to_dataset{R <: Real, S<:AbstractString}( mat::Matrix{R}, names::Vector{S}=alphanumeric_names(size(mat,2)) )
	# converts a Matrix{Int} to a Dataset by writing it to a tempfile and loading it back in
	# each column in mat is one variable's data

	ftempname = tempname()
	writetable(ftempname, mat, names)
	dset = Dataset()
	if !read_file(dset, ftempname)
		error("Could not read data file $ftempname")
	end
	rm(ftempname)
	dset
end

function writetable{R <: Real, S<:AbstractString}(filename::AbstractString, mat::Matrix{R}, names::Vector{S}; separator::Char='\t')
	n,m = size(mat)
	f = open(filename, "w")

	for name in names[1:end-1]
		@printf(f, "%s%c", name, separator)
	end
	@printf(f, "%s\n", names[end])

	if isa(R, Integer)
		for i in 1 : n
			for j in 1 : m-1
				@printf(f, "%d%c", mat[i,j], separator)
			end
			@printf(f, "%d\n", mat[i,m])
		end
	else
		for i in 1 : n
			for j in 1 : m-1
				@printf(f, "%f%c", mat[i,j], separator)
			end
			@printf(f, "%f\n", mat[i,m])
		end
	end

	close(f)
end

function increment_alphanumeric_string(n,i=1)
	if i > length(n)
		n = "a"^i
	elseif n[i] != 'z'
		n = n[1:i-1]*string(n[i]+1)*n[i+1:end]
	else
		n = increment_alphanumeric_string(n[1:i-1]*"a"*n[i+1:end], i+1)
	end
	n
end
function alphanumeric_names(n::Int)
	# generate the first n lexographic names
	# a,b,c,...,z,aa,ab,...
	@assert(n > 0)
	names = Array(AbstractString, n)
	names[1] = "a"
	for i in 2 : n
		names[i] = increment_alphanumeric_string(names[i-1])
	end
	names
end