text:
using Flux
using Random
using Test
using Flux: @functor

"""
    MCLayer(layer, dropout)
A generic Monte Carlo dropout layer. Takes in any "traditional" flux
layer and a function that implements dropout. Performs the usual layer
forward pass and then passes the acitvations through the given dropout function.
Let x be the input array, the forward pass then becomes -- dropout(layer(x))

# Fields
- `layer`: Traditional Flux layer
- `dropout`: A function that implements dropout
"""
struct MCLayer{L,F}
    layer::L
    dropout::F
end

@functor MCLayer

"""
    MCDense(in::Int, out::Int, dropout_rate, σ=identity; bias=true, init=glorot_uniform)
    MCDense(layer, dropout_rate)

Creates a traditional dense layer with MC dropout functionality.
MC Dropout simply means that dropout is activated in both train and test times
Reference - Dropout as a bayesian approximation - https://arxiv.org/abs/1506.02142

The traditional dense layer is a field in the struct MCDense, so all the
arguments required for the dense layer can be provided, or the layer can
be provided too. The forward pass is the affine transformation of the dense
layer followed by dropout applied on the resulting activations.

    y = dropout(σ.(W * x .+ bias), dropout_rate)

# Fields
- `layer`: A traditional dense layer
- `dropout`: A function that implements dropout

# Arguments
- `in::Integer`: Input dimension of features
- `out::Integer`: Output dimension of features
- `dropout_rate::AbstractFloat`: Dropout rate
- `σ::F=identity`: Activation function, defaults to identity
- `init=glorot_normal`: Initialization function, defaults to glorot_normal
"""
function MCDense(in::Integer, out::Integer, dropout_rate, σ = identity, kwargs...)
    layer = Flux.Dense(in, out, σ; kwargs...)
    dropout = (x; k...) -> Flux.dropout(x, dropout_rate; k...)
    return MCLayer(layer, dropout)
end

"""
    MCConv(filter, in => out, σ = identity;
            stride = 1, pad = 0, dilation = 1, groups = 1, [bias, weight, init])
    MCConv(layer, dropout_rate)

Creates a traditional Conv layer with MC dropout functionality.
MC Dropout simply means that dropout is activated in both train and test times

Reference - Dropout as a bayesian approximation - https://arxiv.org/abs/1506.02142

The traditional conv layer is a field in the struct MCConv, so all the
arguments required for the conv layer can be provided, or the layer can
be provided too. The forward pass is the conv operation of the conv
layer followed by dropout applied on the resulting activations.

    y = dropout(Conv(x), dropout_rate)

# Fields
- `layer`: A traditional conv layer
- `dropout_rate::AbstractFloat`: Dropout rate

# Arguments
- `filter::NTuple{N,Integer}`: Kernel dimensions, eg, (5, 5)
- `ch::Pair{<:Integer,<:Integer}`: Input channels => output channels
- `dropout_rate::AbstractFloat`: Dropout rate
- `σ::F=identity`: Activation function, defaults to identity
- `init=glorot_normal`: Initialization function, defaults to glorot_normal
"""
function MCConv(
    k::NTuple{N,Integer},
    ch::Pair{<:Integer,<:Integer},
    dropout_rate,
    σ = identity;
    kwargs...,
) where {N}
    layer = Flux.Conv(k, ch, σ; kwargs...)
    dropout = (x; k...) -> Flux.dropout(x, dropout_rate; k...)
    return MCLayer(layer, dropout)
end

function MCConv(
    w::AbstractArray{T,N},
    bias,
    dropout_rate,
    σ = identity,
    kwargs...,
) where {T,N}
    layer = Flux.Conv(w, bias, σ, kwargs...)
    dropout = (x; k...) -> Flux.dropout(x, dropout_rate; k...)
    return MCLayer(layer, dropout)
end

"""
The forward pass of a MC layer: Passes the input through the
usual layer first and then through a dropout layer.

# Arguments
- `x`: Input tensors
- `dropout=true`: Toggle to control dropout, it's preferred to keep
dropout always on, but just in case if it's needed.
"""
function (mc::MCLayer)(x; dropout = true)
    # Layer forward pass
    # Dropout on activations
    output = mc.dropout(mc.layer(x); active = dropout)
    return output
end


text:
# ENV_name = "gwdg"
# ENV_root = "output_gwdg"
# println("Environment: GWDG")
#
# function jobengine_head(dir::String, Ncores::Integer, gpu::Bool, number::Int64, hours::Int64=48)
#   return """#!/bin/sh
#   #BSUB -L /bin/bash
#   #BSUB -q $(gpu ? "gpu" : "mpi")
#   #BSUB -W $(hours):00
#   #BSUB -e /usr/users/bardenn/reconstruction/$dir/terminal-err
#   #BSUB -o /usr/users/bardenn/reconstruction/$dir/terminal-out
#   #BSUB -n $Ncores
#   #BSUB -R span[hosts=1]
#   #BSUB -R "rusage[ngpus_shared=$Ncores]"
#   #BSUB -a openmp
#   #BSUB -cwd /usr/users/bardenn/reconstruction/$dir
#   #BSUB -J $(jobname(dir,number))
#   #$(number > 1 ? "BSUB -w \"done($(jobname(dir,number-1)))\"" : "" )
#   """
# end
#
# """Launches the script via qsub, in case dir doesnt exist, create it"""
# function launch_job(dir::String, Ncores::Integer, gpu::Bool, julia_script::String, successive_jobs::Integer=1, hours::Int64=48)
#   run(`mkdir -p $dir`)
#
#   #Let's queue some jobs, hold_jid takes care of chaining
#   for n = 1:successive_jobs
#     jobfile = "job$n.sh"
#     head = jobengine_head(dir, Ncores, gpu, n, hours)
#     open("$dir/$jobfile", "w+") do file
#       write(file,"""$head
#
#       julia << EOF
#         addprocs($(Ncores > 1 ? Ncores : 0)) #Let's go into parallel mode if necessary
#         $julia_script
#         flush(STDOUT) #make sure output is flushed to terminal files
#       EOF
#       """)
#       flush(file)
#     end
#     submit_job("$dir/$jobfile")
#   end
# end
#
# """Submits a job to the queue"""
# function submit_job(path)
#   run(pipeline("$path", `bsub`))
# end

ext: jl


text:
# Example of DandelionWebSocket.jl:
# Send some text and binary frames to ws://echo.websocket.org,
# which echoes them back.

using DandelionWebSockets

# Explicitly import the callback functions that we're going to add more methods for.
import DandelionWebSockets: on_text, on_binary,
                            state_connecting, state_open, state_closing, state_closed

# A simple WebSocketHandler which sends a few messages, receives echoes back, and then sends a stop
# signal via a channel when it's done.
mutable struct EchoHandler <: WebSocketHandler
    connection::Union{WebSocketConnection, Nothing}
    stop_channel::Channel{Any}

    EchoHandler(chan::Channel{Any}) = new(nothing, chan)
end

# These are called when you get text/binary frames, respectively.
on_text(::EchoHandler, s::String)  = println("Received text: $s")
on_binary(::EchoHandler, data::AbstractVector{UInt8}) = println("Received data: $(String(data))")

# These are called when the WebSocket state changes.

state_closing(::EchoHandler)    = println("State: CLOSING")
function state_connecting(e::EchoHandler, c::WebSocketConnection)
    println("State: CONNECTING")
    e.connection = c
end

# Called when the connection is open, and ready to send/receive messages.
function state_open(handler::EchoHandler)
    println("State: OPEN")

    # Send some text frames, and a binary frame.
    @async begin
        texts = ["Hello", "world", "!"]

        for text in texts
            println("Sending  text: $text")
            send_text(handler.connection, text)
            sleep(0.5)
        end

        send_binary(handler.connection, b"Hello, binary!")
        sleep(0.5)

        # Signal the script that we're done sending all messages.
        # The script will then close the connection.
        put!(stop_chan, true)
    end
end

function state_closed(::EchoHandler)
    println("State: CLOSED")

    # Signal the script that the connection is closed.
    put!(stop_chan, true)
end

stop_chan = Channel{Any}(3)

# Create a WSClient, which we can use to connect and send frames.
client = WSClient()

handler = EchoHandler(stop_chan)

if length(ARGS) == 0
    uri = "ws://echo.websocket.org"
elseif length(ARGS) == 1
    uri = ARGS[1]
else
    println("Expect zero or one arguments!")
    exit(1)
end

println("Connecting to $uri... ")

wsconnect(client, uri, handler)

println("Connected.")

# The first message on `stop_chan` indicates that all messages have been sent, and we should
# close the connection.
take!(stop_chan)

stop(handler.connection)

# The second message on `stop_chan` indicates that the connection is closed, so we can exit.
take!(stop_chan)


text:
include("status.jl")

const ACTION_FUNCTIONS = Dict(

)

struct Action
    name::String
    prequisits::Dict{Status, Bool}
end

function Action_Activate(action::Action, player::Player)
    ACTION_FUNCTIONS[action.name](player)
end

max_stars_repo_name: Flobii/discord-bot

max_issues_repo_name: Flobii/discord-bot

max_forks_repo_name: Flobii/discord-bot


text:
"""
Zemax paraxial surface
"""
struct Paraxial{T} <: AbstractSurface{T}
    focal_length::T
end

function Paraxial(focal_length, aperture, n, t, id = nothing)
    OpticalSurface(
        Paraxial(focal_length),
        aperture,
        n,
        t,
        id,
    )
end

function sag(radius, s::Paraxial)
    return (zero(x))
end

function sag(x, y, s::Paraxial)
    return (zero(x))
end

function transfer_to_intersection!(ray, t, s::Paraxial)
    return (transfer_to_intersection!(ray, t, Sphere(zero(t))))
end

function refract!(ray, m0, s::Paraxial, m1, λ)
    l = ray.cx
    m = ray.cy
    n = ray.cz
    ux = l / n
    uy = m / n

    n0 = refractive_index(m0, λ)
    n1 = refractive_index(m1, λ)
    phi = one(s.focal_length) / s.focal_length

    ux_prime = (n0 * ux - ray.x * phi) / n1
    uy_prime = (n0 * uy - ray.y * phi) / n1

    n_prime = sqrt(one(ray.x) / (one(ray.x) + ux_prime^2 + uy_prime^2))
    l_prime = ux_prime * n_prime
    m_prime = uy_prime * n_prime

    ray.cx = l_prime
    ray.cy = m_prime
    ray.cz = n_prime
    return ray
end


text:
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::OneToRange) = r
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::AbstractStepRange) = r
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::AbstractStepRangeLen) = r
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::AbstractLinRange) = r

Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepSRange) = srange(-first(r), step=-step(r), length=length(r))
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepMRange) = mrange(-first(r), step=-step(r), length=length(r))
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepSRangeLen) = StepSRangeLen(-r.ref, -r.step, length(r), r.offset)
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepMRangeLen) = StepMRangeLen(-r.ref, -r.step, length(r), r.offset)
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::LinSRange) = LinSRange(-r.start, -r.stop, length(r))
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::LinMRange) = LinMRange(-r.start, -r.stop, length(r))

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Real, r::UnitSRange)
    if is_static(x)
        return srange(x + first(r), length=length(r))
    else
        return mrange(x + first(r), length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Real, r::UnitMRange)
    return mrange(x + first(r), length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::UnitSRange, x::Real)
    if is_static(x)
        return srange(first(r) + x, length=length(r))
    else
        return mrange(first(r) + x, length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::UnitMRange, x::Real)
    return mrange(x + first(r), length=length(r))
end
# For #18336 we need to prevent promotion of the step type: TODO
#Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::AbstractRange, x::Number) = range(first(r) + x, step=step(r), length=length(r))
#Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::AbstractRange) = range(x + first(r), step=step(r), length=length(r))

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::StepSRangeLen{T}, x::Number) where T
    if is_static(x)
        return StepSRangeLen{typeof(T(r.ref)+x)}(r.ref+x, r.step, length(r), r.offset)
    else
        return StepMRangeLen{typeof(T(r.ref)+x)}(r.ref+x, r.step, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::StepMRangeLen{T}, x::Number) where T
    return StepMRangeLen{typeof(T(r.ref)+x)}(r.ref+x, r.step, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::StepSRangeLen{T}) where T
    if is_static(x)
        return StepSRangeLen{typeof(x+T(r.ref))}(x + r.ref, r.step, length(r), r.offset)
    else
        return StepMRangeLen{typeof(x+T(r.ref))}(x + r.ref, r.step, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::StepMRangeLen{T}) where T
    return StepMRangeLen{typeof(x+T(r.ref))}(x + r.ref, r.step, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::LinSRange, x::Number)
    if is_static(x)
        return LinSRange(r.start + x, r.stop + x, length(r))
    else
        return LinMRange(r.start + x, r.stop + x, length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::LinMRange, x::Number)
    return LinMRange(r.start + x, r.stop + x, length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::LinSRange)
    if is_static(x)
        return LinSRange(x + r.start, x + r.stop, length(r))
    else
        return LinMRange(x + r.start, x + r.stop, length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), x::Number, r::LinMRange)
    return LinMRange(x + r.start, x + r.stop, length(r))
end

## TODO
#Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r1::AbstractRange, r2::AbstractRange) = r1 + r2
##

###
### :(-)
###
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::Union{UnitSRange,OneToSRange}, x::Number)
    if is_static(x)
        return range(first(r)-x, length=length(r))
    else
        return range(first(r)-x, length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::Union{UnitMRange,OneToMRange}, x::Number)
    return range(first(r)-x, length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepSRange, x::Number)
    if is_static(x)
        return mrange(first(r)-x, step=step(r), length=length(r))
    else
        return mrange(first(r)-x, step=step(r), length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepMRange, x::Number)
    return mrange(first(r)-x, step=step(r), length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::StepSRange)
    if is_static(x)
        return srange(x-first(r), step=-step(r), length=length(r))
    else
        return mrange(x-first(r), step=-step(r), length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::StepMRange)
    return mrange(x-first(r), step=-step(r), length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepSRangeLen{T}, x::Number) where T
    if is_static(x)
        return StepSRangeLen{typeof(T(r.ref)-x)}(r.ref - x, r.step, length(r), r.offset)
    else
        return StepMRangeLen{typeof(T(r.ref)-x)}(r.ref - x, r.step, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::StepMRangeLen{T}, x::Number) where T
    return StepMRangeLen{typeof(T(r.ref)-x)}(r.ref - x, r.step, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::StepSRangeLen{T}) where T
    if is_static(x)
        return StepSRangeLen{typeof(x-T(r.ref))}(x - r.ref, -r.step, length(r), r.offset)
    else
        return StepMRangeLen{typeof(x-T(r.ref))}(x - r.ref, -r.step, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::StepMRangeLen{T}) where T
    return StepMRangeLen{typeof(x-T(r.ref))}(x - r.ref, -r.step, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::LinSRange, x::Number)
    if is_static(x)
        return LinMRange(r.start - x, r.stop - x, length(r))
    else
        return LinMRange(r.start - x, r.stop - x, length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r::LinMRange, x::Number)
    return LinMRange(r.start - x, r.stop - x, length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::LinSRange)
    if is_static(x)
        return LinMRange(x - r.start, x - r.stop, length(r))
    else
        return LinMRange(x - r.start, x - r.stop, length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), x::Number, r::LinMRange)
    return LinMRange(x - r.start, x - r.stop, length(r))
end
# TODO
#function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(-), r1::AbstractRange, r2::AbstractRange)
#    r1 - r2
#end

###
### :(*)
###
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::Union{StepSRange,UnitSRange,OneToSRange})
    if is_static(x)
        return mrange(x*first(r), step=x*step(r), length=length(r))
    else
        return mrange(x*first(r), step=x*step(r), length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::Union{StepMRange,UnitMRange,OneToMRange})
    return mrange(x*first(r), step=x*step(r), length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::StepSRangeLen{T}) where {T}
    if is_static(x)
        return StepSRangeLen{typeof(x*T(r.ref))}(x*r.ref, x*r.step, length(r), r.offset)
    else
        return StepMRangeLen{typeof(x*T(r.ref))}(x*r.ref, x*r.step, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::StepMRangeLen{T}) where {T}
    return StepMRangeLen{typeof(x*T(r.ref))}(x*r.ref, x*r.step, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::LinSRange)
    if is_static(x)
        return LinSRange(x * r.start, x * r.stop, r.len)
    else
        return LinMRange(x * r.start, x * r.stop, r.len)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), x::Number, r::LinMRange)
    return LinMRange(x * r.start, x * r.stop, r.len)
end
# separate in case of noncommutative multiplication

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::Union{StepSRange,UnitSRange,OneToSRange}, x::Number)
    if is_static(x)
        return srange(first(r)*x, step=step(r)*x, length=length(r))
    else
        return mrange(first(r)*x, step=step(r)*x, length=length(r))
    end
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::Union{StepMRange,UnitMRange,OneToMRange}, x::Number)
    return mrange(first(r)*x, step=step(r)*x, length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::StepSRangeLen{T}, x::Number) where {T}
    if is_static(x)
        return StepSRangeLen{typeof(T(r.ref)*x)}(r.ref*x, r.step*x, length(r), r.offset)
    else
        return StepMRangeLen{typeof(T(r.ref)*x)}(r.ref*x, r.step*x, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::StepMRangeLen{T}, x::Number) where {T}
    return StepMRangeLen{typeof(T(r.ref)*x)}(r.ref*x, r.step*x, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::LinSRange, x::Number)
    if is_static(x)
        return LinMRange(r.start * x, r.stop * x, r.len)
    else
        return LinMRange(r.start * x, r.stop * x, r.len)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(*), r::LinMRange, x::Number)
    return LinMRange(r.start * x, r.stop * x, r.len)
end

###
### :(/)
###
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::Union{StepSRange,UnitSRange,OneToSRange}, x::Number)
    if is_static(x)
        return srange(first(r)/x, step=step(r)/x, length=length(r))
    else
        return mrange(first(r)/x, step=step(r)/x, length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::Union{StepMRange,UnitMRange,OneToMRange}, x::Number)
    return mrange(first(r)/x, step=step(r)/x, length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::StepSRangeLen{T}, x::Number) where {T}
    if is_static(x)
        return StepMRangeLen{typeof(T(r.ref)/x)}(r.ref/x, r.step/x, length(r), r.offset)
    else
        return StepMRangeLen{typeof(T(r.ref)/x)}(r.ref/x, r.step/x, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::StepMRangeLen{T}, x::Number) where {T}
    return StepMRangeLen{typeof(T(r.ref)/x)}(r.ref/x, r.step/x, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::LinSRange, x::Number)
    if is_static(x)
        return LinSRange(r.start / x, r.stop / x, r.len)
    else
        return LinMRange(r.start / x, r.stop / x, r.len)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(/), r::LinMRange, x::Number)
    return LinMRange(r.start / x, r.stop / x, r.len)
end

###
### :(\)
###

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::Union{StepSRange,UnitSRange,OneToSRange})
    if is_static(x)
        return srange(x\first(r), step=x\step(r), length=length(r))
    else
        return mrange(x\first(r), step=x\step(r), length=length(r))
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::Union{StepMRange,UnitMRange,OneToMRange})
    return mrange(x\first(r), step=x\step(r), length=length(r))
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::StepSRangeLen)
    if is_static(x)
        return StepSRangeLen(x\r.ref, x\r.step, length(r), r.offset)
    else
        return StepMRangeLen(x\r.ref, x\r.step, length(r), r.offset)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::StepMRangeLen)
    return StepMRangeLen(x\r.ref, x\r.step, length(r), r.offset)
end

function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::LinSRange)
    if is_static(x)
        return LinSRange(x \ r.start, x \ r.stop, r.len)
    else
        return LinMRange(x \ r.start, x \ r.stop, r.len)
    end
end
function Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(\), x::Number, r::LinMRange)
    return LinMRange(x \ r.start, x \ r.stop, r.len)
end

#= TODO
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::UnitRange) = big(r.start):big(last(r))
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::StepRange) = big(r.start):big(r.step):big(last(r))
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::StepRangeLen) = StepRangeLen(big(r.ref), big(r.step), length(r), r.offset)
Base.broadcasted(::DefaultArrayStyle{1}, ::typeof(big), r::LinRange) = LinRange(big(r.start), big(r.stop), length(r))
=#


text:
# ## Getting started

using MLJ, RDatasets

auto = dataset("ISLR", "Auto")
y, X = unpack(auto, ==(:MPG), col->true)
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=444);

# Note the use of `rng=` to seed the shuffling of indices so that the results are reproducible.

# ### Polynomial regression
#

@load LinearRegressor pkg=MLJLinearModels

# In this part we only build models with the Horsepower feature.
# Let's get a baseline:

lm = LinearRegressor()
mlm = machine(lm, select(X, :Horsepower), y)
fit!(mlm, rows=train)
rms(predict(mlm, rows=test), y[test])^2

# Note that we square the measure to  match the results obtained in the ISL labs where the mean squared error (here we use the `rms` which is the square root of that).
#
# We now want to build three polynomial models of degree 1, 2 and 3 respectively; we start by forming the corresponding feature matrix:

hp = X.Horsepower
Xhp = DataFrame(hp1=hp, hp2=hp.^2, hp3=hp.^3);

# Now we  can write a simple pipeline where the first step selects the features we want (and with it the degree of the polynomial) and the second is the linear regressor:

@pipeline LinMod(fs = FeatureSelector(features=[:hp1]),
                 lr = LinearRegressor());

# Then we can  instantiate and fit 3 models where we specify the features each time:

lrm = LinMod()
lr1 = machine(lrm, Xhp, y) # poly of degree 1 (line)
fit!(lr1, rows=train)

lrm.fs.features = [:hp1, :hp2] # poly of degree 2
lr2 = machine(lrm, Xhp, y)
fit!(lr2, rows=train)

lrm.fs.features = [:hp1, :hp2, :hp3] # poly of degree 3
lr3 = machine(lrm, Xhp, y)
fit!(lr3, rows=train)

# Let's check the performances on the test set

get_mse(lr) = rms(predict(lr, rows=test), y[test])^2

@show get_mse(lr1)
@show get_mse(lr2)
@show get_mse(lr3)

# ## K-Folds Cross Validation
#
# Let's crossvalidate over the degree of the  polynomial.
#
# **Note**: there's a  bit of gymnastics here because MLJ doesn't directly support a polynomial regression; see our tutorial on [tuning models](/pub/getting-started/model-tuning.html) for a gentler introduction to model tuning.

Xhp = DataFrame([hp.^i for i in 1:10])

cases = [[Symbol("x$j") for j in 1:i] for i in 1:10]
r = range(lrm, :(fs.features), values=cases)

tm = TunedModel(model=lrm, ranges=r, resampling=CV(nfolds=10), measure=rms)

# Now we're left with fitting the tuned model

mtm = machine(tm, Xhp, y)
fit!(mtm)
rep = report(mtm)
@show round.(rep.measurements.^2, digits=2)
@show argmin(rep.measurements)

#
# So the conclusion here is that the 5th order polynomial does quite well.
#
# In ISL they use a different seed so the results are a bit different but comparable.

# ## The Bootstrap
#
# _Bootstrapping is not currently supported in MLJ._

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text: using Test, Random, Statistics


text:
using AdaGram

type Context
    position::Int
    words::Array{String}
end

type Sense
    sense::String
    weight::Int
    has_utf::Bool
    contexts::Array{Context}
end

type Word
    plain_word::String
    word::String
    weight::Int
    number_of_senses::Int
    senses::Array{Sense}
end

const MIN_WORD_WEIGHT = 500
const MAX_NUMBER_OF_CONTEXTS = 500

function read_and_process_dataset(dataset_number, vm, dict, stopwords, CONTEXT_SEMIWIDTH, result_filename)
    aris = Float64[]
    f = open("datasets/dataset-$dataset_number", "r")
    while !eof(f)
        tokens = split(strip(readline(f)), '\t')
        plain_word = tokens[1]
        word = lowercase(tokens[2])
        word_weight = int32(tokens[3])
        number_of_senses = int32(tokens[4])
        senses = Array(Sense, number_of_senses)
        for i in 1:number_of_senses
            if eof(f) break end

            tokens = split(strip(readline(f)), '\t')
            has_utf = endof(tokens[1]) != length(tokens[1])
            sense = remove_utf_characters(lowercase(tokens[1]))
            sense_weight = int32(tokens[2])
            number_of_contexts = int32(tokens[3])
            contexts = Array(Context, number_of_contexts)

            for j in 1:number_of_contexts
                position = int32(readline(f)) + 1
                words = split(lowercase(strip(readline(f))), " ")
                if length(words) == 0
                    contexts[j] = Context(position, words)
                end
                left_context = reverse(words[1:min(position,length(words))])
                filter!(w -> !isempty(w) && w in keys(dict.word2id) && !(w in stopwords), left_context)
                left_bound = min(length(left_context), CONTEXT_SEMIWIDTH)
                if position < length(words)
                    right_context = words[(position+1):end]
                    filter!(w -> !isempty(w) && w in keys(dict.word2id) && !(w in stopwords), right_context)
                    right_bound = min(length(right_context), CONTEXT_SEMIWIDTH)
                    words = vcat(left_context[1:left_bound], right_context[1:right_bound])
                else
                    words = left_context[1:left_bound]
                end
                contexts[j] = Context(position, words)
            end

            senses[i] = Sense(sense, sense_weight, has_utf, contexts)
        end
        if endof(word) != length(word) || !(word in keys(dict.word2id))
            # skip non-ASCII words
            continue
        end
        a = process_one_word(Word(plain_word, word, word_weight, number_of_senses, senses), vm, dict)
        if a > -100
            push!(aris, a)
        end
    end
    @printf("%f\n", mean(aris))
    close(f)
end

function remove_utf_characters(s::String)
    if endof(s) == length(s) return s end
    r = Char[]
    for i in 1:endof(s)
        if nextind(s, i) == i + 1 && prevind(s, i + 1) == i
            push!(r, s[i])
        end
    end
    return join(r, "")
end

function get_number_of_senses(vm, dict, word)
    z = Array(Float64,T(vm))
    N = expected_pi!(z, vm, dict.word2id[word])
    return N
end

function get_senses_mask(vm, dict, word)
    z = Array(Float64,T(vm))
    expected_pi!(z, vm, dict.word2id[word])
    return z .> 0
end

function disambiguate_senses(vm, dict, word::Word, N)
    clustering_matrix = zeros(word.number_of_senses, T(vm))
    mask = get_senses_mask(vm, dict, word.word)
    for i in 1:word.number_of_senses
        for j in 1:min(length(word.senses[i].contexts), MAX_NUMBER_OF_CONTEXTS)
            p = disambiguate(vm, dict, word.word, word.senses[i].contexts[j].words)
            clustering_matrix[i, indmax(p)] += 1
        end
    end
    return clustering_matrix[:,mask]
end

function compute_ARI(x)
    a = sum(x, 2)
    b = sum(x, 1)
    index = sum(x.*(x-1)*0.5)
    a1 = sum(a.*(a-1)*0.5)
    b1 = sum(b.*(b-1)*0.5)
    n = sum(x)
    expected_index = a1*b1/(n*(n-1)*0.5)
    max_index = (a1+b1)/2
    return (index - expected_index) / (max_index - expected_index)
end

function process_one_word(word, vm, dict)
    if word.weight >= MIN_WORD_WEIGHT
        N = get_number_of_senses(vm, dict, word.word)
        if N > 1
            clustering_matrix = disambiguate_senses(vm, dict, word, N)
            number_of_all_contexts = sum(clustering_matrix)
            ARI = compute_ARI(clustering_matrix)

#            @printf(f_verbose, "%s\t%d\t%d\n", word.word, word.number_of_senses, N)
#            for i in 1:word.number_of_senses
#                @printf(f_verbose, "%s\n", word.senses[i].sense)
#            end
#            for i in 1:word.number_of_senses
#                for j in 1:N
#                    @printf(f_verbose, "%d ", clustering_matrix[i, j])
#                end
#                @printf(f_verbose, "\n")
#            end
        else
            # hack for faster processing
            ARI = 0
            number_of_all_contexts = 0
            for i in 1:word.number_of_senses
                number_of_all_contexts += length(word.senses[i].contexts)
            end
        end
#        @printf(f, "%s\t%d\t%d\t%d\t%f\n", word.word, N, word.number_of_senses, number_of_all_contexts, ARI)
#        flush(f)
#        flush(f_verbose)
        return ARI
    end
    return -100
end

using ArgParse

s = ArgParseSettings()

@add_arg_table s begin
  "model"
    help = "path to serialized model"
    arg_type = String
    required = true
  "result-filename"
    help = "where to save result"
    arg_type = String
    required = true
  "--window"
    help = "window size"
    arg_type = Int
    default = 5
end

args = parse_args(ARGS, s)

vm, dict = load_model(args["model"])
#println("read stopwords")
stopwords = readdlm("benchmark/english.txt")
#println("process dataset")
read_and_process_dataset(14, vm, dict, stopwords, args["window"], args["result-filename"])


text:
import ..VideoGraphicDataset
import ..anydifferent
import ..load

"""
    VIDEO_EXTS

Supported video extensions.
"""
const VIDEO_EXTS = begin
    raw = split.(VideoIO.FFMPEG.exe("-formats", collect=true)[5:end], " ")
    fields = filter.(!isempty, raw)
    map(f -> f[2], fields)
end

"""
    Frames <: VideoGraphicDataset{C}
    Frames(frames::Array{C,3}) where {C <: Colorant}
    Frames(frames::AbstractVector{Matrix{C}} where {C <: Colorant}

An image stack consists of a collection of 2-D images, running in time.
"""
mutable struct Frames{C <: Colorant} <: VideoGraphicDataset{C}
    frames::Array{C,3}
end

function Frames(frames::AbstractVector{Matrix{C}}) where {C <: Colorant}
    if anydifferent(size.(frames))
        error("frames must all have the same dimension")
    end
    Frames{C}(cat(frames...; dims=3))
end

function Base.show(io::IO, stack::Frames)
    println(io, typeof(stack))
    print(io, "  size: ", join(string.(size(stack)), "×"))
end

frames(stack::Frames) = stack.frames
Base.size(stack::Frames) = size(frames(stack))

"""
    load(::Type{Frames}, path; ext=nothing, ignoredots=true)

Load all frames from a video or image files from a directory into an
`Frames`.  If a directory is provided, only files with extention `ext` will
be loaded. If `isnothing(ext)`, all files will be load. Files whose name begins
with a `.` will be ignored if `ignoredots` is `true`.

When loading from image files, the filenames are sorted lexicographically. The
files must all have the same size, file extention and color space, e.g.
`RGB{N0f8}`, `Gray{N0f8}`, etc..

```julia
julia> load(Frames, "fromtiffs.avi")
Frames{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}}
  size: 20×20×5
```

```julia
julia> readdir()
6-element Array{String,1}:
 ".DS_Store"
 "000001.tiff"
 "000002.tiff"
 "000003.tiff"
 "000004.tiff"
 "000005.tiff"

julia> load(Frames, ".")
Frames{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}}
  size: 20×20×5

julia> load(Frames, "."; ignoredots=false)
ERROR: files must all have the same extension; got ["", ".tiff"]
[...]
```

You can filter files by extension.
```julia
julia> load(Frames, ".")
ERROR: files must all have the same extension; got [".png", ".tiff"]
[...]

julia> load(Frames, ".", ".tiff")
Frames{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}}
  size: 20×20×2
```
"""
function load(::Type{Frames}, path::AbstractString; ext::Union{Nothing,AbstractString}=nothing,
              ignoredots=true)
    if isdir(path)
        files = readdir(path)
        ignoredots && filter!(f -> f[1] != '.', files)
        !isnothing(ext) && filter!(f -> last(splitext(f)) == ext, files)
        if isempty(files)
            error("no files$(isnothing(ext) ? " " : "with extension \"$ext\" ")found in \"$path\"")
        end
        load(Frames, sort!(joinpath.(path, files)))
    elseif isfile(path)
        if lowercase(replace(last(splitext(path)), "." => "")) ∉ VIDEO_EXTS
            error("file extension is not a supported video extension, see CaimCore.VIDEO_EXTS")
        end
        Frames(loadvideo(path))
    else
        error("path is neither a directory nor a file")
    end
end

"""
    load(::Type{Frames}, files)

Load image files into an image stack in the order they are provided.

The files must all have the same size, file extention and color space, e.g.
`RGB{N0f8}`, `Gray{N0f8}`, etc..

```julia
julia> load(Frames, ["000001.tiff", "000002.tiff"])
Frames{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}}
  size: 20×20×2
```
"""
function load(::Type{Frames}, files::AbstractVector{<:AbstractString})
    if isempty(files)
        error("no files provided")
    end

    exts = unique(map(last ∘ splitext, files))
    if length(exts) != 1
        error("files must all have the same extension; got ", exts)
    end

    frames = loadstack(files)
    if anydifferent(typeof.(frames))
        error("files must all have the same color space")
    end

    try
        Frames(frames)
    catch e
        if isa(e, MethodError)
            error("files must all be images")
        else
            rethrow(e)
        end
    end
end

function loadstack(files::AbstractVector{<:AbstractString})
    try
        FileIO.load.(files)
    catch e
        if isa(e, FileIO.UnknownFormat)
            error("unrecognized file format found")
        else
            rethrow(e)
        end
    end
end

function loadvideo(path::AbstractString)
    video = try
        VideoIO.openvideo(path)
    catch
        error("file does not appear to be a properly formatted video")
    end
    @defer close(video)
    frames = [Array(read(video))]
    while !eof(video)
        push!(frames, Array(read(video)))
    end
    frames
end

Base.:(==)(a::Frames, b::Frames) = frames(a) == frames(b)


text:
function __init__udenglish()
    DEPNAME = "UD_English"
    TRAINFILE = "en_ewt-ud-train.conllu"
    DEVFILE   = "en_ewt-ud-dev.conllu"
    TESTFILE  = "en_ewt-ud-test.conllu"

    register(DataDep(
        DEPNAME,
        """
        Dataset: Universal Dependencies - English Dependency Treebank Universal Dependencies English Web Treebank
        Authors: Natalia Silveira and Timothy Dozat and
                 Marie-Catherine de Marneffe and Samuel
                 Bowman and Miriam Connor and John Bauer and
                 Christopher D. Manning
        Website: https://github.com/UniversalDependencies/UD_English-EWT

        The files are available for download at the github
        repository linked above. Note that using the data
        responsibly and respecting copyright remains your
        responsibility. Copyright and License is discussed in
        detail on the Website.
        """,
        "https://raw.githubusercontent.com/UniversalDependencies/UD_English-EWT/master/" .* [TRAINFILE, DEVFILE, TESTFILE],
        "6df6ee25ab3cd1cde3a09ab075dcc6b8c90d18648eef0809f400be4ad8bc81e2"
    ))
end

"""
    UD_English(; split=:train, dir=nothing)
    UD_English(split=; [dir])

A Gold Standard Universal Dependencies Corpus for
English, built over the source material of the
English Web Treebank LDC2012T13
(https://catalog.ldc.upenn.edu/LDC2012T13).

The corpus comprises 254,825 words and 16,621 sentences,
taken from five genres of web media: weblogs, newsgroups, emails, reviews, and Yahoo! answers.
See the LDC2012T13 documentation for more details on the sources of the sentences.
The trees were automatically converted into Stanford Dependencies and then hand-corrected to Universal Dependencies.
All the basic dependency annotations have been single-annotated, a limited portion of them have been double-annotated,
and subsequent correction has been done to improve consistency. Other aspects of the treebank, such as Universal POS,
features and enhanced dependencies, has mainly been done automatically, with very limited hand-correction.

Authors: Natalia Silveira and Timothy Dozat and
            Marie-Catherine de Marneffe and Samuel
            Bowman and Miriam Connor and John Bauer and
            Christopher D. Manning
Website: https://github.com/UniversalDependencies/UD_English-EWT
"""
struct UD_English <: UnsupervisedDataset
    metadata::Dict{String,Any}
    split::Symbol
    features::Vector{Vector{Vector{String}}}
end

UD_English(; split=:train, dir=nothing) = UD_English(split; dir)

function UD_English(split::Symbol; dir=nothing)
    DEPNAME = "UD_English"
    TRAINFILE = "en_ewt-ud-train.conllu"
    DEVFILE   = "en_ewt-ud-dev.conllu"
    TESTFILE  = "en_ewt-ud-test.conllu"

    @assert split ∈ [:train, :test, :dev]

    FILE = split == :train ? TRAINFILE :
           split == :test ? TESTFILE :
           split === :dev ? DEVFILE : error()

    path = datafile(DEPNAME, FILE, dir)

    doc = []
    sent = []
    lines = open(readlines, path)
    for line in lines
        line = chomp(line)
        if length(line) == 0
            length(sent) > 0 && push!(doc, sent)
            sent = []
        elseif line[1] == '#' # comment line
            continue
        else
            items = Vector{String}(Base.split(line,'\t'))
            push!(sent, items)
        end
    end
    length(sent) > 0 && push!(doc, sent)
    T = typeof(doc[1][1])
    features = Vector{Vector{T}}(doc)

    metadata = Dict{String,Any}("n_observations" => length(features))
    return UD_English(metadata, split, features)
end

# DEPRECATED INTERFACE, REMOVE IN v0.7 (or 0.6.x)
function Base.getproperty(::Type{UD_English}, s::Symbol)
    if s === :traindata
        @warn "UD_English.traindata() is deprecated, use `UD_English(split=:train)[]` instead." maxlog=2
        function traindata(; dir=nothing)
            UD_English(; split=:train, dir)[]
        end
        return traindata
    elseif s === :testdata
        @warn "UD_English.testdata() is deprecated, use `UD_English(split=:test)[]` instead."  maxlog=2
        function testdata(; dir=nothing)
            UD_English(; split=:test, dir)[]
        end
        return testdata
    elseif s === :devdata
        @warn "UD_English.devdata() is deprecated, use `UD_English(split=:dev)[]` instead."  maxlog=2
        function devdata(; dir=nothing)
            UD_English(; split=:dev, dir)[]
        end
        return devdata
    else
        return getfield(UD_English, s)
    end
end


text:
function euclidean_distance(set1::Array, set2::Array)

    @assert length(set1) == length(set2) "Sets don't have the same length"

    sum = zero(eltype(set1))

    for i in 1:length(set1)
            sum += (set2[i] - set1[i])^2
    end

    return sqrt(sum)

end


text:
# Simulate a falling segway
using SegwayRobot, Plots

xt = zeros(7)
xt[1:2] .= randn(2)
xt[4] = pi/12

dt = 1/25
t = 0:dt:30

anim = Animation()
for i = 1:length(t)
    # Plot robot
    fig = plot_segway_3d(xt, [-1, 19], [-10, 10])
    plot!(fig, title = "t = $(round(t[i],1))", xlab = "x", ylab = "y", size = (1920, 1080))
    frame(anim)

    # Simulate
    ut = (0.3*xt[1]+50*xt[4]+30*xt[5]+5*xt[6])*ones(2)
	#println(ut)
    segway_rk4!(xt, ut, dt)
end

# Create animation
run(`ffmpeg -v 0 -framerate $(round(Int, 1/dt)) -loop 0 -i $(anim.dir)/%06d.png -y video_pcontrol.mp4`)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using FITSIO
using CFITSIO

# Deal with compatibility issues.
using Test
using Random # for `randstring`

@testset "Images" begin
    # Create a FITS instance and loop over supported types.
    fname = tempname() * ".fits"
    FITS(fname, "w") do f
        for T in (UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64,
                  Float32, Float64)
            indata = reshape(T[1:100;], 5, 20)

            # Test writing the data to a new extension
            write(f, indata)

            # test reading the full array
            outdata = read(f[end])
            @test indata == outdata
            @test eltype(indata) == eltype(outdata)

            # test reading subsets of the array
            @test read(f[end], :, :) == indata
            @test read(f[end], 4, 1:10) == indata[4, 1:10]  # 2-d array
            @test read(f[end], :, 4) == indata[:, 4]  # 1-d array
            @test read(f[end], 2, 3) == indata[2, 3]  # scalar
            @test read(f[end], :, 1:2:10) == indata[:, 1:2:10]
            @test read(f[end], 1:3, :) == indata[1:3, :]

            # test expected errors
            @test_throws DimensionMismatch read(f[end], :)
            @test_throws DimensionMismatch read(f[end], :, :, 1)
            @test_throws BoundsError read(f[end], 1:6, :)
            @test_throws BoundsError read(f[end], 1, 0)

        end

        @test_throws ErrorException f[100]

        # Test representation
        @test repr(f)[end-17:end] == "9          Image  "
        @test repr(f[1])[1:6] == "File: "

        # test iteration
        for hdu in f
            @test size(hdu) == (5, 20)
        end
    end
    rm(fname, force=true)

    # copy_section()
    fname1 = tempname() * ".fits"
    fname2 = tempname() * ".fits"
    try
        f1 = FITS(fname1, "w")
        indata = reshape(Float32[1:400;], 20, 20)
        write(f1, indata)

        f2 = FITS(fname2, "w")
        copy_section(f1[1], f2, 1:10, 1:10)
        copy_section(f1[1], f2, 1:10, 1:2:20)
        outdata = read(f2[1])
        @test outdata == indata[1:10, 1:10]
        outdata = read(f2[2])
        @test outdata == indata[1:10, 1:2:20]
        close(f1)
        close(f2)
    finally
        rm(fname1, force=true)
        rm(fname2, force=true)
    end

    @testset "non-allocating read" begin
        fname1 = tempname() * ".fits"
        try
            FITS(fname1, "r+") do f

                indata = reshape(Float32[1:400;], 20, 20)
                write(f, indata)

                indata3D = reshape([1:30;], 3, 2, 5)
                write(f, indata3D)

                # Read the entire array using read to compare with read!
                a = read(f[1])
                a3D = read(f[2])

                # Read the entire array
                b = zeros(eltype(indata),size(indata))
                read!(f[1],b)
                @test a == b
                read!(f[1],b,:,:)
                @test a == b

                # Read a 2D slice
                b = zeros(eltype(indata),2,2)
                read!(f[1],b,1:2,1:2)
                @test a[1:2,1:2] == b

                # Read an entire 1D slice
                b = zeros(eltype(indata),size(indata,1))
                read!(f[1],b,:,1)
                @test a[:,1] == b

                b = zeros(eltype(indata),size(indata,2))
                read!(f[1],b,1,:)
                @test a[1,:] == b

                # Read a part of a 1D slice
                b = zeros(eltype(indata),1)
                read!(f[1],b,1:1,1)
                @test a[1,1] == b[1]
                read!(f[1],b,1,1:1)
                @test a[1,1] == b[1]

                # Read a single element into a 0-dim array
                b = zeros(eltype(indata))
                read!(f[1],b,1,1)
                @test a[1,1] == first(b)

                # Test for errors
                b = zeros(Float64)
                # Type is checked before dimensions
                @test_throws TypeError read!(f[1],b,1,1)
                @test_throws TypeError read!(f[1],b)

                b = zeros(eltype(indata))
                @test_throws DimensionMismatch read!(f[1],b,1:10,1)
                @test_throws DimensionMismatch read!(f[1],b,1:10,1:10)
                @test_throws DimensionMismatch read!(f[1],b)
                @test_throws DimensionMismatch read!(f[1],b,:,1)
                @test_throws DimensionMismatch read!(f[1],b,1,:)
                @test_throws DimensionMismatch read!(f[1],b,:,:)

                b = zeros(eltype(indata),1)
                @test_throws DimensionMismatch read!(f[1],b,1,1)

                b3D = zeros(eltype(indata3D),size(indata3D))
                read!(f[2],b3D)
                @test a3D == b3D
                read!(f[2],b3D,:,:,:)
                @test a3D == b3D

                b0D = zeros(eltype(indata))
                read!(f[1],b0D,1,1)
                @test a[1,1] == b0D[]

                @testset "read into views" begin

                    b = zeros(eltype(indata),size(indata))

                    # Entire array
                    b_view = @view b[:,:]
                    read!(f[1],b_view)
                    @test a == b

                    b_view = @view b3D[:,:,:]
                    read!(f[2],b_view)
                    @test a3D == b3D

                    @testset "1D slices of a 2D array" begin
                        for ax2 in axes(b,2)
                            b .= zero(eltype(b))
                            b_view = @view b[:,ax2]
                            read!(f[1],b_view,:,ax2)
                            @test a[:,ax2] == b[:,ax2]
                        end

                        # Non-contiguous views can not be read into
                        b .= zero(eltype(b))
                        b_view = @view b[1,:]
                        @test_throws ArgumentError read!(f[1],b_view,1,:)
                    end

                    @testset "1D slices of a 3D array" begin
                        for ax2 in axes(b3D,2), ax3 in axes(b3D,3)
                            b3D .= zero(eltype(b3D))
                            b_view = @view b3D[:,ax2,ax3]
                            read!(f[2],b_view,:,ax2,ax3)
                            @test a3D[:,ax2,ax3] == b3D[:,ax2,ax3]
                        end

                        # Non-contiguous views can not be read into
                        b3D .= zero(eltype(b3D))
                        b_view = @view b3D[1,1,:]
                        @test_throws ArgumentError read!(f[2],b_view,1,1,:)

                        b_view = @view b3D[1,:,1]
                        @test_throws ArgumentError read!(f[2],b_view,1,:,1)
                    end

                    @testset "2D slices of a 3D array" begin
                        for ax3 in axes(b3D,3)
                            b3D .= zero(eltype(b3D))
                            b_view = @view b3D[:,:,ax3]
                            read!(f[2],b_view,:,:,ax3)
                            @test a3D[:,:,ax3] == b3D[:,:,ax3]
                        end

                        # Non-contiguous views can not be read into
                        b3D .= zero(eltype(b3D))
                        b_view = @view b3D[:,1,:]
                        @test_throws ArgumentError read!(f[2],b_view,:,1,:)

                        b_view = @view b3D[1,:,:]
                        @test_throws ArgumentError read!(f[2],b_view,1,:,:)
                    end

                    @testset "0D view" begin
                        b .= zero(eltype(b))
                        b_view = @view b[1,1]
                        read!(f[1],b_view,1,1)
                        @test a[1,1] == b[1,1]
                    end
                end
            end
        finally
            rm(fname1, force=true)
        end
    end

    @testset "reinterpreted complex" begin
        @testset "read" begin
            fname = tempname() * ".fits"

            function _testreadcomplex(arr::AbstractArray{Complex{T}}) where T
                FITS(fname,"w") do f
                    write(f, collect(reinterpret(T,arr)) )
                end
                b = similar(arr)
                FITS(fname,"r") do f
                    read!(f[1], reinterpret(T,b))
                end
                @test b == arr
                c = FITS(fname,"r") do f
                    read(f[1])
                end
                @test c == reinterpret(T,arr)
            end

            function testreadcomplex(arr::Array{Complex{T}}) where T
                # Given a complex array, reinterpret it as float and write to file
                # Read it back and compare
                _testreadcomplex(arr)

                # Write a contiguous subsection instead of the entire array
                region = 1:size(arr,1)
                arrv = @view arr[region]
                _testreadcomplex(arrv)
            end

            try
                testreadcomplex(ones(ComplexF64,3))
                testreadcomplex(ones(ComplexF64,3,4))
                testreadcomplex(ones(ComplexF64,3,4,5))
            finally
                rm(fname,force=true)
            end
        end

        @testset "write" begin
            fname = tempname() * ".fits"

            function _testwritecomplex(arr::AbstractArray{Complex{T}}) where T
                FITS(fname,"w") do f
                    write(f, reinterpret(T,arr) )
                end
                a = FITS(fname,"r") do f
                    read(f[1])
                end
                b = collect(reinterpret(T,arr))
                @test b == a
                c = FITS(fname,"r") do f
                    reinterpret(Complex{T},read(f[1]))
                end
                @test c == arr
            end

            function testwritecomplex(arr::Array{Complex{T}}) where T
                # Given a complex array, reinterpret it as float and write to file
                # Read it back and compare
                _testwritecomplex(arr)

                # Write a subsection instead of the entire array
                region = 1:size(arr,1)
                arrv = @view arr[region]
                _testwritecomplex(arrv)
            end

            try
                testwritecomplex(ones(ComplexF64,3))
                testwritecomplex(ones(ComplexF64,3,4))
                testwritecomplex(ones(ComplexF64,3,4,5))
            finally
                rm(fname,force=true)
            end
        end
    end

    @testset "write views to file" begin
        fname = tempname() * ".fits"
        b = reshape([1:8;],2,2,2)

        @testset "2D slice of a 3D array" begin
            try
                FITS(fname,"r+") do f
                    for (ind,ax3) in enumerate(axes(b,3))
                        b_view = @view b[:,:,ax3]
                        write(f,b_view)
                        @test read(f[ind]) == b_view
                    end
                end
            finally
                rm(fname)
            end
        end

        @testset "1D slice of a 3D array" begin
            try
                FITS(fname,"r+") do f
                    ax23iter = Iterators.product(axes(b)[2:3]...)
                    for (ind,(ax2,ax3)) in enumerate(ax23iter)
                        b_view = @view b[:,ax2,ax3]
                        write(f,b_view)
                        @test read(f[ind]) == b_view
                    end
                end
            finally
                rm(fname)
            end
        end

        @testset "Non-contiguous" begin
            try
                FITS(fname,"r+") do f
                    b_view = @view b[:,1,:]
                    @test_throws ArgumentError write(f,b_view)

                    b_view = @view b[1,:,:]
                    @test_throws ArgumentError write(f,b_view)

                    b_view = @view b[1,1,:]
                    @test_throws ArgumentError write(f,b_view)

                    b_view = @view b[1,:,1]
                    @test_throws ArgumentError write(f,b_view)

                    # write something random to avoid the error on closing
                    write(f,zeros(1))
                end
            finally
                rm(fname)
            end
        end
    end
end

@testset "Write data to an existing image HDU" begin
    @testset "Overwrite an image" begin
        # Create a file with two images
        fname = tempname() * ".fits"
        FITS(fname, "w") do f
            write(f, [[1 2 3]; [4 5 6]])
            write(f, [[7 8 9]; [17 52 10]])
        end

        # Open the file in read-write mode
        FITS(fname, "r+") do f
            @test read(f[1]) == [[1 2 3]; [4 5 6]]
            @test read(f[2]) == [[7 8 9]; [17 52 10]]

            # Write data into the first image HDU
            write(f[1], [[11 12 13]; [14 15 16]])

            @test read(f[1]) == [[11 12 13]; [14 15 16]]  # First HDU is updated
            @test read(f[2]) == [[7 8 9]; [17 52 10]]  # Second HDU is untouched
        end

        rm(fname, force=true)
    end

    @testset "Show size mismatch error" begin
        fname = tempname() * ".fits"
        FITS(fname, "w") do f
            write(f, [[1 2 3 4]; [5 6 7 8]])
            write(f, [[1 2 3]; [4 5 6]])
        end

        FITS(fname, "r+") do f
            second_hdu = f[2]  # Get the second HDU

            # Make the first HDU to be the current HDU
            # Needed to check if `write` method compares data size with the given HDU
            # and not with the current HDU
            fits_movabs_hdu(f.fitsfile, 1)

            # Attempt to write data array of incorrect size to the second HDU
            exception = ErrorException("size of HDU [2, 3] not equal to size of data [2, 4].")
            @test_throws exception write(second_hdu, [[11 12 13 14]; [15 16 17 18]])
        end

        rm(fname, force=true)
    end
end

@testset "Tables" begin
    fname = tempname() * ".fits"
    f = FITS(fname, "w")

    ## Binary table
    indata = Dict{String, Array}()
    for (i, T) in enumerate([UInt8, Int8, UInt16, Int16, UInt32, Int32, Int64,
                             Float32, Float64, ComplexF32, ComplexF64])
        indata["col$i"] = T[1:20;]
    end
    i = length(indata) + 1
    indata["col$i"] = [randstring(10) for j=1:20]  # ASCIIString column
    i += 1
    indata["col$i"] = ones(Bool, 20)  # Bool column
    i += 1
    indata["col$i"] = reshape([1:40;], (2, 20))  # vector Int64 column
    i += 1
    indata["col$i"] = [randstring(5) for j=1:2, k=1:20]  # vector ASCIIString col
    indata["vcol"] = [randstring(j) for j=1:20]  # variable length column
    indata["VCOL"] = [collect(1.:j) for j=1.:20.] # variable length

    # test writing
    write(f, indata; varcols=["vcol", "VCOL"])

    # test reading
    colnames = FITSIO.colnames(f[2])
    for (colname, incol) in indata
        outcol = read(f[2], colname)  # table is in extension 2 (1 = primary hdr)
        @test outcol == incol
        @test eltype(outcol) == eltype(incol)
        @test colname in colnames
    end

    # TODO: remove the tests for deprecation warnings when we don't issue them, but keep the
    # `@test_throws` test.
    @test_throws ErrorException @test_nowarn(read(f[2], "vcol", case_sensitive=false))
    @test_nowarn read(f[2], "col2")
    @test_logs (:warn, r"case_sensitive") read(f[2], "COL2")

    # Test representation
    @test repr(f[2])[end-38:end] == "\n\n         (*) = variable-length column"

    ## ASCII tables

    indata = Dict{String, Array}()
    for (i, T) in enumerate([Int16, Int32, Float32, Float64])
        indata["col$i"] = T[1:20;]
    end
    i = length(indata) + 1
    indata["col$i"] = [randstring(10) for j=1:20]

    write(f, indata; hdutype=ASCIITableHDU)

    # For ASCII tables, the types don't round trip so we need to define the
    # expected output type for each input type.
    expected_type = Dict(Int16=>Int32, Int32=>Int32,
                         Float32=>Float64, Float64=>Float64,
                         String=>String)
    colnames = FITSIO.colnames(f[3])
    for (colname, incol) in indata
        outcol = read(f[3], colname)  # table is in extension 3
        @test outcol == incol
        @test eltype(outcol) == expected_type[eltype(incol)]
        @test colname in colnames
    end
    # test show/repr on ASCIITableHDU by checking that a couple lines are what we expect
    lines = split(repr(f[3]), "\n")
    @test lines[4] == "Rows: 20"
    @test lines[6] == "         col3  Float64  E26.17  "
    close(f)
    rm(fname, force=true)
end

@testset "FITSHeader" begin
    fname = tempname() * ".fits"
    f = FITS(fname, "w")

    # test that show() works on an empty file and that the beginning and end
    # arre what we expect.
    s = repr(f)
    @test s[1:6] == "File: "
    @test s[end-7:end] == "No HDUs."

    @test_throws ErrorException FITSHeader(["KEY1"], [1, 2, 3], ["comment 1", "comment 2"])

    inhdr = FITSHeader(["FLTKEY", "INTKEY", "BOOLKEY", "STRKEY", "COMMENT",
                        "HISTORY"],
                       [1.0, 1, true, "string value", nothing, nothing],
                       ["floating point keyword",
                        "",
                        "boolean keyword",
                        "string value",
                        "this is a comment",
                        "this is a history"])

    @test repr(inhdr) == """
FLTKEY  =                  1.0 / floating point keyword
INTKEY  =                    1
BOOLKEY =                    T / boolean keyword
STRKEY  =       'string value' / string value
COMMENT this is a comment
HISTORY this is a history"""

    inhdr["INTKEY"] = 2  # test setting by key
    inhdr[1] = 2.0  # test settting by index
    set_comment!(inhdr, "INTKEY", "integer keyword") # test setting a comment

    # Test reading possibly missing keyword
    @test_throws KeyError inhdr["BADKEY"]
    @test getkey(inhdr, "BADKEY", nothing) === nothing
    @test getkey(inhdr, "INTKEY", nothing) == inhdr["INTKEY"]

    indata = reshape(Float32[1:100;], 5, 20)
    write(f, indata; header=inhdr)

    # Write a second block.
    inhdr2 = deepcopy(inhdr)
    inhdr2["INTKEY"] = 3 # Set it to a different value.
    write(f, indata; header=inhdr2)

    outhdr = read_header(f[1])
    @test outhdr["FLTKEY"] === 2.0
    @test outhdr["INTKEY"] === 2
    @test outhdr["BOOLKEY"] === true
    @test outhdr["STRKEY"] == "string value"
    @test get_comment(outhdr, 13) == "this is a comment"
    @test get_comment(outhdr, 14) == "this is a history"
    @test length(outhdr) == 14
    @test haskey(outhdr, "FLTKEY")

    # Read entire header as a single string
    s = read_header(f[1], String)
    @test s[1:9] == "SIMPLE  ="  # all headers should start with this.
    @test length(s) == (9 + length(inhdr)) * 80  # 9 lines = 8 default + "END"

    # Test to check that read_header gets the right block even after reading another.
    s_reread = read_header(f[1])
    s_reread = read_header(f[2])
    s_reread = read_header(f[1], String)
    @test s == s_reread

    # update an existing keyword, and read it directly
    write_key(f[1], "FLTKEY", 3.0)
    @test read_key(f[1], 9) == ("FLTKEY", 3.0, "floating point keyword")
    @test read_key(f[1], "FLTKEY") == (3.0, "floating point keyword")

    # Test appending a keyword, then modifying a keyword of different
    # values with write_key()
    for value in [1.0, "string value", 42, false, nothing]
        write_key(f[1], "NEWKEY", value, "new key comment")
        @test read_key(f[1], "NEWKEY") == (value, "new key comment")
        @test read_key(f[1], 15) == ("NEWKEY", value, "new key comment")
    end

    # Test that show() works and that the beginning of output is what we expect.
    @test repr(f)[1:6] == "File: "

    # Clean up from last test.
    close(f)
    rm(fname, force=true)
end

# -----------------------------------------------------------------------------
# parsing non-standard keyword records

# `create_test_file` : Create a simple FITS file for testing, with the
# given header string added after the required keywords. The length of
# `header` must be a multiple of 80.  The purpose of creating such
# files is to test the parsing of non-standard FITS keyword records
# (non-standard files can't be created with cfitsio).

function create_test_file(fname::AbstractString, header::String)
    if length(header) % 80 != 0
        error("length of header must be multiple of 80")
    end

    f = open(fname, "w")

    stdhdr = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                   10 / length of data axis 1                          NAXIS2  =                   10 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            "
    endline = "END                                                                             "
    data = fill(0., (10, 10))  # 10x10 array of big-endian Float64 zeros

    # write header
    write(f, stdhdr)
    write(f, header)
    write(f, endline)

    # add padding
    block_position = (length(stdhdr) + length(header) + length(endline)) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, " "^padding)

    # write data
    write(f, data)

    # add padding
    block_position = sizeof(data) % 2880
    padding = (block_position == 0) ? 0 : 2880 - block_position
    write(f, fill(0x00, (padding,)))

    close(f)
end

@testset "Non-standard keywords" begin
    fname = tempname() * ".fits"
    fname2 = tempname() * ".fits"
    try
        # Create a test file with a few non-standard keyword records
        header = "        Warning: CROTA2 is inaccurate due to considerable skew                  SKEW    =  1.9511305786508E+00,  1.9234924037208E+00 /Measure of skew           "
        create_test_file(fname, header)

        # check that we can read the header (and data).
        f = FITS(fname)
        hdr = read_header(f[1])
        @test hdr["SKEW"] == "1.9511305786508E+00,"
        @test get_comment(hdr, "SKEW") == "1.9234924037208E+00 /Measure of skew"
        data = read(f[1])
        @test data == zeros(10, 10)
        close(f)

        # Test that we can read is at a memory backed file
        f = FITS(read(fname))
        hdr = read_header(f[1])
        data = read(f[1])
        @test data == zeros(10, 10)
        close(f)

        # test that we can write it back out.
        f2 = FITS(fname2, "w")
        write(f2, data; header=hdr)
        close(f2)
    finally
        rm(fname, force=true)
        rm(fname2, force=true)
    end
end

@test_deprecated FITSIO.libcfitsio_version()

# test we can still access Libcfitsio

@test begin
    using FITSIO.Libcfitsio
    Libcfitsio.libcfitsio_version() isa VersionNumber
end


text:
abstract type Optimizer end

"Calculate change in parameters for gradient descent"
update(opt::Optimizer, g_t::AbstractArray{T}) where {T<:Real} = error("not implemented")
update(opt::Optimizer, g_t::Real) = update(opt::Optimizer, [g_t])[1]

"Number of epochs run"
t(opt::Optimizer) = opt.t

optimizer(opt::Optimizer) = opt.opt_type

params(opt::Optimizer) = error("not implemented")

"Print summary"
function Base.show(io::IO, opt::Optimizer)
    print(io,"$(optimizer(opt))(t=$(t(opt::Optimizer)), $(params(opt)))")
end

"Deep copy an optimizer"
function Base.deepcopy(opt::O) where {O<:Optimizer}
    O([deepcopy(getfield(opt,f)) for f in fieldnames(O)]...)
end


text:
using GraphTools
using Base.Test

d = @graph A -> B
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n}\n" == to_dot(d)

d = @graph A -> B -> C
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n2 -> 3 \n3 [\"label\"=\"C\"]\n}\n" == to_dot(d)

d = @graph A -> B 50-> C
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n2 -> 3 [\"label\"=\"50\",\"weight\"=\"50\"]\n3 [\"label\"=\"C\"]\n}\n" == to_dot(d)

da = @graph A -> B 50-> C
db = @graph C 10-> A
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n2 -> 3 [\"label\"=\"50\",\"weight\"=\"50\"]\n3 [\"label\"=\"C\"]\n3 -> 1 [\"label\"=\"10\",\"weight\"=\"10\"]\n}\n" == to_dot(da + db)

da = @graph A
db = @graph B -> C
dc = @graph D -> A
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n2 [\"label\"=\"B\"]\n2 -> 3 \n3 [\"label\"=\"C\"]\n3 -> 2 \n4 [\"label\"=\"D\"]\n4 -> 1 \n}\n" == to_dot(da + db + dc)

d = @graph "A" -> "B"
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n}\n" == to_dot(d)

d = @graph "A" -> "B" -> "C"
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n2 -> 3 \n3 [\"label\"=\"C\"]\n}\n" == to_dot(d)

d = @graph "A" -> "B" 50-> "C"
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n2 -> 3 [\"label\"=\"50\",\"weight\"=\"50\"]\n3 [\"label\"=\"C\"]\n}\n" == to_dot(d)

da = @graph "A" -> "B" 50-> "C"
db = @graph "C" 10-> "A"
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n1 -> 2 \n2 [\"label\"=\"B\"]\n2 -> 3 [\"label\"=\"50\",\"weight\"=\"50\"]\n3 [\"label\"=\"C\"]\n3 -> 1 [\"label\"=\"10\",\"weight\"=\"10\"]\n}\n" == to_dot(da + db)

da = @graph "A"
db = @graph "B" -> "C"
dc = @graph "D" -> "A"
@test "digraph graphname {\n1 [\"label\"=\"A\"]\n2 [\"label\"=\"B\"]\n2 -> 3 \n3 [\"label\"=\"C\"]\n3 -> 2 \n4 [\"label\"=\"D\"]\n4 -> 1 \n}\n" == to_dot(da + db + dc)


text:
using Setfield, Flux, InteractiveUtils

@testset "code2lens & lens2code" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2,"c":"oh"}]}""")
	j2 = JSON.parse("""{"a": [{"a":1,"b":3,"c":"hi"},{"b":2,"a":1,"c":"Mark"}]}""")
	j3 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")
	j4 = JSON.parse("""{"a": []}""")
	j5 = JSON.parse("""{}""")

	sch = JsonGrinder.schema([j1,j2,j3,j4, j5])
	ext = suggestextractor(sch)

    for t in list_traversal(sch)
        @test all(l -> all(t .== lens2code(sch, l)), code2lens(sch, t))
    end
    for t in list_traversal(ext)
        @test all(l -> all(t .== lens2code(ext, l)), code2lens(ext, t))
    end

	l_sch = code2lens(sch, "w") |> only
	l_ext = code2lens(ext, "w") |> only
	@test l_sch == (@lens _.childs[:a].items.childs[:c])
	@test l_ext == (@lens _.dict[:a].item.dict[:c])
	@test get(sch, l_sch) == sch["w"]
	@test get(ext, l_ext) == ext["w"]
	@test lens2code(sch, l_sch) |> only == "w"
	@test lens2code(ext, l_ext) |> only == "w"

	j1 = JSON.parse("""{"a": 4, "b": "birb"}""")
	j2 = JSON.parse("""{"a": { "a": "hello", "b":[5,6]}, "b": "bird"}""")
	j3 = JSON.parse("""{"a": [1, 2, 3, "hi"], "b": "word"}""")

	sch = schema([j1, j2, j3])
	l_sch = code2lens(sch, "S") |> only
	@test l_sch == (@lens _.childs[:a].childs[3].items)
end

@testset "onehot hcat" begin
	X1 = Flux.onehotbatch([1,2,3,4,5], 1:10)
	@test @which(hcat(X1,X1)).module == JsonGrinder
	@test @which(reduce(hcat, [X1,X1])).module == JsonGrinder

	@test hcat(X1,X1) == Flux.onehotbatch([1,2,3,4,5,1,2,3,4,5], 1:10)
	@test reduce(hcat, [X1,X1]) == Flux.onehotbatch([1,2,3,4,5,1,2,3,4,5], 1:10)

	X2 = Flux.onehotbatch([1,2,3,4,5,6], 1:12)
	@test_throws DimensionMismatch hcat(X1,X2)
end


text:
let
    # test outer CuDim3 constructors
    @test CUDA.CuDim3((Cuint(4),Cuint(3),Cuint(2))) == CUDA.CuDim3(Cuint(4),Cuint(3),Cuint(2))
    @test CUDA.CuDim3((Cuint(3),Cuint(2)))          == CUDA.CuDim3(Cuint(3),Cuint(2),Cuint(1))
    @test CUDA.CuDim3((Cuint(2),))                  == CUDA.CuDim3(Cuint(2),Cuint(1),Cuint(1))
    @test CUDA.CuDim3(Cuint(2))                     == CUDA.CuDim3(Cuint(2),Cuint(1),Cuint(1))

    # outer constructor should type convert
    @test CUDA.CuDim3(2)       == CUDA.CuDim3(Cuint(2),Cuint(1),Cuint(1))
    @test_throws InexactError CUDA.CuDim3(typemax(Int64))

    # CuDim type alias should accept conveniently-typed dimensions
    @test isa(2,        CUDA.CuDim)
    @test isa((2,),     CUDA.CuDim)
    @test isa((2,2),    CUDA.CuDim)
    @test isa((2,2,2),  CUDA.CuDim)
    @test isa(Cuint(2), CUDA.CuDim)
end

@testset "device" begin

let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
    dummy = CuFunction(md, "dummy")

    # different cudacall syntaxes
    cudacall(dummy, Tuple{})
    cudacall(dummy, Tuple{}; threads=1)
    cudacall(dummy, Tuple{}; threads=1, blocks=1)
    cudacall(dummy, Tuple{}; threads=1, blocks=1, shmem=0)
    cudacall(dummy, Tuple{}; threads=1, blocks=1, shmem=0, stream=CuDefaultStream())
    cudacall(dummy, Tuple{}; threads=1, blocks=1, shmem=0, stream=CuDefaultStream(), cooperative=false)
    cudacall(dummy, ())
    cudacall(dummy, (); threads=1, blocks=1, shmem=0, stream=CuDefaultStream(), cooperative=false)

    # different launch syntaxes
    CUDA.launch(dummy)
    CUDA.launch(dummy; threads=1)
    CUDA.launch(dummy; threads=1, blocks=1)
    CUDA.launch(dummy; threads=1, blocks=1, shmem=0)
    CUDA.launch(dummy; threads=1, blocks=1, shmem=0, stream=CuDefaultStream())
    CUDA.launch(dummy; threads=1, blocks=1, shmem=0, stream=CuDefaultStream(), cooperative=false)
end

let
    md = CuModuleFile(joinpath(@__DIR__, "ptx/vectorops.ptx"))
    vadd = CuFunction(md, "vadd")
    vmul = CuFunction(md, "vmul")
    vsub = CuFunction(md, "vsub")
    vdiv = CuFunction(md, "vdiv")

    a = rand(Float32, 10)
    b = rand(Float32, 10)
    ad = CuArray(a)
    bd = CuArray(b)

    # Addition
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vadd,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a+b
    end

    # Subtraction
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vsub,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a-b
    end

    # Multiplication
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vmul,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a.*b
    end

    # Division
    let
        c = zeros(Float32, 10)
        c_d = CuArray(c)
        cudacall(vdiv,
                 (CuPtr{Cfloat},CuPtr{Cfloat},CuPtr{Cfloat}), ad, bd, c_d;
                 threads=10)
        c = Array(c_d)
        @test c ≈ a./b
    end
end

end

@testset "host" begin
    c = Condition()
    CUDA.launch() do
        notify(c)
    end
    wait(c)
end

@testset "attributes" begin

md = CuModuleFile(joinpath(@__DIR__, "ptx/dummy.ptx"))
dummy = CuFunction(md, "dummy")

val = attributes(dummy)[CUDA.FUNC_ATTRIBUTE_SHARED_SIZE_BYTES]

if CUDA.version() >= v"9.0"
    attributes(dummy)[CUDA.FUNC_ATTRIBUTE_MAX_DYNAMIC_SHARED_SIZE_BYTES] = val
end

end


text:
# Another way to include functions.
# See Parse_and_Iterate_Example for the Stan way.

using CmdStan, Distributions, Random

ProjDir = @__DIR__
cd(ProjDir)

function full_model(base_model, functions)
  initial_model = open(f->read(f, String), base_model)
  funcs = open(f->read(f, String), functions)
  "functions{\n$(funcs)}\n"*initial_model
end

function full_model(base_model, shared_functions, local_functions)
  initial_model = open(f->read(f, String), base_model)
  shared_funcs = open(f->read(f, String), shared_functions)
  local_funcs = open(f->read(f, String), local_functions)
  "functions{\n$(shared_funcs)\n$(local_funcs)}\n"*initial_model
end

model = full_model("bernoulli.stan", "shared_funcs.stan")

stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli", model=model,
  printsummary=false)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir)

if rc == 0
  show(chn)
end

println("\nOn to model2\n")

model2 = full_model("bernoulli2.stan", "shared_funcs.stan", "local_funcs.stan")

stanmodel = Stanmodel(Sample(save_warmup=true, num_warmup=1000,
  num_samples=2000, thin=1), name="bernoulli2", model=model2,
  printsummary=false)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

rc, chn, cnames = stan(stanmodel, observeddata, ProjDir)

if rc == 0
  read_summary(stanmodel)
end


text:
"""
    arlett_sim([seed])

Run the random walk simulation of Yoojin/Eleanor's experiments.

Simulation controls:
    set_FLOW_BIAS()
    set_CATALASE_ON_WALLS()
    set_NUMBER_OF_WALKS()
    set_MAX_STEPS_PER_WAL()
"""
function arlett_sim(seed::Int=randseed())
    arlett_print(seed)
    data = Dict{String,Integer}()
    data["side wall"] = 0
    data["top wall"] = 0
    data["escape"] = 0
    data["left outer sensor"] = 0
    data["left inner sensor"] = 0
    data["center sensor"] = 0
    data["right inner sensor"] = 0
    data["right outer sensor"] = 0
    data["particles unresolved"] = 0
    output = run_sim!(data, seed)
    present_arlett_data(output)
    nothing
end

"Helper function to start, not exported in Arlett module"
function arlett_print(seed::Int64)
    println("############################")
    println("      Arlett Simulation     ")
    println("############################")
    println("_________Parameters_________")
    println("# of trials:\t\t", string(NUMBER_OF_WALKS))
    println("# of steps (max):\t", string(MAX_STEPS_PER_WALK))
    println("Step size, water:\t", string(WATER_STEP_SIZE))
    println("Step size, enz:\t\t", string(ENZ_STEP_SIZE))
    println("Step size, ppd:\t\t", string(PPD_STEP_SIZE))
    println("Flow bias on:\t\t", string(FLOW_BIAS))
    println("Catalase on walls:\t", string(CATALASE_ON_WALLS))
    println("Random seed:\t\t$seed")
    nothing
end

"Helper function to end, not exported in Arlett module"
function present_arlett_data(data::Dict{String, Integer})
    println("_________Results____________")
    println("Side wall collisions:\t", data["side wall"])
    if CATALASE_ON_WALLS
        println("# caught by catalase:\t", data["top wall"])
    else
        println("Top wall collisions:\t", data["top wall"])
    end
    println("# escaped:\t\t", data["escape"])
    println("*Sensors from left to right*")
    println("# in first:\t\t", data["left outer sensor"])
    println("# in second:\t\t", data["left inner sensor"])
    println("# in third (spawn):\t", data["center sensor"])
    println("# in fourth:\t\t", data["right inner sensor"])
    println("# in fifth:\t\t", data["right outer sensor"])
    println("# unresolved:\t\t", data["particles unresolved"])
    println("Average steps taken:\t", data["avg steps taken"])
    nothing
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
# 2020.02.17, Антонюк В.А.

M(P) = [Int(sign(P[i]-P[j])) for i=1:length(P), j=1:length(P)]

function TrueIndices!(P::Vector{T}) where T<:Int
    i = length(P)-1
    while i > 0
        Curr = P[i]
        for j = i+1:length(P)
            P[j] += (Curr<=P[j])
        end
        i -= 1
    end
end

function HeuristicIterations!(R::Matrix{T}) where T<:Int
    @assert size(R)[1] == size(R)[2]
    Alt = Vector{T}()
    while !isempty(R)
        Size = size(R)[1]
        (RMin,Pos) = findmin([sum(R[i,:]) for i=1:Size])
        push!(Alt, Pos)
        R = reshape([R[i,j] for i=1:Size, j=1:Size
                              if i!=Pos && j!=Pos],(Size-1,Size-1))
    end
    TrueIndices!(Alt)
    return Alt
end

function OptimizeRanking!(P::Vector{T}, R::Matrix{T}) where T<:Int
    N = size(P)[1]
    @assert size(R)[1] == size(R)[2] == N
    k = N-1
    while k > 0
        if R[P[k],P[k+1]] > R[P[k+1],P[k]]    # 2020.04.01
            P[k], P[k+1] = P[k+1], P[k]
            k -= 1         # 2020.04.01
        end
        k -= 1
    end
    return P
end

# 2020.02.27, Антонюк В.А.

# 2020.03.07
#using LinearAlgebra
@static if VERSION >= v"0.7" using LinearAlgebra end
#

function HeuristicIterations(R::Matrix{T}) where T<:Int
    @assert size(R)[1] == size(R)[2]
    (N, _), TMax = size(R), typemax(T)
    E, Dropped = ones(T, N), Vector{T}()
    while length(Dropped) < N
        V = [(E[k]==0 ? TMax : R[k,:]⋅E) for k = 1:N]
        (_, Pos) = findmin(V)
        E[Pos] = 0
        push!(Dropped, Pos)
    end
    return Dropped
end

max_stars_repo_name: vaa-msu/ranking-tests

max_issues_repo_name: vaa-msu/ranking-tests

max_forks_repo_name: vaa-msu/ranking-tests


text:
# Test if we can transform part of a Hessenberg matrix to Schur form.
# Here we look at some edge cases.

using Test, LinearAlgebra
using ArnoldiMethod: eigenvalues, local_schurfact!, is_offdiagonal_small,
                     NotWanted

include("utils.jl")

@testset "Schur factorization" begin
    for T in (Float64, BigFloat)
        let
            # 2-by-2 matrix with distinct eigenvalues while H[2,1] != 0
            H = [one(T) 2one(T); 3one(T) 4one(T)]
            H′ = copy(H)
            Q = Matrix{T}(I, 2, 2)

            @test local_schurfact!(H′, 1, 2, Q, eps(T), 2)
            @test norm(H * Q - Q * H′) < 10eps(T)
            @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
            @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H), by=realimag)
            @test iszero(H′[2,1])
        end

        let
            # 2-by-2 matrix with distinct eigenvalues while H[2,1] = 0
            H = [one(T) 2one(T); zero(T) 4one(T)]
            H′ = copy(H)
            Q = Matrix{T}(I, 2, 2)

            @test local_schurfact!(H′, 1, 2, Q, eps(T), 2)
            @test norm(H * Q - Q * H′) < 10eps(T)
            @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
            @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H), by=realimag)
            @test iszero(H′[2,1])
        end

        let
            # 2-by-2 matrix with conjugate eigenvalues
            H = [1one(T) 4one(T); -5one(T) 3one(T)]
            H′ = copy(H)
            Q = Matrix{T}(I, 2, 2)

            @test local_schurfact!(H′, 1, 2, Q, eps(T), 2)
            @test norm(H * Q - Q * H′) < 10eps(T)
            @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
            @test sort!(eigenvalues(H′), by=realimag) ≈ sort!(eigvals(H), by=realimag)
        end
    end
        # Larger real matrix.
    for T in (Float64,) #randn(BigFloat,n,n) does not yet exist https://github.com/JuliaLang/julia/pull/35111
        let
            n = 10
            # Transforming a 1+i by n-i block of the matrix H into upper triangular form
            for i = 0 : 4
                Q = Matrix{T}(I, n, n)
                H = triu(randn(T, n, n))
                H[1+i:n-i,1+i:n-i] = normal_hessenberg_matrix(T, i+1:n-i)
                H′ = copy(H)

                # Test that the procedure has converged
                @test local_schurfact!(H′, 1+i, n-i, Q)

                for j = 1+i:n-i-1
                    t = H′[j,j] + H′[j+1,j+1]
                    d = H′[j,j] * H′[j+1,j+1] - H′[j+1,j] * H′[j,j+1]

                    # Test if subdiagonal is small. If not, check if conjugate eigenvalues.
                    @test is_offdiagonal_small(H′, j) || t^2 < 4d
                end

                # Test that the elements below the subdiagonal are 0
                @test is_hessenberg(H′)

                # Test that the partial Schur decomposition relation holds
                @test norm(H * Q - Q * H′) < 1000eps(T)

                # Test that the eigenvalues of H are the same before and after transformation
                @test sort!(eigvals(H), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
            end
        end

    # COMPLEX ARITHMETIC

        let
            n = 10
            # Transforming a 1+i by 10-i block of the matrix H into upper triangular form
            for i = 0 : 4
                Q = Matrix{Complex{T}}(I, n, n)
                H = triu(randn(Complex{T}, n, n))
                H[i+1:n-i,i+1:n-i] = normal_hessenberg_matrix(Complex{T}, (i+1:n-i) .* (1 + im))
                H′ = copy(H)

                # Test that the procedure has converged
                @test local_schurfact!(H′, 1+i, n-i, Q)

                # Test if subdiagonal is small.
                for j = 1+i:n-i-1
                    @test iszero(H′[j+1,j])
                end

                # Test that the elements below the subdiagonal are 0
                @test is_hessenberg(H′)

                # Test that the partial Schur decomposition relation holds
                @test norm(H * Q - Q * H′) < 1000eps(T)

                # Test that the eigenvalues of H are the same before and after transformation
                @test sort!(eigvals(H), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
            end
        end
    end
end

@testset "Schur with nearly repeated eigenvalues" begin
    # Matrix with nearly repeated eigenpair could converge very slowly or
    # stagnate completely when there's a tiny perturbation in the computed
    # shift.
    mat(ε) = [2   0    0  ;
              5ε  1-ε  2ε ;
              0   3ε   1+ε]
    for T in (Float64,BigFloat)
        @test local_schurfact!(mat(eps(T)))
    end
end


text:
mutable struct RHS_IIF1M_Scalar{F,CType,tType,P} <: Function
  f::F
  t::tType
  dt::tType
  tmp::CType
  p::P
end

function (f::RHS_IIF1M_Scalar)(resid,u)
  resid[1] = u[1] - f.tmp - f.dt*f.f.f2(u[1],f.p,f.t+f.dt)[1]
end

mutable struct RHS_IIF2M_Scalar{F,CType,tType,P} <: Function
  f::F
  t::tType
  dt::tType
  tmp::CType
  p::P
end

function (f::RHS_IIF2M_Scalar)(resid,u)
  resid[1] = u[1] - f.tmp - 0.5f.dt*f.f.f2(u[1],f.p,f.t+f.dt)[1]
end

@muladd function initialize!(integrator,cache::Union{IIF1MConstantCache,IIF2MConstantCache,IIF1MilConstantCache},f=integrator.f)
  cache.uhold[1] = integrator.uprev
end

@muladd function perform_step!(integrator,cache::Union{IIF1MConstantCache,IIF2MConstantCache,IIF1MilConstantCache},f=integrator.f)
  @unpack t,dt,uprev,u,W,p = integrator
  @unpack uhold,rhs,nl_rhs = cache
  alg = unwrap_alg(integrator, true)
  A = integrator.f.f1(u,p,t)
  if typeof(cache) <: IIF1MilConstantCache
    error("Milstein correction does not work.")
  elseif typeof(cache) <: IIF1MConstantCache
    tmp = exp(A*dt)*(uprev + integrator.g(uprev,p,t)*W.dW)
  elseif typeof(cache) <: IIF2MConstantCache
    tmp = exp(A*dt)*(uprev + 0.5dt*integrator.f.f2(uprev,p,t) + integrator.g(uprev,p,t)*W.dW)
  end

  if integrator.iter > 1 && !integrator.u_modified
    uhold[1] = current_extrapolant(t+dt,integrator)
  end # else uhold is previous value.

  rhs.t = t
  rhs.dt = dt
  rhs.tmp = tmp
  nlres = alg.nlsolve(nl_rhs,uhold)

  u = nlres[1]
  integrator.u = u
end

mutable struct RHS_IIF1{F,uType,tType,DiffCacheType,SizeType,P} <: Function
  f::F
  tmp::uType
  t::tType
  dt::tType
  dual_cache::DiffCacheType
  sizeu::SizeType
  p::P
end
function (f::RHS_IIF1)(resid,u)
  _du = get_du(f.dual_cache, eltype(u))
  du = reinterpret(eltype(u),_du)
  f.f.f2(du,reshape(u,f.sizeu),f.p,f.t+f.dt)
  @.. resid = u - f.tmp - f.dt*du
end

mutable struct RHS_IIF2{F,uType,tType,DiffCacheType,SizeType,P} <: Function
  f::F
  tmp::uType
  t::tType
  dt::tType
  dual_cache::DiffCacheType
  sizeu::SizeType
  p::P
end
function (f::RHS_IIF2)(resid,u)
  _du = get_du(f.dual_cache, eltype(u))
  du = reinterpret(eltype(u),_du)
  f.f.f2(du,reshape(u,f.sizeu),f.p,f.t+f.dt)
  @.. resid = u - f.tmp - 0.5f.dt*du
end

@muladd function perform_step!(integrator,cache::Union{IIF1MCache,IIF2MCache},f=integrator.f)
  @unpack rtmp1,rtmp2,rtmp3,tmp,noise_tmp = cache
  @unpack uhold,rhs,nl_rhs = cache
  @unpack t,dt,uprev,u,W,p = integrator
  alg = unwrap_alg(integrator, true)
  uidx = eachindex(u)

  integrator.g(rtmp2,uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    rmul!(rtmp2,W.dW) # rtmp2 === rtmp3
  else
    mul!(rtmp3,rtmp2,W.dW)
  end

  rtmp3 .+= uprev

  if typeof(cache) <: IIF2MCache
    integrator.f.f2(rtmp1,uprev,p,t)
    dto2 = dt / 2
    @.. rtmp3 = dto2 * rtmp1 + rtmp3
  end

  A = integrator.f.f1(rtmp1,uprev,p,t)
  M = exp(A*dt)
  mul!(tmp,M,rtmp3)

  if integrator.iter > 1 && !integrator.u_modified
    current_extrapolant!(uhold,t+dt,integrator)
  end # else uhold is previous value.

  rhs.t = t
  rhs.dt = dt
  rhs.tmp = tmp
  rhs.sizeu = size(u)
  nlres = alg.nlsolve(nl_rhs,uhold)

  copyto!(uhold,nlres)

end

@muladd function perform_step!(integrator,cache::IIF1MilCache,f=integrator.f)
  @unpack rtmp1,rtmp2,rtmp3,tmp,noise_tmp = cache
  @unpack uhold,rhs,nl_rhs = cache
  @unpack t,dt,uprev,u,W,p = integrator
  alg = unwrap_alg(integrator, true)

  dW = W.dW; sqdt = integrator.sqdt
  f = integrator.f; g = integrator.g

  A = integrator.f.f1(t,uprev,rtmp1)
  M = exp(A*dt)

  uidx = eachindex(u)
  integrator.g(rtmp2,uprev,p,t)
  if typeof(cache) <: Union{IIF1MCache,IIF2MCache}
    if is_diagonal_noise(integrator.sol.prob)
      rmul!(rtmp2,W.dW) # rtmp2 === rtmp3
    else
      mul!(rtmp3,rtmp2,W.dW)
    end
  else #Milstein correction
    rtmp2 = M*rtmp2 # mul!(rtmp2,M,gtmp)
    @unpack gtmp,gtmp2 = cache
    #error("Milstein correction does not work.")
    mul!(rtmp3,rtmp2,W.dW)
    I = zeros(length(dW),length(dW));
    Dg = zeros(length(dW),length(dW)); mil_correction = zeros(length(dW))
    mil_correction .= 0.0
    for i=1:length(dW),j=1:length(dW)
        I[j,i] = 0.5*dW[i]*dW[j]
        j == i && (I[i,i] -= 0.5*dt) # Ito correction
    end
    for j = 1:length(uprev)
      #Kj = uprev .+ dt.*du1 + sqdt*rtmp2[:,j] # This works too
      Kj = uprev .+ sqdt*rtmp2[:,j]
      g(gtmp,Kj,p,t); mul!(gtmp2,M,gtmp)
      Dgj = (gtmp2 - rtmp2)/sqdt
      mil_correction .+= Dgj*I[:,j]
    end
    rtmp3 .+= mil_correction
  end

  if typeof(cache) <: IIF2MCache
    integrator.f.f2(t,uprev,rtmp1)
    dto2 = dt / 2
    @.. rtmp1 = dto2 * rtmp1 + uprev + rtmp3
    mul!(tmp,M,rtmp1)
  elseif !(typeof(cache) <: IIF1MilCache)
    @.. rtmp1 = uprev + rtmp3
    mul!(tmp,M,rtmp1)
  else
    mul!(tmp,M,uprev)
    tmp .+= rtmp3
  end

  if integrator.iter > 1 && !integrator.u_modified
    current_extrapolant!(uhold,t+dt,integrator)
  end # else uhold is previous value.

  rhs.t = t
  rhs.dt = dt
  rhs.tmp = tmp
  rhs.sizeu = size(u)
  nlres = alg.nlsolve(nl_rhs,uhold)

  copyto!(uhold,nlres)
end


text:
# tests for the environment type and related functions.

# using Base.Test
# using PyPlot

# push!(LOAD_PATH, ".")
# include("includes.jl")

function get_test_uavs()
    actions = [-1., 1., 0.]
    joint_actions = ones(3, 3)
    joint_actions[:, 1] = -1

    r = collect(1:10)
    q_table = zeros(10^7, size(joint_actions, 2))
    q_grid = RectangleGrid(r, r, r, r, r, r, r)
    policy = PolicyTable(q_table, q_grid, 10^7)

    c_table = zeros(10, 10, 10, 10, 10)
    c_grid = RectangleGrid(r * 100, r, r, r, r)
    coord = CoordinationTable(c_table, c_grid, 10^5)

    strategy = coordinated_strategy
    act_select = "best"

    rng = MersenneTwister(1)
    sensor_range_std = 15.
    sensor_theta_std = deg2rad(1.)
    sensor_psi_std = deg2rad(1.)
    sensor_velocity_std = 1.
    use_sensor_noise = false

    uav1 = UAV(1, actions, joint_actions, policy, coord, false, strategy, act_select, rng, sensor_range_std, sensor_theta_std, sensor_psi_std, sensor_velocity_std, use_sensor_noise)
    uav2 = UAV(2, actions, joint_actions, policy, coord, false, strategy, act_select, rng, sensor_range_std, sensor_theta_std, sensor_psi_std, sensor_velocity_std, use_sensor_noise)

    return [uav1, uav2]
end

function get_test_env()
    max_steps = 100
    use_noise = true
    uavs = get_test_uavs()
    gauss_sigma = 1.0
    gauss_mu = 0.0
    d = Dynamics(3, 1., 4., 5., 5., false, MersenneTwister(1))
    start_state_generator = start_head_on_encounter
    rng = MersenneTwister(1)
    coord_scheme = 1
    e = Environment(max_steps, use_noise, uavs, gauss_sigma, gauss_mu, d, coord_scheme, start_state_generator, 1/4, 1/3, 100., pi / 32, 153., true, rng, rng)
    return e
end

# Environment tests
function test_env_constructor()
    get_test_env()
end

function test_get_uav_actions()
    env = get_test_env()
    uav_state_own = UAVState(ones(Float64, 8)...)
    uav_state_int = UAVState(zeros(Float64, 8)...)
    s = State([uav_state_own; uav_state_int])
    enc = Encounter()
    actions = get_uav_actions!(env, s, enc)
    @test length(actions.actions) == s.num_uavs
end

function test_step()
    env = get_test_env()
    uav_state_own = UAVState(ones(Float64, 8)...)
    uav_state_int = UAVState(ones(Float64, 8)...)
    s = State([uav_state_own; uav_state_int])
    uav_actions = UAVActions(2)
    uav_actions.actions = [deg2rad(10); deg2rad(10)]
    enc = Encounter()
    next_state = step!(env, s, uav_actions, enc)
end

function test_add_noise()
    env = get_test_env()
    actions = [0., 0.]
    add_noise!(env, actions)
    @test actions != [0., 0.]
end

function test_start_encounter_head_on()
    env = get_test_env()
    env.start_state_generator = start_head_on_encounter
    num_runs = 500
    int_xs = zeros(num_runs)
    int_ys = zeros(num_runs)
    angles = zeros(num_runs)
    for idx in 1:num_runs
        state = start_encounter(env)
        int_xs[idx] = state.uav_states[2].x
        int_ys[idx] = state.uav_states[2].y
        angles[idx] = state.uav_states[2].psi
    end
    scatter(int_xs, int_ys, c=angles)
    title("position and heading of intruder uav at start of episode")
    if SHOW_PLOTS
        show()
    end
end

function test_start_encounter_is_deterministic()
    env_1 = get_test_env()
    state_1 = start_encounter(env_1)
    env_2 = get_test_env()
    state_2 = start_encounter(env_2)
    env_3 = get_test_env()
    state_3 = start_encounter(env_3)
    @test state_1 == state_2 == state_3
    state_12 = start_encounter(env_1)
    state_22 = start_encounter(env_2)
    @test state_12 == state_22
end

function test_start_encounter_towards_origin()
    env = build_environment()
    env.start_state_generator = start_towards_origin_encounter
    num_runs = 300
    own_xs = zeros(num_runs)
    own_ys = zeros(num_runs)
    int_xs = zeros(num_runs)
    int_ys = zeros(num_runs)
    for idx in 1:num_runs
        state = start_encounter(env)
        own_xs[idx] = state.uav_states[1].x
        own_ys[idx] = state.uav_states[1].y
        int_xs[idx] = state.uav_states[2].x
        int_ys[idx] = state.uav_states[2].y

        # xlim([-2000, 2000])
        # ylim([-2000, 2000])
        # scatter(own_xs[idx], own_ys[idx], c="blue")
        # scatter(int_xs[idx], int_ys[idx], c="red")
        # show()
        # readline()
        # close()
    end
    scatter(own_xs, own_ys, c="blue")
    scatter(int_xs, int_ys, c="red")
    title("position and heading of uavs at start of episode")
    if SHOW_PLOTS
        show()
    end
end

function plot_start_states(env, num_runs, c1, c2, extra_random)
    if extra_random
        rand(env.rng, 1:10)
    end

    own_xs = zeros(num_runs)
    own_ys = zeros(num_runs)
    int_xs = zeros(num_runs)
    int_ys = zeros(num_runs)
    for idx in 1:num_runs
        state = start_encounter(env)
        own_xs[idx] = state.uav_states[1].x
        own_ys[idx] = state.uav_states[1].y
        int_xs[idx] = state.uav_states[2].x
        int_ys[idx] = state.uav_states[2].y
    end
    scatter(own_xs, own_ys, c=c1)
    scatter(int_xs, int_ys, c=c2)
    title("position and heading of uavs at start of episode")

end

function test_start_encounter_towards_origin_coordination_vs_greedy()
    num_runs = 100

    env = build_environment()
    env.start_state_generator = start_towards_origin_encounter
    env.uavs[1].strategy = coordinated_strategy
    plot_start_states(env, num_runs, "blue", "red", false)

    env = build_environment()
    env.start_state_generator = start_towards_origin_encounter
    env.uavs[1].strategy = greedy_strategy
    plot_start_states(env, num_runs, "green", "purple", true)

    if SHOW_PLOTS
        show()
    end

end
tic()
# Environment tests
test_env_constructor()

# uav and dynamics tests
test_get_uav_actions()
test_step()
test_add_noise()
toc()


text:
using Test, BFloat16s

@test BFloat16(1)   <  BFloat16(2)
@test BFloat16(1f0) <  BFloat16(2f0)
@test BFloat16(1.0) <  BFloat16(2.0)
@test BFloat16(1)   <= BFloat16(2)
@test BFloat16(1f0) <= BFloat16(2f0)
@test BFloat16(1.0) <= BFloat16(2.0)
@test BFloat16(2)   >  BFloat16(1)
@test BFloat16(2f0) >  BFloat16(1f0)
@test BFloat16(2)   >= BFloat16(1)
@test BFloat16(2f0) >= BFloat16(1f0)
@test BFloat16(2.0) >= BFloat16(1.0)

@test abs(BFloat16(-10)) == BFloat16(10)
@test Float32(BFloat16(10)) == 1f1
@test Float64(BFloat16(10)) == 10.0
@test BFloat16(2) ^ BFloat16(4) == BFloat16(16)
@test sqrt(BFloat16(4f0)) == BFloat16(2f0)


text:
module BasisFunctions

    export Basis, PolynomialBasis, eval_basis, deriv_basis, int_basis

    include("basis_functions/basis_functions.jl")

    export vandermonde_matrix, vandermonde_matrix_inverse

    include("basis_functions/vandermonde_matrix.jl")

    export BernsteinBasis, BernsteinBasisModified, LagrangeBasis,
           LegendreBasis, LegendreBasisHierarchical

    include("basis_functions/bernstein_basis.jl")
    include("basis_functions/bernstein_basis_modified.jl")
    include("basis_functions/lagrange_basis.jl")
    include("basis_functions/legendre_basis.jl")
    include("basis_functions/legendre_basis_hierarchical.jl")

    export Polynomial, BernsteinPolynomial, LagrangePolynomial

    include("basis_functions/polynomial.jl")

    import ..CommonFunctions

    CommonFunctions.nbasis(b::Basis) = nbasis(b)
    CommonFunctions.nnodes(b::Basis) = nnodes(b)
    CommonFunctions.degree(b::Basis) = degree(b)
    CommonFunctions.nodes(b::Basis)  = nodes(b)

    CommonFunctions.evaluate(b::Basis, x, y) = eval_basis(b, x, y)
    CommonFunctions.derivative(b::Basis, x, y) = deriv_basis(b, x, y)
    CommonFunctions.integral(b::Basis, x, y) = int_basis(b, x, y)

    CommonFunctions.evaluate!(pol::Polynomial, x, y) = evaluate!(pol, x, y)
    CommonFunctions.evaluate!(b::LagrangeBasis, c, x, y) = evaluate!(b, c, x, y)

end


text:
#ensure multiple calls to MPIComm works
MPIComm(Bool, Bool, Bool)
MPIComm(UInt64, UInt16, UInt32)

@test 4 == numProc(comm)
@test isa(numProc(comm), UInt16)

@test 1 <= myPid(comm) <= 4
@test isa(myPid(comm), UInt16)

@test [1, 8] == broadcastAll(comm, [myPid(comm), 8], 1)
@test [2, 5] == broadcastAll(comm, [myPid(comm), 5], 2)
@test [3, 7] == broadcastAll(comm, [myPid(comm), 7], 3)
@test [4, 6] == broadcastAll(comm, [myPid(comm), 6], 4)

@test [1, 2, 3, 4] == gatherAll(comm, [myPid(comm)])
@test ([1, 2, 3, 2, 4, 6, 3, 6, 9, 4, 8, 12]
        == gatherAll(comm, [myPid(comm), myPid(comm)*2, myPid(comm)*3]))

#check for hangs and such, hard to test if all processes are at the same spot
barrier(comm)

@test [10] == sumAll(comm, [myPid(comm)])
@test [32, 12, 10, 8] == sumAll(comm, [8, 3, myPid(comm), 2])

@test [4] == maxAll(comm, [myPid(comm)])
@test [4, -1, 8] == maxAll(comm, [myPid(comm), -Int(myPid(comm)), 8])

@test [1] == minAll(comm, [myPid(comm)])
@test [1, -4, 6] == minAll(comm, [myPid(comm), -Int(myPid(comm)), 6])

@test [sum(1:myPid(comm))] == scanSum(comm, [myPid(comm)])
@test ([myPid(comm)*5, sum(-2:-2:(-2*Int(myPid(comm)))), myPid(comm)*3]
        == scanSum(comm, [5, -2*Int(myPid(comm)), 3]))

#test distributor

dist = createDistributor(comm)
@test isa(dist, Distributor{UInt64, UInt16, UInt32})

#check for error when not waiting
@test_throws InvalidStateError resolveWaits(dist)
@test_throws InvalidStateError resolveReverseWaits(dist)

@test 4 == createFromSends(dist, [1, 2, 3, 4])

resolvePosts(dist, [pid, 2*pid, 3*pid, 4*pid])
@test_throws InvalidStateError resolveReverseWaits(dist)
@test pid*[1, 2, 3, 4] == resolveWaits(dist)

#check for error when not waiting
@test_throws InvalidStateError resolveWaits(dist)
@test_throws InvalidStateError resolveReverseWaits(dist)

#test distributor when elements not blocked by processor
dist2 = createDistributor(comm)
@test 8 == createFromSends(dist, [1, 2, 3, 4, 1, 2, 3, 4])

@test (reduce(vcat, [[(pid-1)*5+j, (pid+3)*5+j] for j in 1:4]; init=[])
        == resolve(dist, [pid, 5+pid, 10+pid, 15+pid, 20+pid, 25+pid, 30+pid, 35+pid]))

#test distributore createFromRecvs
@test ([(pid-1)*5+j for j in 1:4], [1, 2, 3, 4]) == createFromRecvs(dist, [pid, 5+pid, 10+pid, 15+pid], [1, 2, 3, 4])