text:
abstract type AbstractMotionTask end

struct SpatialAccelerationTask <: AbstractMotionTask
    path::TreePath{RigidBody{Float64}, Joint{Float64}}
    jacobian::GeometricJacobian{Matrix{Float64}}
    desired::Base.RefValue{SpatialAcceleration{Float64}}

    function SpatialAccelerationTask(
                mechanism::Mechanism, # TODO: would be nice to get rid of this; possible if compact Jacobians were available
                path::TreePath{RigidBody{Float64}, Joint{Float64}};
                frame::CartesianFrame3D = default_frame(target(path)))
        nv = num_velocities(mechanism)
        bodyframe = default_frame(target(path))
        baseframe = default_frame(source(path))
        jacobian = GeometricJacobian(bodyframe, baseframe, frame, zeros(3, nv), zeros(3, nv))
        desired = Ref{SpatialAcceleration{Float64}}(zero(SpatialAcceleration{Float64}, bodyframe, baseframe, frame))
        new(path, jacobian, desired)
    end
end

dimension(task::SpatialAccelerationTask) = 6

function setdesired!(task::SpatialAccelerationTask, desired::SpatialAcceleration)
    @framecheck task.desired[].body desired.body
    @framecheck task.desired[].base desired.base
    @framecheck task.desired[].frame desired.frame
    task.desired[] = desired
    nothing
end

function task_error(task::SpatialAccelerationTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    J = Parameter(task.jacobian, qpmodel) do jac
        world_to_desired = inv(transform_to_root(state, task.desired[].frame))
        geometric_jacobian!(jac, state, task.path, world_to_desired)
    end
    J̇v = Parameter{SpatialAcceleration{Float64}}(qpmodel) do
        bias = -bias_acceleration(state, source(task.path)) + bias_acceleration(state, target(task.path))
        transform(state, bias, task.desired[].frame)
    end
    desired = Parameter{SpatialAcceleration{Float64}}(() -> task.desired[], qpmodel)
    @expression [
        angular(J) * v̇ + angular(J̇v) - angular(desired);
        linear(J) * v̇ + linear(J̇v) - linear(desired)]
end

struct AngularAccelerationTask <: AbstractMotionTask
    path::TreePath{RigidBody{Float64}, Joint{Float64}}
    jacobian::GeometricJacobian{Matrix{Float64}}
    desired::Base.RefValue{FreeVector3D{SVector{3, Float64}}}

    function AngularAccelerationTask(
                mechanism::Mechanism, # TODO: would be nice to get rid of this; possible if compact Jacobians were available
                path::TreePath{RigidBody{Float64}, Joint{Float64}};
                frame::CartesianFrame3D = default_frame(target(path)))
        nv = num_velocities(mechanism)
        bodyframe = default_frame(target(path))
        baseframe = default_frame(source(path))
        jacobian = GeometricJacobian(bodyframe, baseframe, frame, zeros(3, nv), zeros(3, nv))
        desired = Ref(FreeVector3D(frame, 0.0, 0.0, 0.0))
        new(path, jacobian, desired)
    end
end

dimension(task::AngularAccelerationTask) = 3

function setdesired!(task::AngularAccelerationTask, desired::FreeVector3D)
    @framecheck task.desired[].frame desired.frame
    task.desired[] = desired
    nothing
end

function task_error(task::AngularAccelerationTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    J = Parameter(task.jacobian, qpmodel) do jac
        world_to_desired = inv(transform_to_root(state, task.desired[].frame))
        geometric_jacobian!(jac, state, task.path, world_to_desired)
    end
    J̇v = Parameter{SpatialAcceleration{Float64}}(qpmodel) do
        bias = -bias_acceleration(state, source(task.path)) + bias_acceleration(state, target(task.path))
        transform(state, bias, task.desired[].frame)
    end
    desired = Parameter{SVector{3, Float64}}(() -> task.desired[].v, qpmodel)
    @expression angular(J) * v̇ + angular(J̇v) - desired
end

struct LinearAccelerationTask <: AbstractMotionTask
    path::TreePath{RigidBody{Float64}, Joint{Float64}}
    jacobian::GeometricJacobian{Matrix{Float64}}
    desired::Base.RefValue{FreeVector3D{SVector{3, Float64}}}

    function LinearAccelerationTask(
                mechanism::Mechanism,
                path::TreePath{RigidBody{Float64}, Joint{Float64}};
                frame::CartesianFrame3D = default_frame(target(path)))
        nv = num_velocities(mechanism)
        bodyframe = default_frame(target(path))
        baseframe = default_frame(source(path))
        jacobian = GeometricJacobian(bodyframe, baseframe, frame, zeros(3, nv), zeros(3, nv))
        desired = Ref(FreeVector3D(frame, 0.0, 0.0, 0.0))
        new(path, jacobian, desired)
    end
end

dimension(task::LinearAccelerationTask) = 3

function setdesired!(task::LinearAccelerationTask, desired::FreeVector3D)
    @framecheck task.desired[].frame desired.frame
    task.desired[] = desired
    nothing
end

function task_error(task::LinearAccelerationTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    J = Parameter(task.jacobian, qpmodel) do jac
        world_to_desired = inv(transform_to_root(state, task.desired[].frame))
        geometric_jacobian!(jac, state, task.path, world_to_desired)
    end
    J̇v = Parameter{SpatialAcceleration{Float64}}(qpmodel) do
        bias = -bias_acceleration(state, source(task.path)) + bias_acceleration(state, target(task.path))
        transform(state, bias, task.desired[].frame)
    end
    desired = Parameter{SVector{3, Float64}}(() -> task.desired[].v, qpmodel)
    @expression linear(J) * v̇ + linear(J̇v) - desired
end

struct PointAccelerationTask <: AbstractMotionTask
    path::TreePath{RigidBody{Float64}, Joint{Float64}}
    jacobian::PointJacobian{Matrix{Float64}}
    point::Point3D{SVector{3, Float64}}
    desired::Base.RefValue{FreeVector3D{SVector{3, Float64}}}

    function PointAccelerationTask(
                mechanism::Mechanism,
                path::TreePath{RigidBody{Float64}, Joint{Float64}},
                point::Point3D)
        nv = num_velocities(mechanism)
        bodyframe = default_frame(target(path))
        baseframe = default_frame(source(path))
        @framecheck point.frame bodyframe
        jacobian = PointJacobian(baseframe, zeros(3, nv))
        desired = Ref(FreeVector3D(baseframe, 0.0, 0.0, 0.0))
        new(path, jacobian, point, desired)
    end
end

dimension(task::PointAccelerationTask) = 3

function setdesired!(task::PointAccelerationTask, desired::FreeVector3D)
    @framecheck task.desired[].frame desired.frame
    task.desired[] = desired
    nothing
end

function task_error(task::PointAccelerationTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    frame = task.desired[].frame
    state_param = Parameter(identity, state, qpmodel)
    point_in_task_frame = @expression transform(state_param, task.point, frame)
    J = Parameter(task.jacobian, qpmodel) do jac
        point_jacobian!(jac, state, task.path, point_in_task_frame())
        jac
    end
    J̇v = Parameter{SpatialAcceleration{Float64}}(qpmodel) do
        bias = -bias_acceleration(state, source(task.path)) + bias_acceleration(state, target(task.path))
        transform(state, bias, frame)
    end
    desired = Parameter{SVector{3, Float64}}(() -> task.desired[].v, qpmodel)
    T = @expression transform(state_param, relative_twist(state_param, target(task.path), source(task.path)), frame)
    @framecheck T().frame frame
    ω = @expression angular(T)
    ṗ = @expression point_velocity(T, point_in_task_frame)
    @expression ω × ṗ.v + J.J * v̇ + (angular(J̇v) × point_in_task_frame.v + linear(J̇v)) - desired
end

struct JointAccelerationTask{JT<:JointType{Float64}} <: AbstractMotionTask
    joint::Joint{Float64, JT}
    desired::Vector{Float64}

    function JointAccelerationTask(joint::Joint{Float64, JT}) where {JT<:JointType{Float64}}
        new{JT}(joint, zeros(num_velocities(joint)))
    end
end

dimension(task::JointAccelerationTask) = length(task.desired)
setdesired!(task::JointAccelerationTask, desired) = set_velocity!(task.desired, task.joint, desired)

function task_error(task::JointAccelerationTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    desired = Parameter(identity, task.desired, qpmodel)
    v̇joint = v̇[velocity_range(state, task.joint)]
    @expression v̇joint - desired
end

struct MomentumRateTask <: AbstractMotionTask
    momentum_matrix::MomentumMatrix{Matrix{Float64}}
    desired::Base.RefValue{Wrench{Float64}}

    function MomentumRateTask(mechanism::Mechanism, centroidalframe::CartesianFrame3D)
        nv = num_velocities(mechanism)
        momentum_matrix = MomentumMatrix(centroidalframe, zeros(3, nv), zeros(3, nv))
        desired = Ref(zero(Wrench{Float64}, centroidalframe))
        new(momentum_matrix, desired)
    end
end

function momentum_rate_task_params(task, qpmodel, state, v̇)
    # TODO: repeated computation of world_to_centroidal, but running into inference issues if that computation
    # is extracted out into its own Parameter
    centroidalframe = task.momentum_matrix.frame
    A = Parameter(task.momentum_matrix, qpmodel) do A
        com = center_of_mass(state)
        centroidal_to_world = Transform3D(centroidalframe, com.frame, com.v)
        world_to_centroidal = inv(centroidal_to_world)
        momentum_matrix!(A, state, world_to_centroidal)
    end
    Ȧv = Parameter{Wrench{Float64}}(qpmodel) do
        com = center_of_mass(state)
        centroidal_to_world = Transform3D(centroidalframe, com.frame, com.v)
        world_to_centroidal = inv(centroidal_to_world)
        transform(momentum_rate_bias(state), world_to_centroidal)
    end
    A, Ȧv
end

dimension(task::MomentumRateTask) = 6

function setdesired!(task::MomentumRateTask, desired::Wrench)
    @framecheck task.momentum_matrix.frame desired.frame
    task.desired[] = desired
end

function task_error(task::MomentumRateTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    A, Ȧv = momentum_rate_task_params(task, qpmodel, state, v̇)
    desired = Parameter{Wrench{Float64}}(() -> task.desired[], qpmodel)
    @expression [
        angular(A) * v̇ + angular(Ȧv) - angular(desired);
        linear(A) * v̇ + linear(Ȧv) - linear(desired)]
end

struct LinearMomentumRateTask <: AbstractMotionTask
    momentum_matrix::MomentumMatrix{Matrix{Float64}}
    desired::Base.RefValue{FreeVector3D{SVector{3, Float64}}}

    function LinearMomentumRateTask(mechanism::Mechanism, centroidalframe::CartesianFrame3D = CartesianFrame3D())
        nv = num_velocities(mechanism)
        momentum_matrix = MomentumMatrix(centroidalframe, zeros(3, nv), zeros(3, nv))
        desired = Ref(FreeVector3D(centroidalframe, 0.0, 0.0, 0.0))
        new(momentum_matrix, desired)
    end
end

dimension(task::LinearMomentumRateTask) = 3

function setdesired!(task::LinearMomentumRateTask, desired::FreeVector3D)
    @framecheck task.momentum_matrix.frame desired.frame
    task.desired[] = desired
end

function task_error(task::LinearMomentumRateTask, qpmodel, state::MechanismState, v̇::AbstractVector{Parametron.Variable})
    A, Ȧv = momentum_rate_task_params(task, qpmodel, state, v̇)
    desired = Parameter{SVector{3, Float64}}(() -> task.desired[].v, qpmodel)
    @expression linear(A) * v̇ + linear(Ȧv) - desired
end


text:
include("plants.jl")
include("control_systems.jl")
include("feedback_systems.jl")

ext: jl


text:
export Token
export ParsingError

"""
    struct Token

Holds type, lexme and linenumber
"""
struct Token
    type::Symbol
    lexme::String
    line::Int
end

Token(s::Symbol, lexme="", ln=0) = Token(s, lexme, ln)
Token(s::Symbol, ln::Int) = Token(s, "", ln)
Token(s::Symbol, ::Nothing) = Token(s, "", 0)

Base.isapprox(t1::Token, t2::Token) = t1.type == t2.type && t1.lexme == t2.lexme

function Base.show(io::IO, t::Token)
    print(io, t.type)
    isempty(t.lexme) || print(io, "\"", t.lexme, "\"")
    # iszero(t.line) || print(io, " (#", t.line, ")")
end

"""
   struct TokenStream

Holds a stream of tokens (Vector of Channel) together with
a position information.
"""
mutable struct TokenStream{T<:Union{Channel{Token}, Vector{Token}}}
    data::T
    pos::Int
    TokenStream(data) = new{typeof(data)}(data, 0)
end

TokenStream(s::String) = TokenStream(tokenize(s))

function Base.getindex(ts::TokenStream{<:Channel}, idx)
    c = ts.data
    lock(c) # aquire `cond_take` lock, i.e. nobody is allowed to take
    try
        if idx <= length(c.data)
            return c.data[idx]
        end
        # if there arn't enough elements yet just wait
        while isopen(c) && idx > length(c.data)
            wait(c.cond_take) # cond_take gets notified on every put
        end

        return c.data[idx] # finally return the data
    finally
        unlock(c)
    end
end

Base.getindex(ts::TokenStream{<:Vector}, idx) = ts.data[idx]

Base.position(ts::TokenStream) = ts.pos

function setposition!(ts::TokenStream, pos)
    ts.pos = pos
    return ts
end

reset!(ts::TokenStream) = setposition!(ts, 0)

"""
    next!(ts::TokenStream)

Get the next! element in TokenStream (increases position by one)
"""
function next!(ts::TokenStream)
    try
        token = ts[ts.pos + 1]
        ts.pos += 1
        return token
    catch e
        e isa BoundsError && return nothing
        rethrow()
    end
end

"""
    peek(ts::TokenStream)

Get the next element in TokenStream (without advancing position)
"""
function peek(ts::TokenStream)
    try
        token = ts[ts.pos + 1]
        return token
    catch e
        e isa BoundsError && return nothing
        rethrow()
    end
end

"""
    context(ts::TokenStream; before=2, after=5)

Strigify the current state of a TokenStream. Mainly
for debugging and errors.
"""
function context(ts::TokenStream; pos=position(ts), before=3, after=5)
    str = ""
    for i in (pos-before):(pos-1)
        if i == pos
            str *= " 🔥"
        end
        try
            str *= " " * string(ts[i])
        catch err
            err isa BoundsError || rethrow()
        end
    end
    return str
end

"""
    ParsingError <: Exception
"""
struct ParsingError <: Exception
    msg::String
    stream::TokenStream
    streampos::Int
end

ParsingError(msg::String, stream::TokenStream) = ParsingError(msg, stream, position(stream)+1)

function Base.show(io::IO, err::ParsingError)
    print(io, "ParsingError: ", err.msg)
    t = err.stream[err.streampos]
    t===nothing || print(io, " around line #", t.line)
    print(io, "\n", context(err.stream, pos=err.streampos))
end


text:
using NNFEM

NT = 100
T = 1.0
Δt = T/NT

m, n =  15, 15
h = 1/m
# Create a very simple mesh
elements = []
prop = Dict("name"=> "PlaneStrain", "rho"=> 1.0, "E"=> 2.0, "nu"=> 0.35)
coords = zeros((m+1)*(n+1), 2)
for j = 1:n
    for i = 1:m
        idx = (m+1)*(j-1)+i
        elnodes = [idx; idx+1; idx+1+m+1; idx+m+1]
        ngp = 3
        nodes = [
            (i-1)*h (j-1)*h
            i*h (j-1)*h
            i*h j*h
            (i-1)*h j*h
        ]
        coords[elnodes, :] = nodes
        push!(elements, SmallStrainContinuum(nodes, elnodes, prop, ngp))
    end
end

# fixed on the bottom, push on the right
EBC = zeros(Int64, (m+1)*(n+1), 2)
FBC = zeros(Int64, (m+1)*(n+1), 2)
g = zeros((m+1)*(n+1), 2)
f = zeros((m+1)*(n+1), 2)

for i = 2:m
    idx = n*(m+1)+i
    EBC[idx,:] .= -2
    EBC[i,:] .= -2
end

for j = 1:n+1
    idx = (j-1)*(m+1)+1
    EBC[idx,:] .= -2
    idx = (j-1)*(m+1)+m+1
    EBC[idx,:] .= -2
end

ndims = 2
domain = Domain(coords, elements, ndims, EBC, g, FBC, f)

# """
#     init_mesh(domain::Domain)

# Prepares `domain` for use in custom operators.
# """
# function init_mesh(domain::Domain)
#     ccall((:init_mesh, "./build/libdata"), Cvoid, ())

#     for (iele,e) in enumerate(domain.elements)

#         el_eqns = getEqns(domain, iele)
#         el_eqns_active = findall(el_eqns .>= 1)
#         ccall((:create_mesh, "./build/libdata"), Cvoid,
#         (
#             Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint,
#             Ptr{Cdouble}, Cint
#         ), e.elnodes[:], e.coords[:], vcat([x[:] for x in e.dhdx]...), e.weights, hcat([x[:] for x in e.hs]...), length(e.elnodes), length(e.weights),
#         el_eqns_active, length(el_eqns_active))
#     end

#     ccall((:create_domain, "./build/libdata"), Cvoid,
#         (
#             Ptr{Cdouble}, Cint, Cint, Cint
#         ), domain.nodes[:], domain.nnodes, domain.neqs, domain.neles)
#     ccall((:print_mesh, "./build/libdata"), Cvoid, ())
# end

# init_mesh(domain)


text:
RecipesBase.@recipe function f(model::StateSpaceModel, forec::Forecast)
    if !isunivariate(model)
        error("This plot recipe currently works for univariate models only.")
    end
    dists = Vector{Distribution}(undef, length(forec.expected_value))
    for d in eachindex(dists)
        dists[d] = Normal(forec.expected_value[d][1], sqrt(forec.covariance[d][1]))
    end
    n = length(model.system.y)
    forec_idx = collect(n + 1: n + length(forec.expected_value))
    fillrange_color = :steelblue
    # Plot the series
    @series begin
        seriestype := :path
        seriescolor := "black"
        return model.system.y
    end
    # Plot the forecast
    @series begin
        seriescolor := "blue"
        return forec_idx, mean.(dists)
    end
    # Plot the Conf Interval
    @series begin
        seriestype := :path
        seriescolor := fillrange_color
        fillcolor := fillrange_color
        fillalpha := 0.5
        fillrange := quantile.(dists, 0.975)
        label := ""
        return forec_idx, quantile.(dists, 0.025)
    end
    @series begin
        seriestype := :path
        seriescolor := fillrange_color
        fillcolor := fillrange_color
        fillalpha := 0.5
        fillrange := quantile.(dists, 0.90)
        label := ""
        return forec_idx, quantile.(dists, 0.10)
    end
end


text:
import Base: eltype, Pair, Tuple, show, ==, iterate, length
import Graphs: AbstractEdge, SimpleEdge, src, dst, reverse

export AbstractMultipleEdge, MultipleEdge, mul

"""
    AbstractMultipleEdge{T, U} <: AbstractEdge{T}

An abstract type representing multiple edges.
"""
abstract type AbstractMultipleEdge{T, U} <: AbstractEdge{T} end

"""
    MultipleEdge{T, U} <: AbstractMultipleEdge{T, U}

A struct representing multiple edges.

## Examples
```jltestdoc
julia> using Graphs, Multigraphs

julia> me = MultipleEdge(1, 2, 3)
Multiple edge 1 => 2 with multiplicity 3

julia> for e in me println(e) end
Edge 1 => 2
Edge 1 => 2
Edge 1 => 2

```
"""
struct MultipleEdge{T<:Integer, U<:Integer} <: AbstractMultipleEdge{T, U}
    src::T
    dst::T
    mul::U
    function MultipleEdge(src::T, dst::T, mul::U) where {T<:Integer, U<:Integer}
        if mul > 0
            return new{T, U}(src, dst, mul)
        else
            error("a multiple edge should have positive multiplicity")
        end
    end
end

MultipleEdge(src, dst) = MultipleEdge(src, dst, one(Int))
function MultipleEdge(a::Vector{T}) where {T<:Integer}
    l = length(a)
    if l == 2
        return MultipleEdge(a[1], a[2])
    elseif l > 2
        return MultipleEdge(a[1], a[2], a[3])
    end
end
MultipleEdge(t::NTuple{3}) = MultipleEdge(t[1], t[2], t[3])
MultipleEdge(t::NTuple{2}) = MultipleEdge(t[1], t[2], one(Int))
MultipleEdge(p::Pair) = MultipleEdge(p.first, p.second, one(Int))
MultipleEdge(e::T) where {T<:AbstractEdge} = MultipleEdge(src(e), dst(e), one(Int))
eltype(e::T) where {T<:AbstractMultipleEdge} = eltype(src(e))

src(e::MultipleEdge) = e.src
dst(e::MultipleEdge) = e.dst

"""
    mul(e)

Return the multiplicity of the multiple edge `e`.

## Examples
```jltestdoc
julia> using Graphs, Multigraphs

julia> me = MultipleEdge(1, 2, 3)
Multiple edge 1 => 2 with multiplicity 3

julia> mul(me)
3

"""
mul(e::MultipleEdge) = e.mul

show(io::IO, e::AbstractMultipleEdge) = print(io, "Multiple edge $(src(e)) => $(dst(e)) with multiplicity $(mul(e))")

Tuple(e::AbstractMultipleEdge) = (src(e), dst(e), mul(e))
SimpleEdge(e::AbstractMultipleEdge) = SimpleEdge(src(e), dst(e))

reverse(e::T) where {T<:AbstractMultipleEdge} = MultipleEdge(dst(e), src(e), mul(e))
==(e1::AbstractMultipleEdge, e2::AbstractMultipleEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2) && mul(e1) == mul(e2))
==(e1::AbstractMultipleEdge, e2::AbstractEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2) && mul(e1) == 1)
==(e1::AbstractEdge, e2::AbstractMultipleEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2) && mul(e2) == 1)

function iterate(e::MultipleEdge{T, U}, state::U=one(U)) where {T<:Integer, U<:Integer}
    if state > mul(e)
        return nothing
    else
        state += one(U)
        return (SimpleEdge(e), state)
    end
end

length(me::MultipleEdge) = mul(me)


text:
@enum Messages DONE=1

function start_server()
    port, server = listenany(8000)

    @async begin
        while true
            sock = accept(server)
            @async process_client(sock)
        end
    end
    port, server
end

function process_client(sock)
    while isopen(sock)
        data = deserialize(sock)
        response = process_message(data...)
        serialize(sock, response)
    end
end

macro safe(ex)
    esc(quote
        try
            $ex
        catch err
            err
        end
    end)
end

function process_message(data)
    println(data)
    DONE
end

function remotecall_fetch(f::Function, client::TCPSocket,args...)
    serialize(client, (f, args...) )
    deserialize(client)
end
process_message(f::Function, args...) = @safe f(args...)

function remote_eval(client, mod::Module, ex::Expr)
    serialize(client, (mod, ex) )
    deserialize(client)
end
function process_message(mod::Module,ex::Expr)
    @safe eval(mod,ex)
end


text:
module AoTImageSample

using Images
using TestImages, Images, ImageView

function main()
    img = load("julia.jpeg")
    save("aot_result.png",img)
    img = testimage("mandrill")
    imshow(img)
end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    main()
    return 0
end

end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using DandelionWebSockets.UniformResourceIdentifiers

@testset "Requirement 3.2-7" begin
    @testset "Port is not given; Scheme is http; Port is 80" begin
        # Act
        uri = URI("http://host/")

        # Assert
        @test uri.port == 80
    end

    @testset "Port is empty; Scheme is http; Port is 80" begin
        # Act
        uri = URI("http://host:/")

        # Assert
        @test uri.port == 80
    end

    @testset "Port is not given; Scheme is ws; Port is 80" begin
        # Act
        uri = URI("ws://host/")

        # Assert
        @test uri.port == 80
    end

    @testset "Port is not given; Scheme is wss; Port is 443" begin
        # Act
        uri = URI("wss://host/")

        # Assert
        @test uri.port == 443
    end

    @testset "Port is empty; Scheme is wss; Port is 443" begin
        # Act
        uri = URI("wss://host:/")

        # Assert
        @test uri.port == 443
    end

    @testset "Good enough URI; Scheme is ws; Parse is correct" begin
        # Act
        uri = URI("ws://hostname:42/abs/path")

        # Assert
        @test uri.scheme == "ws"
        @test uri.host == "hostname"
        @test uri.port == 42
        @test uri.abs_path == "/abs/path"
    end

    @testset "Good enough URI; No abs_path; abs_path is /" begin
        # Act
        uri = URI("ws://hostname:42")

        # Assert
        @test uri.abs_path == "/"
    end

    @testset "Good enough URI: Scheme is wss; issecure is true" begin
        # Act
        uri = URI("wss://hostname:42/abs/path")

        # Assert
        @test uri.issecure
    end

    @testset "Good enough URI: Scheme is ws; issecure is false" begin
        # Act
        uri = URI("ws://hostname:42/abs/path")

        # Assert
        @test !uri.issecure
    end

    @testset "Good enough URI: Scheme is http; issecure is false" begin
        # Act
        uri = URI("http://hostname:42/abs/path")

        # Assert
        @test !uri.issecure
    end

    @testset "Good enough URI: Scheme is https; issecure is true" begin
        # Act
        uri = URI("https://hostname:42/abs/path")

        # Assert
        @test uri.issecure
    end

    @testset "Good enough URI; URL is http://hostname/abs/path; Host name is hostname" begin
        # Act
        uri = URI("http://hostname/abs/path")

        # Assert
        @test uri.host == "hostname"
        @test uri.abs_path == "/abs/path"
    end
end


text:
include(joinpath("implicit_trapezoid", "integrator.jl"))
include(joinpath("implicit_trapezoid", "body.jl"))
include(joinpath("implicit_trapezoid", "test.jl"))


text:
using Pkg

Pkg.update()

pkg"activate . "

cd(@__DIR__)
ROOT_PATH = pwd()
push!(LOAD_PATH, ROOT_PATH, "src")

pkg"instantiate"
pkg"precompile"


text:
"""
Represent any source (incident) wave

Subtypes may have a symmetry (such as [`PlaneSource`](@ref)) and will contain information about physical medium.
"""
abstract type AbstractSource{T,P} end

"""
    PlaneSource(medium::P, amplitude::SVector, direction::SVector)

Is a struct type which describes a plane-wave source that drives/forces the whole system. It has three fields: a physical `medium`, an `amplitude` of the source, and the direction the propagate in `direction`.

For any given angular frequency ω, the PlaneSource has the value ``e^{i ω/c \\mathbf v \\cdot \\mathbf x }`` at the point ``\\mathbf x``, where ``c`` is the medium wavespeed and ``\\mathbf v`` is the direction.
"""
struct PlaneSource{T<:AbstractFloat,Dim,FieldDim,P<:PhysicalMedium} <: AbstractSource{T,P}
    medium::P
    direction::SVector{Dim,T}
    position::SVector{Dim,T}
    amplitude::SVector{FieldDim,Complex{T}}
    # Check that P has same Dim and FieldDim
    function PlaneSource{T,Dim,FieldDim,P}(medium::P, direction::AbstractArray{T}, position::AbstractArray{T}, amplitude::AbstractArray{CT} where CT<:Union{T,Complex{T}}) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim}}
        normw = sqrt(sum(direction .^2)) # note if direction is complex this is different from norm(direction)
        if !(normw ≈ one(T))
            @warn "The direction will be normalised so that sum(direction .^2) == 1.0"
        end

        if length(direction) != Dim || length(position) != Dim
            throw(DimensionMismatch("The spatial dimensions of the medium do not match dimensions of the direction or position vector."))
        elseif length(amplitude) != FieldDim
            throw(DimensionMismatch("The dimensions of the field in the medium (i.e. always 1 for acoustics) do not match the dimension of the amplitude vector."))
        else
            new{T,Dim,FieldDim,P}(medium,direction ./ normw, position, amplitude)
        end
    end
end

field(s::PlaneSource{T}) where T = function (x::AbstractArray{T}, ω::T)
        s.amplitude * exp(im * (ω / s.medium.c) * dot(x - s.position,s.direction))
    end

field(s::PlaneSource{T,Dim,1}) where {T, Dim} = function (x::AbstractArray{T}, ω::T)
        s.amplitude[1] * exp(im * (ω / s.medium.c) * dot(x - s.position,s.direction))
    end

field(s::PlaneSource{T}, x::AbstractArray{T}, ω::T) where T = field(s)(x,ω)

# Convenience constructor which does not require explicit types/parameters
function PlaneSource(medium::PhysicalMedium{T,Dim,FieldDim},
        direction::AbstractArray{T},
        position::AbstractArray{T} = zeros(T,Dim),
        amplitude::Union{CT,AbstractArray{CT}} where CT<:Union{T,Complex{T}} = SVector{FieldDim}(ones(Complex{T},FieldDim))
    ) where {T,Dim,FieldDim}

    if !(typeof(amplitude) <: AbstractArray)
        amplitude = [amplitude]
    end

    PlaneSource{T,Dim,FieldDim,typeof(medium)}(medium,direction,position,amplitude)
end

function PlaneSource(medium::PhysicalMedium{T,Dim,FieldDim};
        direction::AbstractArray{T} = [one(T); zeros(T,Dim-1)],
        position::AbstractArray{T} = zeros(T,Dim),
        amplitude::Union{CT,AbstractArray{CT}} where CT<:Union{T,Complex{T}} = SVector{FieldDim}(ones(Complex{T},FieldDim))
    ) where {T,Dim,FieldDim}

    PlaneSource(medium,direction,position,amplitude)
end

"""
    Source(medium::P, field::Function, coef::Function)

Is a struct type which describes the source field that drives/forces the whole system. It is also described as an incident wave. It has three fields `Source.medium`, `Source.field`, and `Source.coef`.

The source field at the position 'x' and angular frequency 'ω' is given by
```julia-repl
x = [1.0,0.0]
ω = 1.0
Source.field(x,ω)
```

The field `Source.coef`
regular_basis_function(medium::Acoustic{T,2}, ω::T)
"""
struct Source{T<:AbstractFloat,P<:PhysicalMedium} <: AbstractSource{T,P}
    medium::P
    "Use: field(x,ω)"
    field::Function
    "Use: coefficients(n,x,ω)"
    coefficients::Function
    # Enforce that the Types are the same
    function Source{T,P}(medium::P,field::Function,coefficients::Function) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim}}
        s = new{T,P}(medium,field,coefficients)
        self_test(s)
        return s
    end
end

function Source(medium::P,field::Function,coefficients::Function) where {T,Dim,FieldDim,P<:PhysicalMedium{T,Dim,FieldDim}}
    return Source{T,P}(medium,field,coefficients)
end

"""
    regular_spherical_coefficients(source::Source)

return a function which can calculate the coefficients of a regular spherical wave basis.
"""
function regular_spherical_coefficients(source::Source)
    source.coefficients
end

"""
Check that the source functions return the correct types
"""
function self_test(source::Source{T,P}) where {P,T}

medium, source_field, source_coef

    ω = one(T)

    # choose rand postion, hopefully not the source position/origin
    x = SVector(rand(T.(0.1:0.1:1.0),spatial_dimension(P))...)

    # Check that the result of field has same dimension and type as PhysicalProperty field
    if field_dimension(P) == 1
        source.field(x,ω)::Union{Complex{T},SVector{field_dimension(P),Complex{T}}}
    else
        source.field(x,ω)::SVector{field_dimension(P),Complex{T}}
    end

    # Check that the result of field has the same dimension as field dimension of P
    if field_dimension(P) == 1
        source.coefficients(1,x,ω)::Union{Vector{Complex{T}},Vector{SVector{field_dimension(P),Complex{T}}}}
    # else # this else is not necessarily correct..
    #     source.coefficients(1,x,ω)::SVector{field_dimension(P),Complex{T}}
    end

    return true
end

field(s::Source) = s.field
field(s::Source{T},x::AbstractArray{T}, ω::T) where T = field(s)(x,ω)

function constant_source(medium::P, num::Complex{T} = zero(Float64) * im) where {P,T}
    return Source{T,P}(medium, (x,ω) -> num, (order,x,ω) -> [num])
end

"""
    regular_spherical_source(PhysicalMedium,regular_coefficients; amplitude = one(T), position = zeros(T,Dim))

``c_n = ```regular_coefficients[n]`, ``x_o=```position`, and let ``v_n(kx)`` represent the regular spherical basis with wavenumber ``k`` at position ``x``. The above returns a `source` where `source.field(x,ω) =```\\sum_n c_n v_n(kx -k x_0)``
"""
function regular_spherical_source(medium::PhysicalMedium{T,Dim},regular_coefficients::AbstractVector{CT}; amplitude::Union{T,Complex{T}} = one(T), position::AbstractArray{T} = SVector(zeros(T,Dim)...)) where {T,Dim,CT<:Union{T,Complex{T}}}

    coeff_order = basislength_to_basisorder(typeof(medium),length(regular_coefficients))

    function source_field(x,ω)
        vs = regular_basis_function(medium, ω)(coeff_order,x-position)
        return amplitude * sum(vs .* regular_coefficients)

    end

    function source_coef(order,centre,ω)
        k = ω / medium.c

        V = regular_translation_matrix(medium, max(order,coeff_order), ω, centre - position)

        len2 = basisorder_to_basislength(typeof(medium),order)
        len1 = length(regular_coefficients)

        return amplitude .* sum(regular_coefficients[n] .* V[n,1:len2] for n in 1:len1)
    end

    return Source(medium, source_field, source_coef)
end

"""
    source_expand(Source, centre; basis_order = 4)

Returns a function of `(x,ω)` which approximates the value of the source at `(x,ω)`. That is, the source is written in terms of a regular basis expansion centred at `centre`.
"""
function source_expand(source::AbstractSource{T}, centre::AbstractVector{T}; basis_order::Int = 4) where T

    # Convert to SVector for efficiency and consistency
    centre = SVector(centre...)

    return function (x::AbstractVector{T}, ω::T)
        vs = regular_basis_function(source.medium, ω)
        regular_coefficients = regular_spherical_coefficients(source)

        sum(regular_coefficients(basis_order,centre,ω) .* vs(basis_order, x - centre))
    end
end

import Base.(+)
function +(s1::Source{T,P},s2::Source{T,P})::Source{T,P} where {P,T}
    if s1.medium != s2.medium
        error("Can not add sources from different physical mediums.")
    end
    field(x,ω) = s1.field(x,ω) + s2.field(x,ω)
    coef(n,centre,ω) = s1.coefficients(n,centre,ω) + s2.coefficients(n,centre,ω)

    Source{T,P}(s1.medium,field,coef)
end

import Base.(*)
function *(a,s::Source{T,P})::Source{T,P} where {P,T}
    if typeof(one(Complex{T})*a) != Complex{T}
        error("Multiplying source by $a would cause source type to change, please explicitly cast $a to same type as Source ($T)")
    end

    field(x,ω) = a * s.field(x,ω)
    coef(n,centre,ω) = a .* s.coefficients(n,centre,ω)

    Source{T,P}(s.medium,field,coef)
end

*(s::Source,a) = *(a,s)


text:
function likelihoods!(L::AbstractMatrix, hmm::AbstractHMM{Univariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(L) == (T, K)
    @inbounds for i in OneTo(K), t in OneTo(T)
        L[t,i] = pdf(hmm.B[i], observations[t])
    end
end

function loglikelihoods!(LL::AbstractMatrix, hmm::AbstractHMM{Univariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, K)
    @inbounds for i in OneTo(K), t in OneTo(T)
        LL[t,i] = logpdf(hmm.B[i], observations[t])
    end
end

function likelihoods!(L::AbstractMatrix, hmm::AbstractHMM{Multivariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(L) == (T, K)
    @inbounds for i in OneTo(K), t in OneTo(T)
        L[t,i] = pdf(hmm.B[i], view(observations, t, :))
    end
end

function loglikelihoods!(LL::AbstractMatrix, hmm::AbstractHMM{Multivariate}, observations)
    T, K = size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, K)
    @inbounds for i in OneTo(K), t in OneTo(T)
        LL[t,i] = logpdf(hmm.B[i], view(observations, t, :))
    end
end


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

# timing

# time() in libc.jl

# high-resolution relative time, in nanoseconds

"""
    time_ns()

Get the time in nanoseconds. The time corresponding to 0 is undefined, and wraps every 5.8 years.
"""
time_ns() = ccall(:jl_hrtime, UInt64, ())

# This type must be kept in sync with the C struct in src/gc.h
struct GC_Num
    allocd      ::Int64 # GC internal
    deferred_alloc::Int64 # GC internal
    freed       ::Int64 # GC internal
    malloc      ::UInt64
    realloc     ::UInt64
    poolalloc   ::UInt64
    bigalloc    ::UInt64
    freecall    ::UInt64
    total_time  ::UInt64
    total_allocd::UInt64 # GC internal
    since_sweep ::UInt64 # GC internal
    collect     ::Csize_t # GC internal
    pause       ::Cint
    full_sweep  ::Cint
end

gc_num() = ccall(:jl_gc_num, GC_Num, ())

# This type is to represent differences in the counters, so fields may be negative
struct GC_Diff
    allocd      ::Int64 # Bytes allocated
    malloc      ::Int64 # Number of GC aware malloc()
    realloc     ::Int64 # Number of GC aware realloc()
    poolalloc   ::Int64 # Number of pool allocation
    bigalloc    ::Int64 # Number of big (non-pool) allocation
    freecall    ::Int64 # Number of GC aware free()
    total_time  ::Int64 # Time spent in garbage collection
    pause       ::Int64 # Number of GC pauses
    full_sweep  ::Int64 # Number of GC full collection
end

gc_total_bytes(gc_num::GC_Num) =
    (gc_num.allocd + gc_num.deferred_alloc +
     Int64(gc_num.collect) + Int64(gc_num.total_allocd))

function GC_Diff(new::GC_Num, old::GC_Num)
    # logic from `src/gc.c:jl_gc_total_bytes`
    old_allocd = gc_total_bytes(old)
    new_allocd = gc_total_bytes(new)
    return GC_Diff(new_allocd - old_allocd,
                   Int64(new.malloc       - old.malloc),
                   Int64(new.realloc      - old.realloc),
                   Int64(new.poolalloc    - old.poolalloc),
                   Int64(new.bigalloc     - old.bigalloc),
                   Int64(new.freecall     - old.freecall),
                   Int64(new.total_time   - old.total_time),
                   new.pause              - old.pause,
                   new.full_sweep         - old.full_sweep)
end

function gc_alloc_count(diff::GC_Diff)
    diff.malloc + diff.realloc + diff.poolalloc + diff.bigalloc
end

# total time spend in garbage collection, in nanoseconds
gc_time_ns() = ccall(:jl_gc_total_hrtime, UInt64, ())

# total number of bytes allocated so far
gc_bytes() = ccall(:jl_gc_total_bytes, Int64, ())

# print elapsed time, return expression value
const _mem_units = ["byte", "KiB", "MiB", "GiB", "TiB", "PiB"]
const _cnt_units = ["", " k", " M", " G", " T", " P"]
function prettyprint_getunits(value, numunits, factor)
    if value == 0 || value == 1
        return (value, 1)
    end
    unit = ceil(Int, log(value) / log(factor))
    unit = min(numunits, unit)
    number = value/factor^(unit-1)
    return number, unit
end

function padded_nonzero_print(value,str)
    if value != 0
        blanks = "                "[1:18-length(str)]
        println("$str:$blanks$value")
    end
end

function format_bytes(bytes)
    bytes, mb = prettyprint_getunits(bytes, length(_mem_units), Int64(1024))
    if mb == 1
        Printf.@sprintf("%d %s%s", bytes, _mem_units[mb], bytes==1 ? "" : "s")
    else
        Printf.@sprintf("%.3f %s", bytes, _mem_units[mb])
    end
end

function time_print(elapsedtime, bytes, gctime, allocs)
    Printf.@printf("%10.6f seconds", elapsedtime/1e9)
    if bytes != 0 || allocs != 0
        allocs, ma = prettyprint_getunits(allocs, length(_cnt_units), Int64(1000))
        if ma == 1
            Printf.@printf(" (%d%s allocation%s: ", allocs, _cnt_units[ma], allocs==1 ? "" : "s")
        else
            Printf.@printf(" (%.2f%s allocations: ", allocs, _cnt_units[ma])
        end
        print(format_bytes(bytes))
        if gctime > 0
            Printf.@printf(", %.2f%% gc time", 100*gctime/elapsedtime)
        end
        print(")")
    elseif gctime > 0
        Printf.@printf(", %.2f%% gc time", 100*gctime/elapsedtime)
    end
    println()
end

function timev_print(elapsedtime, diff::GC_Diff)
    allocs = gc_alloc_count(diff)
    time_print(elapsedtime, diff.allocd, diff.total_time, allocs)
    print("elapsed time (ns): $elapsedtime\n")
    padded_nonzero_print(diff.total_time,   "gc time (ns)")
    padded_nonzero_print(diff.allocd,       "bytes allocated")
    padded_nonzero_print(diff.poolalloc,    "pool allocs")
    padded_nonzero_print(diff.bigalloc,     "non-pool GC allocs")
    padded_nonzero_print(diff.malloc,       "malloc() calls")
    padded_nonzero_print(diff.realloc,      "realloc() calls")
    padded_nonzero_print(diff.freecall,     "free() calls")
    padded_nonzero_print(diff.pause,        "GC pauses")
    padded_nonzero_print(diff.full_sweep,   "full collections")
end

"""
    @time

A macro to execute an expression, printing the time it took to execute, the number of
allocations, and the total number of bytes its execution caused to be allocated, before
returning the value of the expression.

See also [`@timev`](@ref), [`@timed`](@ref), [`@elapsed`](@ref), and
[`@allocated`](@ref).

```julia-repl
julia> @time rand(10^6);
  0.001525 seconds (7 allocations: 7.630 MiB)

julia> @time begin
           sleep(0.3)
           1+1
       end
  0.301395 seconds (8 allocations: 336 bytes)
2
```
"""
macro time(ex)
    quote
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        local diff = GC_Diff(gc_num(), stats)
        time_print(elapsedtime, diff.allocd, diff.total_time,
                   gc_alloc_count(diff))
        val
    end
end

"""
    @timev

This is a verbose version of the `@time` macro. It first prints the same information as
`@time`, then any non-zero memory allocation counters, and then returns the value of the
expression.

See also [`@time`](@ref), [`@timed`](@ref), [`@elapsed`](@ref), and
[`@allocated`](@ref).

```julia-repl
julia> @timev rand(10^6);
  0.001006 seconds (7 allocations: 7.630 MiB)
elapsed time (ns): 1005567
bytes allocated:   8000256
pool allocs:       6
malloc() calls:    1
```
"""
macro timev(ex)
    quote
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        timev_print(elapsedtime, GC_Diff(gc_num(), stats))
        val
    end
end

"""
    @elapsed

A macro to evaluate an expression, discarding the resulting value, instead returning the
number of seconds it took to execute as a floating-point number.

See also [`@time`](@ref), [`@timev`](@ref), [`@timed`](@ref),
and [`@allocated`](@ref).

```julia-repl
julia> @elapsed sleep(0.3)
0.301391426
```
"""
macro elapsed(ex)
    quote
        local t0 = time_ns()
        local val = $(esc(ex))
        (time_ns()-t0)/1e9
    end
end

# measure bytes allocated without *most* contamination from compilation
# Note: This reports a different value from the @time macros, because
# it wraps the call in a function, however, this means that things
# like:  @allocated y = foo()
# will not work correctly, because it will set y in the context of
# the local function made by the macro, not the current function
"""
    @allocated

A macro to evaluate an expression, discarding the resulting value, instead returning the
total number of bytes allocated during evaluation of the expression. Note: the expression is
evaluated inside a local function, instead of the current context, in order to eliminate the
effects of compilation, however, there still may be some allocations due to JIT compilation.
This also makes the results inconsistent with the `@time` macros, which do not try to adjust
for the effects of compilation.

See also [`@time`](@ref), [`@timev`](@ref), [`@timed`](@ref),
and [`@elapsed`](@ref).

```julia-repl
julia> @allocated rand(10^6)
8000080
```
"""
macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = gc_bytes()
                $(esc(ex))
                gc_bytes() - b0
            end
            f()
        end
    end
end

"""
    @timed

A macro to execute an expression, and return the value of the expression, elapsed time,
total bytes allocated, garbage collection time, and an object with various memory allocation
counters.

See also [`@time`](@ref), [`@timev`](@ref), [`@elapsed`](@ref), and
[`@allocated`](@ref).

```julia-repl
julia> val, t, bytes, gctime, memallocs = @timed rand(10^6);

julia> t
0.006634834

julia> bytes
8000256

julia> gctime
0.0055765

julia> fieldnames(typeof(memallocs))
9-element Array{Symbol,1}:
 :allocd
 :malloc
 :realloc
 :poolalloc
 :bigalloc
 :freecall
 :total_time
 :pause
 :full_sweep

julia> memallocs.total_time
5576500
```
"""
macro timed(ex)
    quote
        local stats = gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = time_ns() - elapsedtime
        local diff = GC_Diff(gc_num(), stats)
        val, elapsedtime/1e9, diff.allocd, diff.total_time/1e9, diff
    end
end

## printing with color ##

function with_output_color(f::Function, color::Union{Int, Symbol}, io::IO, args...; bold::Bool = false)
    buf = IOBuffer()
    iscolor = get(io, :color, false)
    try f(IOContext(buf, io), args...)
    finally
        str = String(take!(buf))
        if !iscolor
            print(io, str)
        else
            bold && color == :bold && (color = :nothing)
            enable_ansi  = get(text_colors, color, text_colors[:default]) *
                               (bold ? text_colors[:bold] : "")
            disable_ansi = (bold ? disable_text_style[:bold] : "") *
                               get(disable_text_style, color, text_colors[:default])
            first = true
            for line in split(str, '\n')
                first || print(buf, '\n')
                first = false
                isempty(line) && continue
                print(buf, enable_ansi, line, disable_ansi)
            end
            print(io, String(take!(buf)))
        end
    end
end

"""
    printstyled([io], xs...; bold::Bool=false, color::Union{Symbol,Int}=:normal)

Print `xs` in a color specified as a symbol or integer, optionally in bold.

`color` may take any of the values $(Base.available_text_colors_docstring)
or an integer between 0 and 255 inclusive. Note that not all terminals support 256 colors.
If the keyword `bold` is given as `true`, the result will be printed in bold.
"""
printstyled(io::IO, msg...; bold::Bool=false, color::Union{Int,Symbol}=:normal) =
    with_output_color(print, color, io, msg...; bold=bold)
printstyled(msg...; bold::Bool=false, color::Union{Int,Symbol}=:normal) =
    printstyled(STDOUT, msg...; bold=bold, color=color)

function julia_cmd(julia=joinpath(Sys.BINDIR, julia_exename()))
    opts = JLOptions()
    cpu_target = unsafe_string(opts.cpu_target)
    image_file = unsafe_string(opts.image_file)
    compile = if opts.compile_enabled == 0
                  "no"
              elseif opts.compile_enabled == 2
                  "all"
              elseif opts.compile_enabled == 3
                  "min"
              else
                  "yes"
              end
    depwarn = if opts.depwarn == 0
                  "no"
              elseif opts.depwarn == 2
                  "error"
              else
                  "yes"
              end
    `$julia -C$cpu_target -J$image_file --compile=$compile --depwarn=$depwarn`
end

function julia_exename()
    if ccall(:jl_is_debugbuild, Cint, ()) == 0
        return @static Sys.iswindows() ? "julia.exe" : "julia"
    else
        return @static Sys.iswindows() ? "julia-debug.exe" : "julia-debug"
    end
end

"""
    securezero!(o)

`securezero!` fills the memory associated with an object `o` with zeros.
Unlike `fill!(o,0)` and similar code, which might be optimized away by
the compiler for objects about to be discarded, the `securezero!` function
will always be called.
"""
function securezero! end
@noinline securezero!(a::AbstractArray{<:Number}) = fill!(a, 0)
securezero!(s::String) = unsafe_securezero!(pointer(s), sizeof(s))
@noinline unsafe_securezero!(p::Ptr{T}, len::Integer=1) where {T} =
    ccall(:memset, Ptr{T}, (Ptr{T}, Cint, Csize_t), p, 0, len*sizeof(T))
unsafe_securezero!(p::Ptr{Cvoid}, len::Integer=1) = Ptr{Cvoid}(unsafe_securezero!(Ptr{UInt8}(p), len))

if Sys.iswindows()
function getpass(prompt::AbstractString)
    print(prompt)
    flush(STDOUT)
    p = Vector{UInt8}(uninitialized, 128) # mimic Unix getpass in ignoring more than 128-char passwords
                          # (also avoids any potential memory copies arising from push!)
    try
        plen = 0
        while true
            c = ccall(:_getch, UInt8, ())
            if c == 0xff || c == UInt8('\n') || c == UInt8('\r')
                break # EOF or return
            elseif c == 0x00 || c == 0xe0
                ccall(:_getch, UInt8, ()) # ignore function/arrow keys
            elseif c == UInt8('\b') && plen > 0
                plen -= 1 # delete last character on backspace
            elseif !iscntrl(Char(c)) && plen < 128
                p[plen += 1] = c
            end
        end
        return unsafe_string(pointer(p), plen) # use unsafe_string rather than String(p[1:plen])
                                               # to be absolutely certain we never make an extra copy
    finally
        securezero!(p)
    end

    return ""
end
else
getpass(prompt::AbstractString) = unsafe_string(ccall(:getpass, Cstring, (Cstring,), prompt))
end

"""
    prompt(message; default="", password=false) -> Union{String, Nothing}

Displays the `message` then waits for user input. Input is terminated when a newline (\\n)
is encountered or EOF (^D) character is entered on a blank line. If a `default` is provided
then the user can enter just a newline character to select the `default`. Alternatively,
when the `password` keyword is `true` the characters entered by the user will not be
displayed.
"""
function prompt(message::AbstractString; default::AbstractString="", password::Bool=false)
    if Sys.iswindows() && password
        error("Command line prompt not supported for password entry on windows. Use `Base.winprompt` instead")
    end
    msg = !isempty(default) ? "$message [$default]:" : "$message:"
    if password
        # `getpass` automatically chomps. We cannot tell an EOF from a '\n'.
        uinput = getpass(msg)
    else
        print(msg)
        uinput = readline(keep=true)
        isempty(uinput) && return nothing  # Encountered an EOF
        uinput = chomp(uinput)
    end
    isempty(uinput) ? default : uinput
end

# Windows authentication prompt
if Sys.iswindows()
    struct CREDUI_INFO
        cbSize::UInt32
        parent::Ptr{Cvoid}
        pszMessageText::Ptr{UInt16}
        pszCaptionText::Ptr{UInt16}
        banner::Ptr{Cvoid}
    end

    const CREDUIWIN_GENERIC                 = 0x0001
    const CREDUIWIN_IN_CRED_ONLY            = 0x0020
    const CREDUIWIN_ENUMERATE_CURRENT_USER  = 0x0200

    const CRED_PACK_GENERIC_CREDENTIALS     = 0x0004

    const ERROR_SUCCESS                     = 0x0000
    const ERROR_CANCELLED                   = 0x04c7

    function winprompt(message, caption, default_username; prompt_username = true)
        # Step 1: Create an encrypted username/password bundle that will be used to set
        #         the default username (in theory could also provide a default password)
        credbuf = Vector{UInt8}(uninitialized, 1024)
        credbufsize = Ref{UInt32}(sizeof(credbuf))
        succeeded = ccall((:CredPackAuthenticationBufferW, "credui.dll"), stdcall, Bool,
            (UInt32, Cwstring, Cwstring, Ptr{UInt8}, Ptr{UInt32}),
             CRED_PACK_GENERIC_CREDENTIALS, default_username, "", credbuf, credbufsize)
        @assert succeeded

        # Step 2: Create the actual dialog
        #      2.1: Set up the window
        messageArr = Base.cwstring(message)
        captionArr = Base.cwstring(caption)
        pfSave = Ref{Bool}(false)
        cred = Ref{CREDUI_INFO}(CREDUI_INFO(sizeof(CREDUI_INFO), C_NULL, pointer(messageArr), pointer(captionArr), C_NULL))
        dwflags = CREDUIWIN_GENERIC | CREDUIWIN_ENUMERATE_CURRENT_USER
        if !prompt_username
            # Disable setting anything other than default_username
            dwflags |= CREDUIWIN_IN_CRED_ONLY
        end
        authPackage = Ref{Culong}(0)
        outbuf_data = Ref{Ptr{Cvoid}}(C_NULL)
        outbuf_size = Ref{Culong}(0)

        #      2.2: Do the actual request
        code = ccall((:CredUIPromptForWindowsCredentialsW, "credui.dll"), stdcall, UInt32, (Ptr{CREDUI_INFO}, UInt32, Ptr{Culong},
            Ptr{Cvoid}, Culong, Ptr{Ptr{Cvoid}}, Ptr{Culong}, Ptr{Bool}, UInt32), cred, 0, authPackage, credbuf, credbufsize[],
            outbuf_data, outbuf_size, pfSave, dwflags)

        #      2.3: If that failed for any reason other than the user canceling, error out.
        #           If the user canceled, just return nothing
        if code == ERROR_CANCELLED
            return nothing
        elseif code != ERROR_SUCCESS
            error(Base.Libc.FormatMessage(code))
        end

        # Step 3: Convert encrypted credentials back to plain text
        passbuf = Vector{UInt16}(uninitialized, 1024)
        passlen = Ref{UInt32}(length(passbuf))
        usernamebuf = Vector{UInt16}(uninitialized, 1024)
        usernamelen = Ref{UInt32}(length(usernamebuf))
        # Need valid buffers for domain, even though we don't care
        dummybuf = Vector{UInt16}(uninitialized, 1024)
        succeeded = ccall((:CredUnPackAuthenticationBufferW, "credui.dll"), Bool,
            (UInt32, Ptr{Cvoid}, UInt32, Ptr{UInt16}, Ptr{UInt32}, Ptr{UInt16}, Ptr{UInt32}, Ptr{UInt16}, Ptr{UInt32}),
            0, outbuf_data[], outbuf_size[], usernamebuf, usernamelen, dummybuf, Ref{UInt32}(1024), passbuf, passlen)
        if !succeeded
            error(Base.Libc.FormatMessage())
        end

        # Step 4: Free the encrypted buffer
        # ccall(:SecureZeroMemory, Ptr{Cvoid}, (Ptr{Cvoid}, Csize_t), outbuf_data[], outbuf_size[]) - not an actual function
        unsafe_securezero!(outbuf_data[], outbuf_size[])
        ccall((:CoTaskMemFree, "ole32.dll"), Cvoid, (Ptr{Cvoid},), outbuf_data[])

        # Done.
        passbuf_ = passbuf[1:passlen[]-1]
        result = (String(transcode(UInt8, usernamebuf[1:usernamelen[]-1])),
                  String(transcode(UInt8, passbuf_)))
        securezero!(passbuf_)
        securezero!(passbuf)

        return result
    end

end

unsafe_crc32c(a, n, crc) = ccall(:jl_crc32c, UInt32, (UInt32, Ptr{UInt8}, Csize_t), crc, a, n)

_crc32c(a::Union{Array{UInt8},FastContiguousSubArray{UInt8,N,<:Array{UInt8}} where N}, crc::UInt32=0x00000000) =
    unsafe_crc32c(a, length(a) % Csize_t, crc)

_crc32c(s::String, crc::UInt32=0x00000000) = unsafe_crc32c(s, sizeof(s) % Csize_t, crc)

function _crc32c(io::IO, nb::Integer, crc::UInt32=0x00000000)
    nb < 0 && throw(ArgumentError("number of bytes to checksum must be ≥ 0"))
    # use block size 24576=8192*3, since that is the threshold for
    # 3-way parallel SIMD code in the underlying jl_crc32c C function.
    buf = Vector{UInt8}(uninitialized, min(nb, 24576))
    while !eof(io) && nb > 24576
        n = readbytes!(io, buf)
        crc = unsafe_crc32c(buf, n, crc)
        nb -= n
    end
    return unsafe_crc32c(buf, readbytes!(io, buf, min(nb, length(buf))), crc)
end
_crc32c(io::IO, crc::UInt32=0x00000000) = _crc32c(io, typemax(Int64), crc)
_crc32c(io::IOStream, crc::UInt32=0x00000000) = _crc32c(io, filesize(io)-position(io), crc)

"""
    @kwdef typedef

This is a helper macro that automatically defines a keyword-based constructor for the type
declared in the expression `typedef`, which must be a `struct` or `mutable struct`
expression. The default argument is supplied by declaring fields of the form `field::T =
default`. If no default is provided then the default is provided by the `kwdef_val(T)`
function.

# Examples
```jldoctest
julia> struct Bar end

julia> Base.@kwdef struct Foo
           a::Cint            # implied default Cint(0)
           b::Cint = 1        # specified default
           z::Cstring         # implied default Cstring(C_NULL)
           y::Bar             # implied default Bar()
       end
Foo

julia> Foo()
Foo(0, 1, Cstring(0x0000000000000000), Bar())
```
"""
macro kwdef(expr)
    expr = macroexpand(__module__, expr) # to expand @static
    T = expr.args[2]
    params_ex = Expr(:parameters)
    call_ex = Expr(:call, T)
    _kwdef!(expr.args[3], params_ex, call_ex)
    quote
        Base.@__doc__($(esc(expr)))
        $(esc(Expr(:call,T,params_ex))) = $(esc(call_ex))
    end
end

# @kwdef helper function
# mutates arguments inplace
function _kwdef!(blk, params_ex, call_ex)
    for i in eachindex(blk.args)
        ei = blk.args[i]
        isa(ei, Expr) || continue
        if ei.head == :(=)
            # var::Typ = defexpr
            dec = ei.args[1]  # var::Typ
            var = dec.args[1] # var
            def = ei.args[2]  # defexpr
            push!(params_ex.args, Expr(:kw, var, def))
            push!(call_ex.args, var)
            blk.args[i] = dec
        elseif ei.head == :(::)
            dec = ei # var::Typ
            var = dec.args[1] # var
            def = :(Base.kwdef_val($(ei.args[2])))
            push!(params_ex.args, Expr(:kw, var, def))
            push!(call_ex.args, dec.args[1])
        elseif ei.head == :block
            # can arise with use of @static inside type decl
            _kwdef!(ei, params_ex, call_ex)
        end
    end
    blk
end

"""
    kwdef_val(T)

The default value for a type for use with the `@kwdef` macro. Returns:

 - null pointer for pointer types (`Ptr{T}`, `Cstring`, `Cwstring`)
 - zero for integer types
 - no-argument constructor calls (e.g. `T()`) for all other types

# Examples
```jldoctest
julia> struct Foo
           i::Int
       end

julia> Base.kwdef_val(::Type{Foo}) = Foo(42)

julia> Base.@kwdef struct Bar
           y::Foo
       end
Bar

julia> Bar()
Bar(Foo(42))
```
"""
function kwdef_val end

kwdef_val(::Type{Ptr{T}}) where {T} = Ptr{T}(C_NULL)
kwdef_val(::Type{Cstring}) = Cstring(C_NULL)
kwdef_val(::Type{Cwstring}) = Cwstring(C_NULL)

kwdef_val(::Type{T}) where {T<:Integer} = zero(T)

kwdef_val(::Type{T}) where {T} = T()

# testing

"""
    Base.runtests(tests=["all"]; ncores=ceil(Int, Sys.CPU_CORES / 2),
                  exit_on_error=false, [seed])

Run the Julia unit tests listed in `tests`, which can be either a string or an array of
strings, using `ncores` processors. If `exit_on_error` is `false`, when one test
fails, all remaining tests in other files will still be run; they are otherwise discarded,
when `exit_on_error == true`.
If a seed is provided via the keyword argument, it is used to seed the
global RNG in the context where the tests are run; otherwise the seed is chosen randomly.
"""
function runtests(tests = ["all"]; ncores = ceil(Int, Sys.CPU_CORES / 2),
                  exit_on_error=false,
                  seed::Union{BitInteger,Nothing}=nothing)
    if isa(tests,AbstractString)
        tests = split(tests)
    end
    exit_on_error && push!(tests, "--exit-on-error")
    seed != nothing && push!(tests, "--seed=0x$(hex(seed % UInt128))") # cast to UInt128 to avoid a minus sign
    ENV2 = copy(ENV)
    ENV2["JULIA_CPU_CORES"] = "$ncores"
    try
        run(setenv(`$(julia_cmd()) $(joinpath(Sys.BINDIR,
            Base.DATAROOTDIR, "julia", "test", "runtests.jl")) $tests`, ENV2))
    catch
        buf = PipeBuffer()
        versioninfo(buf)
        error("A test has failed. Please submit a bug report (https://github.com/JuliaLang/julia/issues)\n" *
              "including error messages above and the output of versioninfo():\n$(read(buf, String))")
    end
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
module ErrorAnalysis
  using ...Statistics
  using Distances

  function generateConfusionMatrix(test_set,predicted_set)
    maxval = convert(Int,max(maximum(test_set),maximum(predicted_set)))
    minval = 1

    if(maxval==1)
      minval=0

      confusionmatrix = zeros(Int,maxval+1,maxval+1)

      for i=1:length(test_set)
        #println(test_set[i]," ",predicted_set[i])
        confusionmatrix[test_set[i]+1, predicted_set[i]+1]+=1
      end
    else
      confusionmatrix = zeros(Int,maxval,maxval)

      for i=1:length(test_set)
        #println(test_set[i]," ",predicted_set[i])
        confusionmatrix[test_set[i], predicted_set[i]]+=1
      end
    end
    topRow = collect(minval:maxval)
    sideRow = [0; topRow]

    return [sideRow [topRow';confusionmatrix]]

  end

  function calculateErrorRate(test_set,predicted_set)
    incorrect=0
    for i=1:length(test_set)
      if test_set[i]!=predicted_set[i]
        incorrect+=1
      end
    end
    println(incorrect/length(test_set))
    return incorrect/length(test_set)
  end

  function calculateWeightedSuccessRate(confusionmatrix)
    confusionmatrix = confusionmatrix[2:end,:]
    confusionmatrix = confusionmatrix[:,2:end]

    totals = Statistics.sumColumn(confusionmatrix)
    correct = diag(confusionmatrix)'

    averages = filter(x->!isnan(x),correct.*(1.0./totals))
    return mean(averages)
  end

  function calculateRootMeanStandardError(observed,actual)
    y = actual-observed;
    return sqrt(y'*y/length(y))[1]
  end

  function calculateWCSS(centroids,dataset, y)
    wcss = 0;
    for i=1:size(centroids)[1]
      cent = centroids[i,:]
      newdataset =dataset[find(x->x==i,y),:]
      for j=1:size(newdataset)[1]
        wcss+=(euclidean(cent,newdataset[j,:]))^2
      end
    end
    return wcss
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using InvertibleTransformations
using Test

@testset "InvertibleTransformations.jl" begin
    # Write your tests here.
end


text:
"""
    BufferedFiles

Implementations of a better Filehandles.

Includes a Read-Only Filehandle and a better Write-Only Filehandle

Uses Julias Filesystem.File for Portability.

Main Functions: [`open`](@ref) and [`@om_str`](@ref)
"""
module BufferedFiles
using Base.Filesystem

export @om_str, bufferedopen

function cmemcpy(dest, src, nbytes)
    @ccall memcpy(dest::Ptr{Cvoid}, src::Ptr{Cvoid}, nbytes::Csize_t)::Ptr{Cvoid}
    nothing
end

function cmemmove(dest, src, nbytes)
    @ccall memmove(dest::Ptr{Cvoid}, src::Ptr{Cvoid}, nbytes::Csize_t)::Ptr{Cvoid}
    nothing
end

"""
    const DEFAULT_BUFSIZE = 16 * (1024) * (1024)

Default size of the read/write-buffer of [`BufferedReadFile`](@ref) and [`BufferedWriteFile`](@ref).
"""
const DEFAULT_BUFSIZE = 16 * (1024) * (1024)

"""
    abstract type BufferedOpenMode end

Types Implemented by [`@om_str`](@ref), used in [`open`](@ref)
"""
abstract type BufferedOpenMode end

struct Read <: BufferedOpenMode end

struct WriteTrunc <: BufferedOpenMode end

struct Append <: BufferedOpenMode end

struct ReadWrite <: BufferedOpenMode end

struct ReadWriteTrunc <: BufferedOpenMode end

struct ReadAppend <: BufferedOpenMode end

struct MmapRead <: BufferedOpenMode end

struct MmapReadWrite <: BufferedOpenMode end

"""
    BufferedOpenMode(pattern::AbstractString)

Convert `pattern` into its corresponding subtype of `BufferedOpenMode`.

Recommended Usage of this is with the [`@om_str`](@ref) macro at parse-time.
"""
function BufferedOpenMode(pattern::AbstractString)
    if pattern == "r"
        return Read()
    elseif pattern == "w"
        return WriteTrunc()
    elseif pattern == "a"
        return Append()
    elseif pattern == "r+"
        return ReadWrite()
    elseif pattern == "w+"
        return ReadWriteTrunc()
    elseif pattern == "a+"
        return ReadAppend()
    elseif pattern == "mr"
        return MmapRead()
    elseif pattern == "mr+"
        return MmapReadWrite()
    end
    return error("Error, no matching string found!")
end

"""
    @om_str "pattern"
    om"pattern"

Convert `pattern` into its corresponding subtype of [`BufferedOpenMode`](@ref) at parse-time.
"""
macro om_str(pattern)
    return BufferedOpenMode(pattern)
end

"""
    abstract type BufferedFile <: IO end

Abstract Type for my Filehandles, to overwrite some of Julias functions from io.jl that cause allocations.
"""
abstract type BufferedFile <: IO end

const bufferedopen = Base.open

# Overwrites two @noinline julia functions which cause allocations in a lot of situations
function Base.unsafe_read(s::BufferedFile, p::Ref{T}, n::Integer) where {T}
    GC.@preserve p unsafe_read(s, Base.unsafe_convert(Ref{T}, p)::Ptr, n)
end
function Base.unsafe_write(s::BufferedFile, p::Ref{T}, n::Integer) where {T}
    GC.@preserve p unsafe_write(s, Base.unsafe_convert(Ref{T}, p)::Ptr, n)
end

function Base.unsafe_read(s::BufferedFile, p::Ptr, n::Integer)
    return unsafe_read(s, convert(Ptr{UInt8}, p), convert(UInt, n))
end
function Base.unsafe_write(s::BufferedFile, p::Ptr, n::Integer)
    return unsafe_write(s, convert(Ptr{UInt8}, p), convert(UInt, n))
end

"""
    BufferedReadFile

Buffered Filehandle for more performant Reading.

Opened via [`open(filename, om"r")`](@ref Base.open(::AbstractString, ::Read)).
"""
mutable struct BufferedReadFile <: BufferedFile
    rf::File
    buf::Vector{UInt8}
    filepos::Int
    fileend::Int
    bufpos::Int
    buflastread::Int
    function BufferedReadFile(rf, buf)
        f = new(rf, buf, position(rf), filesize(rf), 0, 0)
        finalizer(close, f)
        refresh(f)
        return f
    end
end

Base.isopen(p::BufferedReadFile) = isopen(p.rf)
Base.isreadable(p::BufferedReadFile) = isopen(p)
Base.iswritable(::BufferedReadFile) = false

"""
    open(file::AbstractString, om"r"[, bufsize::Integer])

Open a File in Read-Only Mode.

Returns a [`BufferedReadFile`](@ref).
"""
function Base.open(file::AbstractString, ::Read, bufsize::Integer=DEFAULT_BUFSIZE)
    rf = Filesystem.open(file, JL_O_RDONLY)
    buf = Vector{UInt8}(undef, bufsize)
    f = BufferedReadFile(rf, buf)
    return f
end

function refresh(f::BufferedReadFile)
    GC.@preserve f begin
        cmemmove(pointer(f.buf), pointer(f.buf) + f.bufpos, bytesavailable(f))
        f.buflastread = bytesavailable(f)
        f.bufpos = 0
        to_read::UInt = min(length(f.buf) - f.buflastread, max(0, f.fileend - f.filepos))
        f.filepos += to_read
        unsafe_read(f.rf, pointer(f.buf) + f.buflastread, to_read)
        f.buflastread += to_read
    end
    return f
end

function Base.position(f::BufferedReadFile)::Int
    return f.filepos - (f.buflastread - f.bufpos)
end

function Base.eof(f::BufferedReadFile)::Bool
    return (f.buflastread - f.bufpos) == 0
end

function Base.skip(f::BufferedReadFile, n::Integer)
    todo::UInt = bytesavailable(f)
    if n >= 0
        nb::UInt = n
        if todo == 0
            if nb > 0
                throw(EOFError())
            end
        elseif todo > nb
            f.bufpos += nb
        elseif todo <= nb
            f.bufpos = f.buflastread
            refresh(f)
            skip(f, (nb - todo)::UInt)
        end
    else
        u::UInt = abs(n)
        if f.bufpos < u
            dif = u - f.bufpos
            skip(f.rf, -Int(f.buflastread + dif))
            f.filepos -= (f.buflastread + dif)
            f.bufpos = 0
            f.buflastread = 0
            refresh(f)
        else
            f.bufpos -= u
        end
    end
    return f
end

function Base.seek(f::BufferedReadFile, n::Integer)
    if position(f) != n
        seek(f.rf, n)
        f.filepos = n
        f.buflastread = 0
        f.bufpos = 0
        refresh(f)
    end
    return f
end

function Base.close(f::BufferedReadFile)
    if isopen(f)
        close(f.rf)
        resize!(f.buf, 0)
        sizehint!(f.buf, 0)
    end
    return nothing
end

Base.bytesavailable(f::BufferedReadFile) = f.buflastread - f.bufpos

function Base.read(s::BufferedReadFile, ::Type{UInt8})
    r = bytesavailable(s)
    if r >= 2
        s.bufpos += 1
        return @inbounds s.buf[s.bufpos]
    elseif r == 1
        s.bufpos += 1
        cnt = @inbounds s.buf[s.bufpos]
        refresh(s)
        return cnt
    elseif r == 0
        throw(EOFError())
    else
        error("this should never happen!")
    end
end

function Base.unsafe_read(f::BufferedReadFile, p::Ptr{UInt8}, nb::UInt)
    todo::UInt = bytesavailable(f)
    if todo == 0
        if nb > 0
            throw(EOFError())
        else
            return nothing
        end
    end
    if todo > nb
        GC.@preserve f begin
            cmemcpy(p, pointer(f.buf) + f.bufpos, nb)
        end
        f.bufpos += nb
        return nothing
    elseif todo < nb
        GC.@preserve f begin
            cmemcpy(p, pointer(f.buf) + f.bufpos, todo)
        end
        f.bufpos = f.buflastread
        refresh(f)
        return unsafe_read(f, p + todo, nb - todo)
    elseif todo == nb
        GC.@preserve f begin
            cmemcpy(p, pointer(f.buf) + f.bufpos, todo)
        end
        f.bufpos = f.buflastread
        refresh(f)
        return nothing
    else
        error("This should never happen!")
    end
end

function Base.resize!(b::BufferedReadFile, ns::Integer)
    s = Int(ns)
    if b.buflastread < s
        resize!(b.buf, s)
    elseif b.buflastread >= s
        p = position(b)
        seek(b.rf, p)
        b.filepos = p
        b.buflastread = 0
        b.bufpos = 0
        resize!(b.buf, s)
        refresh(b)
    end
    return b
end

"""
    BufferedWriteFile

Buffered Filehandle for more performant sequential Writing.

Opened via [`open(filename, om"r")`](@ref open(::AbstractString, ::Read)).
"""
mutable struct BufferedWriteFile <: BufferedFile
    rf::File
    buf::Vector{UInt8}
    filepos::Int
    bufpos::Int
    function BufferedWriteFile(rf, buf)
        f = new(rf, buf, position(rf), 0)
        finalizer(close, f)
        return f
    end
end

function Base.resize!(b::BufferedWriteFile, ns::Integer)
    s = Int(ns)
    if s <= b.bufpos
        flush(b)
    end
    resize!(b.buf, s)
    return b
end

Base.isopen(p::BufferedWriteFile) = isopen(p.rf)
Base.isreadable(::BufferedWriteFile) = false
Base.iswritable(p::BufferedWriteFile) = isopen(p)

"""
    open(file::AbstractString, om"w"[, bufsize::Integer])

Open a File in Write-Only Mode.

Truncates the File when opening it.

Returns a [`BufferedWriteFile`](@ref).
"""
function Base.open(file::AbstractString,
                   ::WriteTrunc,
                   bufsize::Integer=DEFAULT_BUFSIZE)
    buf = Vector{UInt8}(undef, bufsize)
    rf = Filesystem.open(file, JL_O_CREAT | JL_O_TRUNC | JL_O_WRONLY,
                         S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH)
    return BufferedWriteFile(rf, buf)
end

function Base.flush(f::BufferedWriteFile)
    GC.@preserve f begin
        unsafe_write(f.rf, pointer(f.buf), f.bufpos)
    end
    f.filepos += f.bufpos
    f.bufpos = 0
    return f
end

function Base.position(f::BufferedWriteFile)::Int
    return f.filepos + f.bufpos
end

function Base.skip(f::BufferedWriteFile, n::Integer)
    flush(f)
    skip(f.rf, n)
    f.filepos += n
    return f
end

function Base.seek(f::BufferedWriteFile, n::Integer)
    if position(f) != n
        flush(f)
        seek(f.rf, n)
    end
    f.filepos = n
    return f
end

function Base.close(f::BufferedWriteFile)
    if isopen(f)
        flush(f)
        close(f.rf)
        resize!(f.buf, 0)
        sizehint!(f.buf, 0)
    end
    return nothing
end

function Base.unsafe_write(f::BufferedWriteFile, p::Ptr{UInt8}, nb::UInt)::UInt
    freespace = (length(f.buf) - f.bufpos)
    @assert freespace > 0
    if freespace > nb
        GC.@preserve f begin
            cmemcpy(pointer(f.buf) + f.bufpos, p, nb)
        end
        f.bufpos += nb
        return nb
    elseif freespace <= nb
        GC.@preserve f begin
            cmemcpy(pointer(f.buf) + f.bufpos, p, freespace)
        end
        f.bufpos += freespace
        flush(f)
        return freespace + unsafe_write(f, p + freespace, nb - freespace)
    else
        error("This should never happen!")
    end
end

function Base.write(s::BufferedWriteFile, u::UInt8)
    freespace = (length(s.buf) - s.bufpos)
    if freespace >= 2
        s.bufpos += 1
        @inbounds s.buf[s.bufpos] = u
    elseif freespace == 1
        s.bufpos += 1
        @inbounds s.buf[s.bufpos] = u
        flush(s)
    else
        error("This should never happen!")
    end
    return 1
end

end # module


text:
import Base: read, write, close, eof

# TODO: Documentation

"An exception thrown into a task in order to stop it."
mutable struct StopTaskException <: Exception end

abstract type AbstractServerReader end

"Reading from a network socket and placing the resulting frame on a channel."
struct ServerReader <: AbstractServerReader
    s::IO
    task::Task
end

"Read frames from the network, until an exception is thrown in this task."
function do_reader(s::IO, logic::AbstractClientProtocol)
    # Requirement
    # @6_2-1 Incoming data
    # @7_1_1-2 Discard trailing bytes

    try
        while true
            frame = read(s, Frame)
            # This is a proxy, so the actual `handle` call made on the logic object is done in a
            # separate coroutine.
            handle(logic, FrameFromServer(frame))
        end
    catch ex
        # TODO: Handle errors better.
    end
    try
        handle(logic, SocketClosed())
    catch ex
        # TODO: Log this error.
    end
end

function start_reader(s::IO, logic::AbstractClientProtocol)
    t = @async do_reader(s, logic)
    ServerReader(s, t)
end

function stop(t::ServerReader)
    try
        Base.throwto(t.task, StopTaskException())
    catch
        # Ignore any stop errors, at least for the time being.
    end
end

#
#
#

"""
TLSBufferedIO adapts a TLS socket so we can do byte I/O.

The stream returned by MbedTLS when using a TLS socket does not support the byte I/O used when
reading a frame. It only supports reading a chunk of data. This is a fake stream that buffers some
data and lets us do byte I/O.

Note: This should have been done by the BufferedStreams.jl package. However, I couldn't get it to
work with the MbedTLS stream, for reasons unknown. If we can investigate and fix that problem, then
we should really replace this type with a BufferedInputStream.
"""
struct TLSBufferedIO <: IO
    tls_stream::IO
    buf::IOBuffer

    TLSBufferedIO(tls_stream::IO) = new(tls_stream, IOBuffer())
end

"Read all available data, and block until we have enough to fulfíll the next read."
function fill_buffer(s::TLSBufferedIO, n::Int)
    begin_ptr = mark(s.buf)
    while s.buf.size - begin_ptr < n
        if eof(s.tls_stream)
            throw(EOFError())
        end
        write(s.buf, readavailable(s.tls_stream))
    end
    reset(s.buf)
end

function read(s::TLSBufferedIO, t::Type{UInt8})
    fill_buffer(s, sizeof(t))
    read(s.buf, t)
end

function read(s::TLSBufferedIO, t::Type{UInt16})
    fill_buffer(s, sizeof(t))
    read(s.buf, t)
end

function read(s::TLSBufferedIO, t::Type{UInt64})
    fill_buffer(s, sizeof(t))
    read(s.buf, t)
end

function read(s::TLSBufferedIO, t::Type{UInt8}, n::Int)
    fill_buffer(s, sizeof(t) * n)
    read(s.buf, t, n)
end

write(s::TLSBufferedIO, t::UInt8) = write(s.tls_stream, t)
write(s::TLSBufferedIO, t::UInt16) = write(s.tls_stream, t)
write(s::TLSBufferedIO, t::UInt64) = write(s.tls_stream, t)

close(s::TLSBufferedIO) = close(s.tls_stream)
eof(s::TLSBufferedIO) = eof(s.tls_stream)


text:
module TestingCartesianBoxes

using Test
using CartesianBoxes
using CartesianBoxes: indices, ranges

# Deal with compatibility issues.
import Base: CartesianIndices, axes
const CartesianRegions = Union{CartesianIndices, CartesianBox}

@inline _add(a,b) = a + b
@inline _sub(a,b) = a - b
@inline _mul(a,b) = a * b
@inline _div(a,b) = a / b

for op in (:_add, :_sub, :_mul, :_div)
    slow = Symbol(:slow,op,:!)
    fast = Symbol(:fast,op,:!)
    @eval begin
        function $slow(dst::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::AbstractArray{T,N}) where {T,N}
            inds = axes(dst)
            @assert axes(A) == inds
            @assert axes(B) == inds
            for i in eachindex(dst, A, B)
                dst[i] = $op(A[i], B[i])
            end
            return dst
        end

        function $slow(::Type{R},
                       dst::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::AbstractArray{T,N}) where {R<:CartesianRegions,T,N}
            inds = axes(dst)
            @assert axes(A) == inds
            @assert axes(B) == inds
            for i in R(inds)
                dst[i] = $op(A[i], B[i])
            end
            return dst
        end

        function $slow(R::CartesianRegions,
                       dst::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::AbstractArray{T,N}) where {T,N}
            inds = axes(dst)
            @assert axes(A) == inds
            @assert axes(B) == inds
            @assert isnonemptypartof(CartesianBox(R), A)
            for i in R
                dst[i] = $op(A[i], B[i])
            end
            return dst
        end

        # Implementation using what is provided by Julia base (probably linear
        # indexing for regular arrays).
        function $fast(dst::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::AbstractArray{T,N}) where {T,N}
            inds = axes(dst)
            @assert axes(A) == inds
            @assert axes(B) == inds
            @inbounds @simd for i in eachindex(dst, A, B)
                dst[i] = $op(A[i], B[i])
            end
            return dst
        end

        # Implementation using a specific iterator.
        function $fast(::Type{R},
                       dst::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::AbstractArray{T,N}) where {R<:CartesianRegions,T,N}
            inds = axes(dst)
            @assert axes(A) == inds
            @assert axes(B) == inds
            @inbounds @simd for i in R(inds)
                dst[i] = $op(A[i], B[i])
            end
            return dst
        end

        # Operation over a sub-region.
        function $fast(R::CartesianRegions,
                       dst::AbstractArray{T,N},
                       A::AbstractArray{T,N},
                       B::AbstractArray{T,N}) where {T,N}
            inds = axes(dst)
            @assert axes(A) == inds
            @assert axes(B) == inds
            @assert isnonemptypartof(CartesianBox(R), A)
            @inbounds @simd for i in R
                dst[i] = $op(A[i], B[i])
            end
            return dst
        end
    end
end

maxabsdif(A, B) = maximum(abs.(A .- B))

# This method is to exercise the loop and check the computed length.
function stupidcount(iter)
    n = 0
    for i in iter
        n += 1
    end
    return n
end

const SIZES = ((), (45,), (21,22), (5,6,7,8))
const TYPES = (Float32, )

@testset "CartesianBoxes" begin
    @testset "Basic operations" begin
        for inds in (VERSION < v"1.6.0-beta1"
                     ? ((Base.OneTo(5), 2:8),
                        (-1:7, 0:8, Base.OneTo(4)))
                     : ((Base.OneTo(5), 2:8),
                        (-1:2:7, 0:8, Base.OneTo(4))))
            B = CartesianBox(inds...)
            @test ndims(B) == length(inds)
            @test Tuple(B) === inds
            @test axes(B) === map(x -> Base.OneTo(length(x)), inds)
            @test size(B) == map(length, inds)
            @test length(B) == prod(map(length, inds))
            @test isempty(B) == (length(B) == 0)
        end
    end

    @testset "Cartesian boxes of arrays" for dims in SIZES, T in TYPES
        A = Array{T}(undef, dims)
        r = CartesianIndices(A)
        b = CartesianBox(A)
        @test ndims(b) == ndims(r)
        @test eltype(b) === eltype(r)
        @test length(b) == stupidcount(b)
        @test length(b) == length(r)
        @test size(b) == size(r)
        @test axes(b) == axes(r)
        @test indices(b) == indices(r)
        for d in 1:ndims(b)
            @test size(b,d) == size(r,d)
            @test axes(b,d) == axes(r,d)
        end
        @test first(b) === first(r)
        @test last(b) === last(r)
        buf = IOBuffer();
        show(buf, MIME"text/plain"(), b);
        @test String(take!(buf)) == repr(b)

        @test CartesianBox(r) == b
        @test CartesianIndices(b) == r
        @test CartesianIndices(CartesianBox(r)) === r
        @test CartesianBox(CartesianIndices(b)) === b

        @test indices(b) === Tuple(b)
        @test CartesianBox(indices(b)) === b
        @test CartesianBox(indices(b)...) === b
        @test CartesianBox(first(b), last(b)) == b
        flag = true
        for i in eachindex(r, b)
            flag &= (r[i] == b[i])
        end
        @test flag

        I0 = CartesianIndex(ntuple(i -> 0, ndims(b)))
        I1 = CartesianIndex(ntuple(i -> i, ndims(b)))
        @test indices(I1) === I1.I
        @test (b + I0) == b
        @test (b - I0) == b
        @test (b + I0.I) == b
        @test (b - I0.I) == b
        for f in (+, -)
            c = f(b,I1)
            @test f(b,I1) == CartesianBox(map((r,i) -> f(first(r),i):f(last(r),i),
                                              indices(b), indices(I1)))
            @test CartesianBox(indices(c)) === c
            @test CartesianBox(indices(c)...) === c
            @test CartesianBox(first(c), last(c)) == c
        end
    end

    @testset "Array indexing" for dims in SIZES, T in TYPES
        A = randn(T, dims)
        B = randn(T, dims)
        C1 = Array{T}(undef, dims)
        C2 = Array{T}(undef, dims)
        C3 = Array{T}(undef, dims)

        @testset "Operation $oper" for oper in (slow_add!, slow_sub!,
                                                slow_mul!)
            oper(C1, A, B)
            oper(CartesianBox, fill!(C2, 0), A, B)
            @test maxabsdif(C1, C2) == 0
            if ndims(A) ≥ 1
                # Take a sub-region.
                S = CartesianBox(map(n -> (min(2,n) : max(n-1,1)), size(A)))
                oper(CartesianBox(S),     fill!(C2, 0), A, B)
                oper(CartesianIndices(S), fill!(C3, 0), A, B)
                @test maxabsdif(C2, C3) == 0
            end
        end

        @testset "Operation $oper" for oper in (fast_add!, fast_sub!,
                                                fast_mul!)
            oper(C1, A, B)
            oper(CartesianBox, fill!(C2, 0), A, B)
            @test maxabsdif(C1, C2) ≤ 4*eps(T)
            if ndims(A) ≥ 1
                # Take a sub-region.
                S = CartesianBox(map(n -> (min(2,n) : max(n-1,1)), size(A)))
                oper(CartesianBox(S),     fill!(C2, 0), A, B)
                oper(CartesianIndices(S), fill!(C3, 0), A, B)
                @test maxabsdif(C2, C3) == 0
            end
         end
    end

    # FIXME: for now, we skip tests with ndims = 0
    @testset "Array functions" for dims in SIZES[2:end], T in (TYPES..., Bool)
        A = rand([zero(T), one(T)], dims)
        B = boundingbox(A)
        fill!(A, B, 0)
        C = boundingbox(A)
        @test isempty(C)
        @test length(C) == 0
        B = CartesianBox(map(r -> first(r)+1:last(r)-1, axes(A)))
        C = CartesianBox(map(r -> first(r)+2:last(r)-2, axes(A)))
        fill!(A, zero(T))
        fill!(A, B, one(T))
        @test boundingbox(A) === B
        A[B] .= zero(T)
        @test isempty(boundingbox(A))
        if T !== Bool
            X = rand(T, size(C))
            fill!(A, zero(T))
            fill!(A, B, typemax(T))
            A[C] = X
            @test all(map(isequal, A[C], X))
            @test boundingbox(x -> x != zero(T) && x < typemax(T), A) == C
            @test boundingbox(x -> x != typemax(T), A, B) == C
            inds = CartesianBoxes.indices(B)
            @test (@test_deprecated ranges(B)) === inds
            @test boundingbox(x -> x != typemax(T), A, inds) == C
            V = view(A,C)
            @test all(map(isequal, V, X))
            # Indices starting at 1.
            @. A = randn(T)
            inds = map(r -> first(r):last(r)-1, axes(A))
            B = CartesianBox(inds)
            @test maxabsdif(A[inds...], A[B]) == 0
        end
    end

    @testset "Miscellaneous" begin
        B = CartesianBox(1:3, -7:6)
        @test eachindex(IndexCartesian(), B) == B
        Z = CartesianBox(())
        @test stupidcount(Z) == 1
    end
end

end # module


text:
##################
## JFinM_Charts ##
##################

## handle data
##------------

import JFinM_Charts.writeData

@doc doc"""
Write data as .csv file. Hence, intended name needs not be modified.
"""->
function writeData(tm::AbstractTimenum,
                   chrt::MLineChart,
                   intendedFname::Array{ASCIIString, 1})

    if length(intendedFname) > 1
        error("Only a single data path is required.")
    end
    println("used data name: $(intendedFname[1])")
    writeTimedata(intendedFname[1], tm)
end