text:
export AuditLog

const AUDIT_LOG_CHANGE_TYPES = Dict(
    "name"                          => (String, Guild),
    "icon_hash"                     => (String, Guild) ,
    "splash_hash"                   => (String, Guild),
    "owner_id"                      => (Snowflake, Guild),
    "region"                        => (String, Guild),
    "afk_channel_id"                => (Snowflake, Guild),
    "afk_timeout"                   => (Int, Guild),
    "mfa_level"                     => (MFALevel, Guild),
    "verification_level"            => (VerificationLevel, Guild),
    "explicit_content_filter"       => (ExplicitContentFilterLevel, Guild),
    "default_message_notifications" => (MessageNotificationLevel, Guild),
    "vanity_url_code"               => (String, Guild),
    "\$add"                         => (Vector{Role}, Guild),
    "\$remove"                      => (Vector{Role}, Guild),
    "prune_delete_days"             => (Int, Guild),
    "widget_enabled"                => (Bool, Guild),
    "widget_channel_id"             => (Snowflake, Guild),
    "position"                      => (Int, DiscordChannel),
    "topic"                         => (String, DiscordChannel),
    "bitrate"                       => (Int, DiscordChannel),
    "permission_overwrites"         => (Vector{Overwrite}, DiscordChannel),
    "nsfw"                          => (Bool, DiscordChannel),
    "application_id"                => (Snowflake, DiscordChannel),
    "permissions"                   => (Int64, Role),
    "color"                         => (Int, Role),
    "hoist"                         => (Bool, Role),
    "mentionable"                   => (Bool, Role),
    "allow"                         => (Int, Role),
    "deny"                          => (Int, Role),
    "code"                          => (String, Invite),
    "channel_id"                    => (Snowflake, Invite),
    "inviter_id"                    => (Snowflake, Invite),
    "max_uses"                      => (Int, Invite),
    "uses"                          => (Int, Invite),
    "max_age"                       => (Int, Invite),
    "temporary"                     => (Bool, Invite),
    "deaf"                          => (Bool, User),
    "mute"                          => (Bool, User),
    "nick"                          => (String, User),
    "avatar_hash"                   => (String, User),
    "id"                            => (Snowflake, Any),
    "type"                          => (Any, Any),
    # Undocumented.
    "rate_limit_per_user"           => (Int, DiscordChannel),
)

"""
[`AuditLog`](@ref) action types.
More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-entry-object-audit-log-events).
"""
@enum ActionType  begin
    AT_GUILD_UPDATE=1
    AT_CHANNEL_CREATE=10
    AT_CHANNEL_UPDATE=11
    AT_CHANNEL_DELETE=12
    AT_CHANNEL_OVERWRITE_CREATE=13
    AT_CHANNEL_OVERWRITE_UPDATE=14
    AT_CHANNEL_OVERWRITE_DELETE=15
    AT_MEMBER_KICK=20
    AT_MEMBER_PRUNE=21
    AT_MEMBER_BAN_ADD=22
    AT_MEMBER_BAN_REMOVE=23
    AT_MEMBER_UPDATE=24
    AT_MEMBER_ROLE_UPDATE=25
    AT_ROLE_CREATE=30
    AT_ROLE_UPDATE=31
    AT_ROLE_DELETE=32
    AT_INVITE_CREATE=40
    AT_INVITE_UPDATE=41
    AT_INVITE_DELETE=42
    AT_WEBHOOK_CREATE=50
    AT_WEBHOOK_UPDATE=51
    AT_WEBHOOK_DELETE=52
    AT_EMOJI_CREATE=60
    AT_EMOJI_UPDATE=61
    AT_EMOJI_DELETE=62
    AT_MESSAGE_DELETE=72
end
@boilerplate ActionType :export :lower

"""
A change item in an [`AuditLogEntry`](@ref).

The first type parameter is the type of `new_value` and `old_value`. The second is the type
of the entity that `new_value` and `old_value` belong(ed) to.

More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-change-object).
"""
struct AuditLogChange{T, U}
    new_value::Optional{T}
    old_value::Optional{T}
    key::String
    type::Type{U}
end
@boilerplate AuditLogChange :docs :mock

AuditLogChange(d::Dict{Symbol, Any}) = AuditLogChange(; d...)
function AuditLogChange(; kwargs...)
    return if haskey(AUDIT_LOG_CHANGE_TYPES, kwargs[:key])
        T, U = AUDIT_LOG_CHANGE_TYPES[kwargs[:key]]
        func = if T === Any
            identity
        elseif T === Snowflake
            snowflake
        elseif T <: Vector
            eltype(T)
        else
            T
        end

        new_value = if haskey(kwargs, :new_value)
            if kwargs[:new_value] isa Vector
                func.(kwargs[:new_value])
            else
                func(kwargs[:new_value])
            end
        else
            missing
        end
        old_value = if haskey(kwargs, :old_value)
            if kwargs[:old_value] isa Vector
                func.(kwargs[:old_value])
            else
                func(kwargs[:old_value])
            end
        else
            missing
        end

        AuditLogChange{T, U}(new_value, old_value, kwargs[:key], U)
    else
        AuditLogChange{Any, Any}(
            get(kwargs, :new_value, missing),
            get(kwargs, :old_value, missing),
            kwargs[:key],
            Any,
        )
    end
end

function JSON.lower(x::AuditLogChange)
    d = Dict{Symbol, Any}()
    if !ismissing(x.new_value)
        d[:new_value] = x.new_value
    end
    if !ismissing(x.old_value)
        d[:old_value] = x.new_value
    end
    d[:key] = x.key
    return d
end

"""
Optional information in an [`AuditLogEntry`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-entry-object-optional-audit-entry-info).
"""
struct AuditLogOptions
    delete_member_days::Optional{Int}
    members_removed::Optional{Int}
    channel_id::Optional{Snowflake}
    count::Optional{Int}
    id::Optional{Snowflake}
    type::Optional{OverwriteType}
    role_name::Optional{String}
end
@boilerplate AuditLogOptions :docs :merge :mock

AuditLogOptions(d::Dict{Symbol, Any}) = AuditLogOptions(; d...)
function AuditLogOptions(; kwargs...)
    dmd = if haskey(kwargs, :delete_member_days)
        parse(Int, kwargs[:delete_member_days])
    else
        missing
    end
    return AuditLogOptions(
        dmd,
        haskey(kwargs, :members_removed) ? parse(Int, kwargs[:members_removed]) : missing,
        haskey(kwargs, :channel_id) ? snowflake(kwargs[:channel_id]) : missing,
        haskey(kwargs, :count) ? parse(Int, kwargs[:count]) : missing,
        haskey(kwargs, :id) ? snowflake(kwargs[:id]) : missing,
        haskey(kwargs, :type) ? OverwriteType(kwargs[:type]) : missing,
        get(kwargs, :role_name, missing),
    )
end

function JSON.lower(x::AuditLogOptions)
    d = Dict{Symbol, Any}()
    if !ismissing(x.delete_member_days)
        d[:delete_member_days] = string(d.delete_member_days)
    end
    ismissing(x.members_removed) || (d[:members_removed] = string(d.members_removed))
    ismissing(x.channel_id) || (d[:channel_id] = string(d.members_removed))
    ismissing(x.count) || (d[:count] = string(x.count))
    ismissing(x.id) || (d[:id] = string(x.id))
    ismissing(x.type) || (d[:type] = string(x.type))
    ismissing(x.role_name) || (d[:role_name] = x.role_name)
    return d
end

"""
An entry in an [`AuditLog`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-entry-object).
"""
struct AuditLogEntry
    target_id::Nullable{Snowflake}
    changes::Optional{Vector{AuditLogChange}}
    user_id::Snowflake
    id::Snowflake
    action_type::ActionType
    options::Optional{AuditLogOptions}
    reason::Optional{String}
end
@boilerplate AuditLogEntry :constructors :docs :lower :merge :mock

"""
An audit log.
More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-object).
"""
struct AuditLog
    webhooks::Vector{Webhook}
    users::Vector{User}
    audit_log_entries::Vector{AuditLogEntry}
end
@boilerplate AuditLog :constructors :docs :lower :merge :mock


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ROCmOpenCLRuntime"
version = v"4.2.0"

# Collection of sources required to build
sources = [
    ArchiveSource("https://github.com/ROCm-Developer-Tools/ROCclr/archive/rocm-$(version).tar.gz",
                  "c57525af32c59becf56fd83cdd61f5320a95024d9baa7fb729a01e7a9fcdfd78"),
    ArchiveSource("https://github.com/RadeonOpenCompute/ROCm-OpenCL-Runtime/archive/rocm-$(version).tar.gz",
                  "18133451948a83055ca5ebfb5ba1bd536ed0bcb611df98829f1251a98a38f730"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
OPENCL_SRC=$(realpath $WORKSPACE/srcdir/ROCm-OpenCL-Runtime-*)
ROCCLR_SRC=$(realpath $WORKSPACE/srcdir/ROCclr-*)

cd $ROCCLR_SRC
atomic_patch -p1 $WORKSPACE/srcdir/patches/rocclr-install-prefix.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/rocclr-link-lrt.patch
if [[ "${target}" == *-musl* ]]; then
atomic_patch -p1 $WORKSPACE/srcdir/patches/musl-rocclr.patch
atomic_patch -p1 $WORKSPACE/srcdir/patches/rocclr-disable-initial-exec.patch
fi
mkdir build && cd build
cmake -DCMAKE_PREFIX_PATH=${prefix} \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DOPENCL_DIR=$OPENCL_SRC \
      -DLLVM_DIR="${prefix}/lib/cmake/llvm" \
      -DClang_DIR="${prefix}/lib/cmake/clang" \
      -DLLD_DIR="${prefix}/lib/cmake/ldd" \
      ..
make -j${nproc}
make install

cd $OPENCL_SRC
if [[ "${target}" == *-musl* ]]; then
atomic_patch -p1 $WORKSPACE/srcdir/patches/musl-opencl.patch
fi
mkdir build && cd build
cmake -DCMAKE_PREFIX_PATH=${prefix} \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_BUILD_TYPE=Release \
      -DUSE_COMGR_LIBRARY=ON \
      -DBUILD_TESTS:BOOL=OFF \
      -DBUILD_TESTING:BOOL=OFF \
      -DLLVM_DIR="${prefix}/lib/cmake/llvm" \
      -DClang_DIR="${prefix}/lib/cmake/clang" \
      -DLLD_DIR="${prefix}/lib/cmake/ldd" \
      ..
make -j${nproc}
make install
install_license ../LICENSE.txt

# fix paths to point from src dir to install dir
sed -i "s#$ROCCLR_SRC/build#$WORKSPACE/destdir/lib#" $WORKSPACE/destdir/lib/cmake/rocclr/rocclr-targets.cmake
sed -i "s#$ROCCLR_SRC/include\;##" $WORKSPACE/destdir/lib/cmake/rocclr/rocclr-targets.cmake
sed -i "s#$ROCCLR_SRC#$WORKSPACE/destdir/include#g;s#$OPENCL_SRC.*\;##g" $WORKSPACE/destdir/lib/cmake/rocclr/rocclr-targets.cmake
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi="cxx11"),
    Platform("x86_64", "linux"; libc="musl", cxxstring_abi="cxx11"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libamdocl64"], :libamdocl),
    LibraryProduct(["libOpenCL"], :libOpenCL),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("hsa_rocr_jll", v"4.2.0"),
    Dependency("hsakmt_roct_jll", v"4.2.0"),
    Dependency("ROCmDeviceLibs_jll", v"4.2.0"),
    Dependency("ROCmCompilerSupport_jll", v"4.2.0"),
    BuildDependency(PackageSpec(; name="ROCmLLVM_jll", version="4.2.0")),
    Dependency("Libglvnd_jll"),
    Dependency("Xorg_libX11_jll"),
    Dependency("Xorg_xorgproto_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               julia_compat="1.7",
               preferred_gcc_version=v"8", preferred_llvm_version=v"11")

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the generation of random rotation representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/random.jl
# =====================

@testset "DCM (Float64)" begin
    D_exp = DCM(
         0.1770772447942165,
        -0.24937398754945297,
         0.9520799670772275,
        -0.16371156336807055,
        -0.9613536278168439,
        -0.22135430039413143,
         0.9704855348719379,
        -0.11666969022448989,
        -0.21105925705178208
    )

    D = rand(MersenneTwister(1986), DCM)
    @test D ≈ D_exp atol = 1e-15

    D = rand(MersenneTwister(1986), DCM{Float64})
    @test D ≈ D_exp atol = 1e-15
end

@testset "DCM (Float32)" begin
    D_exp = DCM(
         0.9071183f0,
        -0.38511035f0,
         0.1697833f0,
        -0.18077055f0,
         0.0077917147f0,
         0.98349446f0,
        -0.38007677f0,
        -0.9228377f0,
        -0.06254859f0
    )

    D = rand(MersenneTwister(1986), DCM{Float32})

    @test D ≈ D_exp atol = 1e-15
end

@testset "Euler angle and axis (Float64)" begin
    aa_exp = EulerAngleAxis(
        3.0732832539795845,
        @SVector [
             0.7668495688612252,
            -0.13482690249452656,
             0.6275057331220841
        ]
    )

    aa = rand(MersenneTwister(1986), EulerAngleAxis)
    @test aa ≈ aa_exp atol = 1e-15

    aa = rand(MersenneTwister(1986), EulerAngleAxis{Float64})
    @test aa ≈ aa_exp atol = 1e-15
end

@testset "Euler angle and axis (Float32)" begin
    aa_exp = EulerAngleAxis(
        1.6446829f0,
        @SVector [
            -0.9557738f0,
             0.2756822f0,
             0.102449425f0
        ]
    )

    aa = rand(MersenneTwister(1986), EulerAngleAxis{Float32})
    @test aa ≈ aa_exp atol = 1e-15
end

@testset "Euler angles (Float64)" begin
    ea_exp = EulerAngles(
        2.5852880822263606,
        3.1861195131389626,
        2.9299487570635576,
        :YXY
    )

    ea = rand(MersenneTwister(1986), EulerAngles)
    @test ea ≈ ea_exp atol = 1e-15

    ea = rand(MersenneTwister(1986), EulerAngles{Float64})
    @test ea ≈ ea_exp atol = 1e-15
end

@testset "Euler angles (Float32)" begin
    ea_exp = EulerAngles(
        0.29179642f0,
        5.512145f0,
        4.3565836f0,
        :YXY
    )

    ea = rand(MersenneTwister(1986), EulerAngles{Float32})
    @test ea ≈ ea_exp atol = 1e-15
end

@testset "Quaternion (Float64)" begin
    q_exp = Quaternion(
        -0.03414805970179898,
        -0.7664023306434491,
         0.1347482694144174,
        -0.6271397623279131
    )

    q = rand(MersenneTwister(1986), Quaternion)
    @test q ≈ q_exp atol = 1e-15

    q = rand(MersenneTwister(1986), Quaternion{Float64})
    @test q ≈ q_exp atol = 1e-15
end

@testset "Quaternion (Float32)" begin
    q_exp = Quaternion(
        -0.6805074f0,
         0.7003348f0,
        -0.2020037f0,
        -0.07506891f0
    )

    q = rand(MersenneTwister(1986), Quaternion{Float32})
    @test q ≈ q_exp atol = 1e-15
end


text:
using NestedEnvironments
using DifferentialEquations
using Plots
using Random
using Transducers

function dynamics(env::NestedEnvironments.InputAffineQuadraticCostEnv)
    return function (x, p, t)
        u = command(env, x)
        NestedEnvironments.ẋ(env, x, t, u)
    end
end
command(env, x) = NestedEnvironments.u_optimal(env, x)

initial_condition(env::NestedEnvironments.InputAffineQuadraticCostEnv) = 2*(rand(2) .- 0.5)

__env = NestedEnvironments.InputAffineQuadraticCostEnv()
__x0 = initial_condition(__env)
@reg_env __env __x0

## main
function single()
    Random.seed!(1)
    env = NestedEnvironments.InputAffineQuadraticCostEnv()
    x0 = initial_condition(env)
    t0 = 0.0
    tf = 10.0
    Δt = 0.01
    tspan = (t0, tf)
    ts = t0:Δt:tf
    prob = ODEProblem(env, dynamics(env), x0, tspan)
    @time sol = solve(prob, Tsit5(), saveat=ts)
    plot(sol)
end

function parallel()
    Random.seed!(1)
    env = NestedEnvironments.InputAffineQuadraticCostEnv()
    num = 10
    x0s = 1:num |> Map(i -> initial_condition(env))
    t0 = 0.0
    tf = 10.0
    Δt = 0.01
    tspan = (t0, tf)
    ts = t0:Δt:tf
    probs = x0s |> Map(x0 -> ODEProblem(env, dynamics(env), x0, tspan)) |> collect
    sols = probs |> Map(prob -> solve(prob, Tsit5(), saveat=ts)) |> tcollect
    p = plot()
    _ = sols |> Map(sol -> plot!(sol, label=nothing)) |> collect
    display(p)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
module Ising

struct Board{S, T, U}
    state::S
    β::T
    prob::U
end

function randboard(; n = 200, β = log(1 + √2)/2)
    state = rand(Int8.((-1, 1)), n, n)
    prob = Tuple(exp(-2k*β) for k in -4:4)
    Board(state, β, prob)
end

using Plots
@recipe function f(board::Board)
    seriestype := :heatmap
    board.state
end

P(i, m) = ifelse(i == m, 1, i + 1)
Q(i, m) = ifelse(i == 1, m, i - 1)

struct SingleThread end
function update!(board::Board, niters, ::SingleThread)
    s, prob = board.state, board.prob
    m, n = size(s)
    for _ in 1:niters
        for j in 1:n
            @inbounds for i in 1:m
                sij = s[i, j]
                k = sij * (s[Q(i, m), j] + s[P(i, m), j] + s[i, Q(j, n)] + s[i, P(j, n)])
                s[i, j] -= 2(rand() < prob[k + 5]) * sij
            end
        end
    end
end

struct MultiThread end
function update!(board::Board, niters, ::MultiThread)
    s, prob = board.state, board.prob
    m, n = size(s)
    for _ in 1:niters, j0 in 1:2
        Threads.@threads for j in j0:2:n
            @inbounds for i in 1:m
                sij = s[i, j]
                k = sij * (s[Q(i, m), j] + s[P(i, m), j] + s[i, Q(j, n)] + s[i, P(j, n)])
                s[i, j] -= 2(rand() < prob[k + 5]) * sij
            end
        end
    end
end

update!(board::Board, niters) = update!(board::Board, niters, MultiThread())

end

using Plots
default(size=(240, 240), colorbar=false, ticks=false, axis=false)

using BenchmarkTools

println("Julia v", VERSION)
println("nthreads() = ", Threads.nthreads())

# %%
board = Ising.randboard(n = 200)
@time anim = @animate for _ in 1:400
    Ising.update!(board, 1, Ising.SingleThread())
    plot(board)
end
gif(anim, "ising_singlethread.gif"; fps = 20)

# %%
board = Ising.randboard(n = 200)
@time anim = @animate for _ in 1:400
    Ising.update!(board, 1, Ising.MultiThread())
    heatmap(board)
end
gif(anim, "ising_multithread.gif"; fps = 20)

# %%
@show Threads.nthreads()
board = Ising.randboard(n = 200)
print("SingleThread:")
@btime Ising.update!($board, 1, $(Ising.SingleThread()))
print("MultiThread:  ")
@btime Ising.update!($board, 1, $(Ising.MultiThread()))

# %%
@show Threads.nthreads()
board = Ising.randboard(n = 632)
print("SingleThread:")
@btime Ising.update!($board, 1, $(Ising.SingleThread()))
print("MultiThread:   ")
@btime Ising.update!($board, 1, $(Ising.MultiThread()))

# %%

max_stars_repo_path: 0015/2D Ising multi-threaded.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/2D Ising multi-threaded.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/2D Ising multi-threaded.jl

max_forks_repo_name: genkuroki/public


text:
# classification
building(w::Way) = haskey(w.tags, "building")
highway(w::Way) = haskey(w.tags, "highway")

# predicates
visible{T <: OSMElement}(obj::T) = (get(obj.tags, "visible", "") != "false")
services(w::Way) = (get(w.tags,"highway", "") == "services")
reverse(w::Way) = (get(w.tags,"oneway", "") == "-1")

function haslanes(w::Way)
    v = get(w.tags, "lanes", "")
    length(v)==1 && '1' <= v[1] <= '9'
end

function oneway(w::Way)
    v = get(w.tags,"oneway", "")

    if v == "false" || v == "no" || v == "0"
        return false
    elseif v == "-1" || v == "true" || v == "yes" || v == "1"
        return true
    end

    highway = get(w.tags,"highway", "")
    junction = get(w.tags,"junction", "")

    return (highway == "motorway" ||
            highway == "motorway_link" ||
            junction == "roundabout")
end

# attributes
lanes(w::Way) = parse(Int, w.tags["lanes"])

# feature class
roadway(w::Way) = highway(w) * get(ROAD_CLASSES, w.tags["highway"], 0)

function walkway(way::Way)
    sidewalk = get(way.tags, "sidewalk", "")
    if sidewalk != "no"
        if haskey(PED_CLASSES, "sidewalk:$(sidewalk)")
            return PED_CLASSES["sidewalk:$(sidewalk)"]
        elseif haskey(PED_CLASSES, way.tags["highway"])
            return PED_CLASSES[way.tags["highway"]]
        end
    end
    return 0
end

function cycleway(w::Way)
    bicycle = get(w.tags, "bicycle", "")
    cycleway = get(w.tags, "cycleway", "")
    highway = get(w.tags, "highway", "")

    cycleclass = "cycleway:$(cycleway)"
    bikeclass = "bicycle:$(bicycle)"

    if bicycle != "no"
        if haskey(CYCLE_CLASSES, cycleclass)
            return CYCLE_CLASSES[cycleclass]
        elseif haskey(CYCLE_CLASSES, bikeclass)
            return CYCLE_CLASSES[bikeclass]
        elseif haskey(CYCLE_CLASSES, highway)
            return CYCLE_CLASSES[highway]
        end
    end
    return 0
end


text:
@testset "projectors" begin

#test project_bounds! with scalar bounds
  x = randn(100)
  l = -0.11
  u = 0.01;
  project_bounds!(x,l,u)
  @test maximum(x) <= u
  @test minimum(x) >= l

#test project_bounds! with vector bounds
  x = 100.0 .* randn(100)
  l = randn(100) .- 10.0
  u = randn(100) .+ 10.0
  project_bounds!(x,l,u)
  [@test x[i] <= u[i] for i=1:100]
  [@test x[i] >= l[i] for i=1:100]

#test project l1_Duchi!
  x=randn(100); tau=norm(x,1)*2; #see if projector returns untouched vector
  y=deepcopy(x)
  project_l1_Duchi!(x,tau)
  @test x == y

  x=randn(100); tau=norm(x,1)*0.234; #choose smaller constraint than vector l1 norm
  project_l1_Duchi!(x,tau)
  @test isapprox(norm(x,1),tau,rtol=10*eps())

  #test project l1_Duchi! with a complex vector
  x=randn(100)+im*randn(100); tau=norm(x,1)*2; #see if projector returns untouched vector
  y=deepcopy(x)
  project_l1_Duchi!(x,tau)
  @test x == y

  x=randn(100)+im*randn(100); tau=norm(x,1)*0.234;
  project_l1_Duchi!(x,tau)
  @test isapprox(norm(x,1),tau,rtol=10*eps())

#test project_cardinality!

  #test vector mode
  x=randn(100)
  project_cardinality!(x,5)
  @test count(!iszero, x)==5

  #test on closed form solution (vector)
  x=[0;0;1;2;3]
  project_cardinality!(x,2)
  @test x==[0;0;0;2;3]

  #test on closed form solution with negative numbers (vector)
  x=[0;0;-1;2;-3]
  project_cardinality!(x,2)
  @test x==[0;0;0;2;-3]

  #test mode that projects each column of a matrix separately
  X=randn(50,100)
  project_cardinality!(X,7,("fiber","x"))
  [@test count(!iszero, X[:,i]) == 7 for i=1:size(X,2)]

  #test mode that projects each row of a matrix separately
  X=randn(50,100)
  project_cardinality!(X,11,("fiber","z"))
  [@test count(!iszero, X[i,:]) == 11 for i=1:size(X,1)]

  #test mode that projects each fiber of a tensor separately (either in x,y or z direction)
  X=randn(50,60,30)
  project_cardinality!(X,7,("fiber","x"))
  [@test count(!iszero, X[:,i,j]) == 7 for i=1:size(X,2) for j=1:size(X,3)]

  X=randn(50,60,30)
  project_cardinality!(X,6,("fiber","y"))
  [@test count(!iszero, X[i,:,j]) == 6 for i=1:size(X,1) for j=1:size(X,3)]

  X=randn(50,60,30)
  project_cardinality!(X,4,("fiber","z"))
  [@test count(!iszero, X[i,j,:]) == 4 for i=1:size(X,1) for j=1:size(X,2)]

  #test mode that projects onto each slice of a tensor separately (either in x,y or z direction)
  X=randn(50,60,30)
  test=project_cardinality!(X,7,("slice","x"),false)
  [@test count(!iszero, test[i,:,:]) == 7 for i=1:size(X,1)]

  X=randn(50,60,30)
  test=project_cardinality!(X,6,("slice","y"),false)
  [@test count(!iszero, test[:,i,:]) == 6 for i=1:size(X,2)]

  X=randn(50,60,30)
  project_cardinality!(X,5,("slice","z"),false)
  [@test count(!iszero, X[:,:,i]) == 5 for i=1:size(X,3)]

# test project_l2!
  #test for sigma smaller than 2-norm of vector
  x=randn(100)
  project_l2!(x,0.123)
  @test isapprox(norm(x,2),0.123,rtol=eps()*10)

  #test for sigma larger than 2-norm of vector
  x=randn(100)
  y=deepcopy(x)
  project_l2!(x,1.234*norm(x,2))
  @test x==y

#test project_rank! (codes assumes input and output are vectors)
    X=randn(100,20)*randn(20,40) #randn(100,40) #test tall matrix
    grid_X=compgrid((1, 1),(100,40))
    Y=randn(40,50)*randn(50,100)#randn(40,100) #test flat matrix
    grid_Y=compgrid((1, 1),(40,100))
    Z=randn(100,30)*randn(30,100) #randn(100,100) #test square matrix
    grid_Z=compgrid((1, 1),(100,100))

    # test if input is returned untouched if input matrix rank is less than the constraint
    XX=vec(deepcopy(X))
    YY=vec(deepcopy(Y))
    ZZ=vec(deepcopy(Z))
    X=vec(X)
    Y=vec(Y)
    Z=vec(Z)
    project_rank!(reshape(X,grid_X.n),30,"")
    project_rank!(reshape(Y,grid_Y.n),40,"")
    project_rank!(reshape(Z,grid_Z.n),35,"")
    @test isapprox(XX,X,rtol=20*eps())
    @test isapprox(YY,Y,rtol=20*eps())
    @test isapprox(ZZ,Z,rtol=20*eps())

    # test if output matrix has correct rank when constraint is smaller than matrix rank
    X=randn(100,40) #test tall matrix
    Y=randn(40,100) #test flat matrix
    Z=randn(100,100) #test square matrix
    X=vec(X); Y=vec(Y); Z=vec(Z)
    project_rank!(reshape(X,grid_X.n),12,"")
    project_rank!(reshape(Y,grid_Y.n),3,"")
    project_rank!(reshape(Z,grid_Z.n),4,"")
    X=reshape(X,grid_X.n[1],grid_X.n[2])
    Y=reshape(Y,grid_Y.n[1],grid_Y.n[2])
    Z=reshape(Z,grid_Z.n[1],grid_Z.n[2])
    r_X=rank(X); r_Y=rank(Y); r_Z=rank(Z);
    @test r_X==12
    @test r_Y==3
    @test r_Z==4

    #test rank projection on slices of a 3D tensor
    X=randn(50,60,30)
    project_rank!(X,7,("slice","x"))
    [@test rank(X[i,:,:]) == 7 for i=1:size(X,1)]

    X=randn(50,60,30)
    project_rank!(X,6,("slice","y"))
    [@test rank(X[:,i,:]) == 6 for i=1:size(X,2)]

    X=randn(50,60,30)
    project_rank!(X,5,("slice","z"))
    [@test rank(X[:,:,i]) == 5 for i=1:size(X,3)]

#test project nuclear! (codes assumes input and output are vectors)
  X=randn(100,40) #test tall matrix
  grid_X=compgrid((1, 1),(100,40))
  Y=randn(40,100) #test flat matrix
  grid_Y=compgrid((1, 1),(40,100))
  Z=randn(100,100) #test square matrix
  grid_Z=compgrid((1, 1),(100,100))

  nn_X = norm(svdvals(X),1)
  nn_Y = norm(svdvals(Y),1)
  nn_Z = norm(svdvals(Z),1)

  # test if input is returned untouched if nuclear norm is less than constraint
  X=vec(X)
  Y=vec(Y)
  Z=vec(Z)
  XX=deepcopy(X)
  YY=deepcopy(Y)
  ZZ=deepcopy(Z)
  project_nuclear!(reshape(X,grid_X.n),nn_X*1.1,"")
  project_nuclear!(reshape(Y,grid_Y.n),nn_Y*1.1,"")
  project_nuclear!(reshape(Z,grid_Z.n),nn_Z*1.1,"")
  @test isapprox(XX,X,rtol=20*eps())
  @test isapprox(YY,Y,rtol=20*eps())
  @test isapprox(ZZ,Z,rtol=20*eps())

  # test if output matrix has correct nuclear norm when constraint is smaller than matrix nuclear norm
  X=randn(100,40) #test tall matrix
  Y=randn(40,100) #test flat matrix
  Z=randn(100,100) #test square matrix
  nn_X = norm(svdvals(X),1)
  nn_Y = norm(svdvals(Y),1)
  nn_Z = norm(svdvals(Z),1)
  X=vec(X)
  Y=vec(Y)
  Z=vec(Z)
  project_nuclear!(reshape(X,grid_X.n),nn_X*0.5,"")
  project_nuclear!(reshape(Y,grid_Y.n),nn_Y*0.5,"")
  project_nuclear!(reshape(Z,grid_Z.n),nn_Z*0.5,"")
  X=reshape(X,grid_X.n[1],grid_X.n[2])
  Y=reshape(Y,grid_Y.n[1],grid_Y.n[2])
  Z=reshape(Z,grid_Z.n[1],grid_Z.n[2])
  nn_Xp = norm(svdvals(X),1)
  nn_Yp = norm(svdvals(Y),1)
  nn_Zp = norm(svdvals(Z),1)
  @test isapprox(nn_Xp,nn_X*0.5,rtol=20*eps())
  @test isapprox(nn_Yp,nn_Y*0.5,rtol=20*eps())
  @test isapprox(nn_Zp,nn_Z*0.5,rtol=20*eps())

  #test nuclear norm projection onto slices of a 3D tensor
  X=randn(50,60,30)
  project_nuclear!(X,1.234,("slice","x"))
  [@test isapprox(norm(svdvals(X[i,:,:]),1), 1.234,rtol=100*eps()) for i=1:size(X,1)]

  X=randn(50,60,30)
  project_nuclear!(X,1.234,("slice","y"))
  [@test isapprox(norm(svdvals(X[:,i,:]),1), 1.234,rtol=100*eps()) for i=1:size(X,2)]

  X=randn(50,60,30)
  project_nuclear!(X,1.234,("slice","z"))
  [@test isapprox(norm(svdvals(X[:,:,i]),1), 1.234,rtol=100*eps()) for i=1:size(X,3)]

#test project_subspace!
  #on a vector inputs
  #orthogonal subspace
  M=randn(100,50)
  F=svd(M)
  x=randn(100)
  y=deepcopy(x)
  project_subspace!(x,F.U,true)
  @test isapprox(x,F.U*(F.U'*y),rtol=eps()*10)

  #non-orthogonal subspace
  M=randn(100,50)
  x=randn(100)
  y=deepcopy(x)
  project_subspace!(x,M,false)
  @test isapprox(x,M*((M'*M)\(M'*y)),rtol=eps()*10)

  #on matrix input: project every column onto the subspace
  M=randn(100,50)
  x=randn(100,5)
  y=deepcopy(x)
  project_subspace!(x,M,false,("fiber","x"))
  @test isapprox(x,M*((M'*M)\(M'*y)),rtol=eps()*10)

  #on matrix input: project every row onto the subspace
  M=randn(23,50)
  x=randn(5,23)
  y=deepcopy(x)
  project_subspace!(x,M,false,("fiber","z"))
  @test isapprox(x,(M*((M'*M)\(M'*y')))',rtol=eps()*10)

  #on tensor input: project every z slice onto the subspace
  x=randn(23,50,11)
  y=deepcopy(x)
  M=randn(23*50,5)
  project_subspace!(x,M,false,("slice","z"))
  x=reshape(x,size(y))
  [@test isapprox(vec(x[:,:,i]),M*((M'*M)\(M'*vec(y[:,:,i]))),rtol=eps()*100) for i=1:11]

  #on tensor input: project every y slice onto the subspace
  x=randn(23,11,50)
  y=deepcopy(x)
  M=randn(23*50,5)
  x=project_subspace!(x,M,false,("slice","y"))
  x=reshape(x,size(y))
  [@test isapprox(vec(x[:,i,:]),M*((M'*M)\(M'*vec(y[:,i,:]))),rtol=eps()*100) for i=1:11]

  #on tensor input: project every x slice onto the subspace
  x=randn(11,23,50)
  y=deepcopy(x)
  M=randn(23*50,5)
  x=project_subspace!(x,M,false,("slice","x"))
  x=reshape(x,size(y))
  [@test isapprox(vec(x[i,:,:]),M*((M'*M)\(M'*vec(y[i,:,:]))),rtol=eps()*100) for i=1:11]

#test projection onto relaxed histogram
  #first test exact histogram projection
  ref = sort(randn(100))
  x   = randn(100)
  project_histogram_relaxed!(x,ref,ref)
  @test ref == sort(x)

  #relaxed histogram projection
  LB = sort(randn(100))
  UB = LB.+0.7
  x   = randn(100)
  project_histogram_relaxed!(x,LB,UB)
  x = sort(x)
  [@test x[i] <= UB[i] for i=1:100]
  [@test x[i] >= LB[i] for i=1:100]

end


text:
function check_compatible_sizes(C, A, B)
    n, m = size(C)
    a, k = size(A)
    b, c = size(B)
    @assert (n == a) && (m == c) && (k == b) "matrices of size $(size(C)), $(size(A)), $(size(B)) are incompatible"
    nothing
end

function blocked_mul(A::MatTypes, B::MatTypes)
    T = promote_type(eltype(A), eltype(B))
    C = Matrix{T}(undef, size(A,1), size(B,2))
    blocked_mul!(C, A, B)
    C
end

choose_block_size(C::AbstractMatrix, ::Nothing) = size(C, 1) >= (3DEFAULT_BLOCK_SIZE) >>> 1 ? DEFAULT_BLOCK_SIZE : 32
choose_block_size(::Any, block_size::Integer) = block_size

function blocked_mul!(C::MatTypes{T}, A::MatTypes{T}, B::MatTypes{T};
              block_size = nothing, sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C, A, B)

    _block_size = choose_block_size(C, block_size)

    GC.@preserve C A B _mul!(PtrMatrix(C), PtrMatrix(A), PtrMatrix(B), _block_size)
    C
end

function blocked_mul(A::StructArray{Complex{T}, 2}, B::StructArray{Complex{T}, 2}) where {T <: Eltypes}
    C = StructArray{Complex{T}}((Matrix{T}(undef, size(A, 1), size(B,2)),
                                 Matrix{T}(undef, size(A, 1), size(B,2))))
    blocked_mul!(C, A, B)
    C
end

function blocked_mul!(C::StructArray{Complex{T}, 2}, A::StructArray{Complex{T}, 2}, B::StructArray{Complex{T}, 2};
              block_size = DEFAULT_BLOCK_SIZE, sizecheck=true) where {T <: Eltypes}
    sizecheck && check_compatible_sizes(C, A, B)

    _block_size = choose_block_size(C, block_size)

    GC.@preserve C A B begin
        Cre, Cim = PtrMatrix(C.re), PtrMatrix(C.im)
        Are, Aim = PtrMatrix(A.re), PtrMatrix(A.im)
        Bre, Bim = PtrMatrix(B.re), PtrMatrix(B.im)
        # Cre, Cim = C.re, C.im
        # Are, Aim = A.re, A.im
        # Bre, Bim = B.re, B.im
        _mul!(    Cre, Are, Bre, _block_size)          # C.re = A.re * B.re
        _mul_add!(Cre, Aim, Bim, _block_size, Val(-1)) # C.re = C.re - A.im * B.im
        _mul!(    Cim, Are, Bim, _block_size)          # C.im = A.re * B.im
        _mul_add!(Cim, Aim, Bre, _block_size)          # C.im = C.im + A.im * B.re
    end
    C
end

function _mul!(C, A, B, sz)
    n, k, m = size(C, 1), size(A, 2), size(C, 2)
    if n >= sz+8 && m >= sz+8 && k >= sz+8
        block_mat_mat_mul!(C, A, B, sz)
    elseif n >= sz+8 && k >= sz+8 && m <  sz+8
        block_mat_vec_mul!(C, A, B, sz)
    elseif n <  sz+8 && k >= sz+8 && m >= sz+8
        block_covec_mat_mul!(C, A, B, sz)
    elseif n >= sz+8 && k <  sz+8 && m >= sz+8
        block_vec_covec_mul!(C, A, B, sz)
    elseif n <  sz+8 && k >= sz+8 && m <  sz+8
        block_covec_vec_mul!(C, A, B, sz)
    else
        gemm_kernel!(C, A, B)
    end
end

function _mul_add!(C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    n, k, m = size(C, 1), size(A, 2), size(C, 2)
    if n >= sz+8 && m >= sz+8 && k >= sz+8
        block_mat_mat_mul_add!(C, A, B, sz, Val(factor))
    elseif n >= sz+8 && k >= sz+8 && m <  sz+8
        block_mat_vec_mul_add!(C, A, B, sz, Val(factor))
    elseif n <  sz+8 && k >= sz+8 && m >= sz+8
        block_covec_mat_mul_add!(C, A, B, sz, Val(factor))
    elseif n >= sz+8 && k <  sz+8 && m >= sz+8
        block_vec_covec_mul_add!(C, A, B, sz, Val(factor))
    elseif n <  sz+8 && k >= sz+8 && m <  sz+8
        block_covec_vec_mul_add!(C, A, B, sz, Val(factor))
    else
        add_gemm_kernel!(C, A, B, Val(factor))
    end
end


text:
__precompile__()
module ByteParsers

using Base.Dates: daysinmonth, UTD, totaldays
import Base: getindex, isequal, length, show, unsafe_get

using ArgCheck: @argcheck
using DocStringExtensions: SIGNATURES
using MacroTools: @capture
using Parameters: @unpack

export
    # generic
    parsedtype,
    parsenext,
    MaybeParsed,
    isparsed,
    getpos,
    # parsers
    PosInteger,
    PosFixedInteger,
    DateYYYYMMDD,
    Skip,
    ViewBytes,
    Line

# generic interface

const ByteVector = Vector{UInt8}

"""
A parser that parses to `MaybeParsed{T}` values.

Supports [`parsedtype`](@ref) and [`parsenext`](@ref) methods.
"""
abstract type AbstractParser{T} end

"""
    parsedtype(parser)

Return type of a parser, see [`parsenext`](@ref).

See [`Line`](@ref) for the semantics of skipping types which parse into `Void`.
"""
parsedtype(::AbstractParser{T}) where {T} = T

"""
    parsenext(parser, str, pos, sep)::MaybeParsed{parsedtype(parser)}

Parse `str` starting at `pos` using `parser`. Variable-length fields are parsed
until terminated by `sep`. Return a [`MaybeParsed`](@ref) object, the position
is that of the *next* character after parsing.
"""
function parsenext end

function parsenext(parser, str::AbstractString, pos, sep::Char)
    @argcheck Int(sep) < 0x80
    parsenext(parser, convert(ByteVector, str), pos, UInt8(sep))
end

"""
    MaybeParsed(pos, [value])

When `pos > 0`, it is the position of the character **after** parsing, and
`value` holds the value. When `pos < 0`, it wraps the error code (see
[`pos_to_error`](@ref)).
"""
struct MaybeParsed{T}
    pos::Int
    value::T
    MaybeParsed{T}(pos) where {T} = new{T}(pos)
    MaybeParsed{T}(pos, value::T) where {T} = new{T}(pos, value)
    MaybeParsed(pos, value::T) where {T} = new{T}(pos, value)
end

@inline unsafe_get(x::MaybeParsed) = x.value

"""
    $SIGNATURES

When [`isparsed(x)`](@ref), return the position of the *next* byte after parsing
the value that was returned as `x`.

Otherwise, return the position of the byte where parsing failed (may be outside
the length of the string for end of line errors).
"""
@inline getpos(x::MaybeParsed) = abs(x.pos)

function isequal(x::MaybeParsed{Tx}, y::MaybeParsed{Ty}) where {Tx,Ty}
    if isbits(Tx) && isbits(Ty)
        (x.pos == y.pos) & ifelse(x.pos > 0, isequal(x.value, y.value), true)
    else
        x.pos == y.pos && (x.pos > 0 ? isequal(x.value, y.value) : true)
    end
end

"""
    $SIGNATURES

Convert a position to an error code.
"""
@inline pos_to_error(pos::Int) = -pos

"""
    $SIGNATURES

Test if the object has been parsed (otherwise, it contains an error).
"""
@inline isparsed(x::MaybeParsed) = x.pos > 0

"""
    pos, value = @checkpos x [label]

Macro for checking a [`MaybeParsed`](@ref) object (eg returned by
[`parsenext`](@ref). If there is an error, `@goto label` is invoked, otherwise
the position and the parsed value are extracted into `pos` and `value`.
"""
macro checkpos(ex, label = :error)
    @capture ex (pos_, value_) = rhs_
    quote
        maybeparsed = $(esc(rhs))
        $(esc(pos)) = maybeparsed.pos
        if $(esc(pos)) < 0
            @goto $label
        else
            $(esc(value)) = maybeparsed.value
        end
    end
end

# integer parsing

struct PosInteger{T <: Integer, S <: Integer} <: AbstractParser{T} end

"""
    $SIGNATURES

Parse a positive integer until a delimiter as the first type. The second type is
used for the internal accumulation of digits, which is *not checked for
overflows*, but conversions to `T` are.
"""
function PosInteger(T::Type{<:Integer} = Int, S::Type{<:Integer} = Int64)
    if !(promote_type(T, S) ≡ S)
        warn("Using narrower type for parsing than return value.")
    end
    PosInteger{T, S}()
end

show(io::IO, parser::PosInteger{T, S}) where {T, S} =
     print(io, "Parse positive integer into $T (accumulating in $S)")

function parsenext(parser::PosInteger{T,S}, str::ByteVector, start::Int,
                   sep::UInt8) where {T,S}
    n = zero(S)
    z = UInt8('0')
    pos = start
    len = length(str)
    @inbounds while true
        chr = str[pos]
        if chr == sep
            if pos == start
                pos = pos_to_error(pos) # empty field
            else
                pos += 1
            end
            break
        end
        maybe_digit = chr - z
        if 0 ≤ maybe_digit ≤ 9
            n = n*10 + maybe_digit
            pos += 1
        else
            pos = pos_to_error(pos) # invalid character
            break
        end
        if pos > len            # EOL without separator
            pos = pos_to_error(pos)
            break
        end
    end
    MaybeParsed(pos, T(n))
end

struct PosFixedInteger{T <: Integer} <: AbstractParser{T}
    width::Int
    function PosFixedInteger{T}(width::Int) where {T}
        @argcheck width ≥ 1
        new{T}(width)
    end
end

show(io::IO, parser::PosFixedInteger{T}) where T =
     print(io, "Parse positive integer of $(parser.width) digits into $T")

"""
    $SIGNATURES

Parse a positive integer of *fixed width* as type `T`.

!!! NOTE
    Does not check for overflows, make sure `T` has enough digits.
"""
PosFixedInteger(width, T::Type{<: Integer} = Int) =
    PosFixedInteger{T}(width)

function parsenext(parser::PosFixedInteger{T}, str::ByteVector, start,
                   sep) where {T}
    n = zero(T)
    z = UInt8('0')
    pos = start
    stop = start + parser.width - 1
    len = length(str)
    stop > len && return MaybeParsed{Int}(pos_to_error(len+1))
    @inbounds while pos ≤ stop
        chr = str[pos]
        maybe_digit = chr - z
        if 0 ≤ maybe_digit ≤ 9
            n = n*10 + maybe_digit
            pos += 1
        else
            pos = pos_to_error(pos)
            break
        end
    end
    MaybeParsed(pos, n)
end

"""
    DateYYYYMMDD()

Parse a date in YYYYMMDD format.
"""
struct DateYYYYMMDD <: AbstractParser{Date} end

show(io::IO, ::DateYYYYMMDD) = print(io, "Parse date in YYYYMMDD format")

function parsenext(parser::DateYYYYMMDD, str::ByteVector, pos, sep)
    len = length(str)
    len ≥ pos + 8 || return MaybeParsed{Date}(pos_to_error(pos + 9)) # not enough characters
    @checkpos (pos, year) = parsenext(PosFixedInteger(4), str, pos, sep)
    @checkpos (pos, month) = parsenext(PosFixedInteger(2), str, pos, sep)
    1 ≤ month ≤ 12 || (pos -= 2; @goto error)
    @checkpos (pos, day) = parsenext(PosFixedInteger(2), str, pos, sep)
    1 ≤ day ≤ daysinmonth(year, month) || (pos -= 2; @goto error)
    @inbounds str[pos] == sep || @goto error
    pos += 1
    return MaybeParsed(pos, Date(UTD(totaldays(year, month, day))))
    @label error
    MaybeParsed{Date}(pos_to_error(pos))
end

# skip and view fields

"""
    Skip()

Skip a delimited field.
"""
struct Skip <: AbstractParser{Void} end

show(io::IO, ::Skip) = print(io, "skip field (parsed as `nothing`)")

function parsenext(parser::Skip, str::ByteVector, pos, sep::UInt8)
    len = length(str)
    @inbounds while pos ≤ len
        str[pos] == sep && return MaybeParsed(pos + 1, nothing)
        pos += 1
    end
    MaybeParsed{Void}(pos_to_error(pos))
end

"""
    ViewBytes()

Return the delimited field as a view (`SubArray`) into the parsed vector,
without the delimiter.
"""
struct ViewBytes <: AbstractParser{SubArray{UInt8,1,Array{UInt8,1},
                                            Tuple{UnitRange{Int64}},
                                            true}}
end

show(io::IO, ::ViewBytes) = print(io, "return field as bytes (SubArray)")

function parsenext(parser::ViewBytes, str::ByteVector, start, sep::UInt8)
    pos = start
    @checkpos (pos, value) = parsenext(Skip(), str, pos, sep)
    return MaybeParsed(pos, @view str[start:(pos-2)])
    @label error
    MaybeParsed{parsedtype(parser)}(pos)
end

# capture known values WIP

# struct Capture{T, S, P <: AbstractParser{S}} <: AbstractParser{Union{T,S}}
#     dict::Dict{ByteVector, S}
#     inner_parser::P
# end

# function _captured(dict::Dict{ByteVector, S}, str, start, sep) where S
#     for (key, value) in dict
#         len = length(key)
#         stop = start + len
#         if length(str) ≥ stop &&
#             str[stop] == sep &&
#             all(@view str[pos:(pos+len-1)] .== key)
#             return MaybeParsed(stop + 1, key)
#         end
#     end
#     MaybeParsed{S}(start)
# end

# function show(io::IO, parser::Capture)
#     @unpack dict, inner_parser = parser
#     print(io, "parser\n    capturing $(parser.dict)\n    inner parser ")
#     print(io, parser.inner)
# end

# function parsenext(parser::Capture, str::ByteVector, start, sep::UInt8)
#     @unpack values, inner_parser = parser
#     ifparsed
# end

# empty fields

struct FixEmpty{T, P <: AbstractParser{T}} <: AbstractParser{T}
    empty_value::T
    inner_parser::P
end

function show(io::IO, parser::FixEmpty)
    @unpack empty_value, inner_parser = parser
    print(io, "parse empty fields as $(empty_value), otherwise " * repr(inner_parser))
end

function parsenext(parser::FixEmpty, str::ByteVector, start, sep::UInt8)
    @unpack empty_value, inner_parser = parser
    if str[start] == sep
        MaybeParsed(start + 1, empty_value)
    else
        parsenext(inner_parser, str, start, sep)
    end
end

# lines

"""
    Line(parsers...)

Parse a string using `parsers`, separated by the separator given to
`parsenext`. Parsed values are returned as a tuple, wrapped in a `MaybeParsed`,
except for parsers which have `parsedtype(parser) == Void`, which are skipped.

# implementation

The return type is calculated and saved in parameter `S`, while the indices of
parsers for which the values are kept are in `K`. For example, for

```julia
Line(PosInteger(), Skip(), DateYYYYMMDD())
 ```

`S == Tuple{Int, Date}` and `K == (1, 3)`.
"""
struct Line{T <: Tuple, S, K} <: AbstractParser{S}
    parsers::T
end

function show(io::IO, parser::Line)
    print(io, "parse line as the following:")
    for p in parser.parsers
        print(io, "\n    ")
        print(io, p)
    end
end

function Line{T <: Tuple}(parsers::T)
    parsedtypes = map(parsedtype, parsers)
    keep(T) = T ≠ Void
    S = Tuple{Iterators.filter(keep, parsedtypes)...}
    K = tuple(find(keep, parsedtypes)...)
    Line{T, S, K}(parsers)
end

Line(parsers::AbstractParser...) = Line(parsers)

const LineN{N, S, K} = Line{T, S, K} where T <: NTuple{N, Any}

function _tuplefor(itr, ex)
    vars = [ Base.Cartesian.inlineanonymous(ex, i) for i in itr ]
    Expr(:escape, Expr(:tuple, vars...))
end

getindex(line::Line{T, S, K}, i) where {T,S,K} = line.parsers[K[i]]

length(line::Line{T, S, K}) where {T,S,K} = length(K)

"""
    @tuplefor itr expr

Generates a tuple using `Base.Cartesian.inlineanonymous`, using elements in
`itr` (cf `Base.Cartesian.@ntuple`, which is a special case of this for `itr =
1:N`). For example,

```julia @tuplefor (1,2,9) k ``` would generate ```julia (k_1, k_2, k_9) ```

!!! NOTE

    internal use, not exported.
"""
macro tuplefor(itr, ex)
    _tuplefor(itr, ex)
end

@generated function parsenext(line::LineN{N, S, K}, str::ByteVector,
                              pos, sep) where {N, S, K}
    quote
        R = MaybeParsed{S}
        Base.@nexprs $N j -> begin
            @checkpos (pos, value_j) = parsenext(line.parsers[j], str, pos, sep)
        end
        return R(pos, @tuplefor $K value)
        @label error
        R(pos)
    end
end

end # module


text:
module DFUtils
using DataFrames,
Parsers,
Missings,
Transducers

export complete,
    realtype,
    toReal,
    fixnothing!,
    readtypes,
    dictstodf

"""
    complete(df, cols...; replace_missing = missing)

Creates a `DataFrame` that has rows all combinations of `cols...`. Optionally fill missing values in other columns with `replace_missing`.
"""
function complete(df, cols...; replace_missing = missing)
    T = eltype(cols)
    dfo = Iterators.product([unique(df[:, col]) for col in cols]...) |> DataFrame
    rename!(dfo, Symbol.(1:length(cols)) .=> cols)
    dfo = rightjoin(df, dfo, on = [cols...])
    if !ismissing(replace_missing)
        for col in filter(n->T(n) ∉ cols, names(df))
            dfo[!, col] = coalesce.(dfo[!, col], replace_missing)
        end
    end
    return dfo
end

function realtype(s::S, r) where S <: AbstractString
    s = strip(s)
    o = Parsers.tryparse(Int, s)
    !isnothing(o) && return Int
    o = Parsers.tryparse(BigInt, s)
    !isnothing(o) && return BigInt
    p = Parsers.tryparse(Float64, s)
    if isnothing(p)
        return r ? Missing : S
    else
        return realtype(p)
    end
end

function realtype(n::T, r = true) where T <: Union{Integer, AbstractFloat}
    if round(n) == n
        return T <: Union{BigFloat, BigInt} ? BigInt : Int
    else
        return T
    end
end

realtype(m::Missing, r = true) = Missing

realtype(t::T, r = true) where T = T

"""
    realtype(v; replace_string = true)

Finds a subtype of `Real` that `v` can be converted to. Automatically handles `missing`s in vectors and optionally treats `String`s that cannot be parsed to any `Real` with `missing`
"""
function realtype(v::Vector; replace_string = true)
    T = Int
    for e in v
        T = promote_type(T, realtype(e, replace_string))
    end
    return T
end

toReal(::Type{T}, v::V, r=true) where {T<:Union{Integer, AbstractFloat}, V<:Union{Integer, AbstractFloat}} = convert(T,v)

toReal(::Type{T}, v::S, r) where {T<: AbstractString, S<: AbstractString} = r ? missing : v

function toReal(::Type{T}, v::S, r) where T<: Union{Integer, AbstractFloat} where S<:AbstractString
    if r
        o = Parsers.tryparse(T, v)
        isnothing(o) ? missing : o
    else
        return parse(T, v)
    end
end

toReal(::Type{Real}, v, r) = toReal(realtype(v), v, r)

function toReal(::Type{T}, v::S, r) where T<: BigFloat where S<:AbstractString
    if r
        try
            o = parse(T, v)
            return o
        catch e
            return missing
        end
    else
        return parse(T, v)
    end
end

toReal(::Type{T}, v::Missing, r = true) where T = missing
toReal(::Type{Missing}, v, r = true) = missing
toReal(::Type{Missing}, v::Missing, r = true) = missing

function toReal(::Type{T}, s::S, r) where T<:Integer where S<:AbstractString
    o = Parsers.tryparse(T, s)
    if !isnothing(o)
        return o
    else
        f = T <: BigInt ? toReal(BigFloat, s, r) : toReal(Float64, s, r)
        i = toReal(T, f)
        return i
    end
end

"""
    toReal(v; replace_string = true, threads = length(v) > 500)

Converts `v` to a `Vector{T}` where `T <: Real`. Automatically adds `Union{T, Missing}` if necessary and optionally replaces `Strings`s that cannot be parsed with `missing`.
"""
function toReal(v::Vector; replace_string = true, threads = length(v)>500)
    op = v |> Map(x -> toReal(realtype(x, replace_string), x, replace_string))
    if threads
        return tcollect(op)
    else
        return collect(op)
    end
end

toReal(x::Union{S, T, Missing}; replace_string = true) where S <: AbstractString where T <: Real = toReal(realtype(x), x, replace_string)

"""
    readtypes(U::Union, types = DataType[])

Get a vector of the types of a `Union` type
"""
function readtypes(U::Union, types = DataType[])
    push!(types, getfield(U, :a))
    if isa(getfield(U, :b), DataType)
        push!(types, getfield(U, :b))
        return types
    else
        readtypes(U.b, types)
    end
end

readtypes(T::DataType) = [T]

"""
    fixnothing!(df::DataFrame, col::Symbol)

Replace all `nothing` with `missing` in column and remove `Nothing` type.
"""
function fixnothing!(df::DataFrame, col)
    if Nothing ∈ readtypes(eltype(df[!, col]))
        allowmissing!(df, col)
        replace!(df[!, col], nothing => missing)
        df[!, col] = convert(Vector{Core.Compiler.typesubtract(eltype(df[!, col]), Nothing)}, df[!, col])
    else
        return nothing
    end
end

function fixnothing!(df::DataFrame)
    [fixnothing!(df, Symbol(col)) for col in names(df)]
    return nothing
end

"""
    dictstodf(dicts::Vector{Dict})

Create a DataFrame from a list of `Dict`s.
"""
function dictstodf(dicts)
    out = DataFrame()
    for l in dicts
        push!(out, l, cols = :union)
    end
    return out
end

end


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

JLLWrappers.@generate_main_file_header("Hwloc")
JLLWrappers.@generate_main_file("Hwloc", UUID("e33a78d0-f292-5ffc-b300-72abe9b543c8"))
end  # module Hwloc_jll


text:
@testset "Test How to Stop I" begin
###############################################################################
#
# The Stopping structure eases the implementation of algorithms and the
# stopping criterion.
# We illustrate here the basic features of Stopping.
#
# -> the case where a Stopping is a sub-Stopping is treated in the next tuto.
#
###############################################################################
#using Test, Stopping

x0 = ones(2)
pb = nothing

###############################################################################
#I. Initialize a Stopping
#The lazy way to initialize the stopping is to provide an initial point:
stop1 = GenericStopping(pb, x0, rtol = 1e-1)

#The more sophisticated way is to first build a State:
state1 = GenericState(ones(2))
#then, use it to create a Stopping:
stop2 = GenericStopping(pb, state1, rtol = 1e-1)

#Both ways give the same result:
@test stop1.current_state.x == stop2.current_state.x
@test isnan(stop1.current_state.current_time) && isnan(stop2.current_state.current_time)
#Keywords given in the Stopping creator are forwarded to the StoppingMeta.
@test stop1.meta.rtol == 1e-1

###############################################################################
#II. Check the status
#To ask the Stopping what is the current situation, we have the status function:
@test status(stop1) == :Unknown #nothing happened yet.
#The status function check the boolean values in the Meta:
#unbounded, unbounded_pb, tired, stalled, iteration_limit, resources, optimal,
#infeasible, main_pb, domainerror, suboptimal
stop1.meta.unbounded  = true
stop1.meta.suboptimal = true
#By default the status function prioritizes a status:
@test status(stop1) == :SubOptimal
#while you can access the list of status by turning the keyword list as true:
@test status(stop1, list =true) == [:SubOptimal, :Unbounded]

###############################################################################
#III. Analyze the situation: start!
#Two functions are designed to ask Stopping to analyze the current situation
#mainly described by the State: start!, stop!
#start! is designed to be used right at the beginning of the algorithm:
start!(stop1) #we will compare with stop2

#this call initializes a few entries:
#a) start_time in the META
@test isnan(stop2.meta.start_time)
@test !isnan(stop1.meta.start_time)
#b) optimality0 in the META (used to check the relative error)
@test stop2.meta.optimality0 == 1.0 #default value was 1.0
@test stop1.meta.optimality0 == 1.0 #GenericStopping has no specified measure, but get 1. if optimality is Inf
#c) the time measured is also updated in the State (if void)
@test stop1.current_state.current_time != nothing
#d) in the case where optimality0 is NaN, meta.domainerror becomes true
@test stop1.meta.domainerror == false
#e) the problem would be already solved if optimality0 pass a _null_test
#Since optimality0 is 1., any value would pass the relative error check:
@test Stopping._null_test(stop1, Inf) == false
@test stop1.meta.optimal == false
@test :SubOptimal in status(stop1, list = true)
#The Stopping determines the optimality by testing a score at zero.
#The test at zero is controlled by the function meta.tol_check which
#takes 3 arguments: atol, rtol, optimality0. By default it check if the score
#is less than: max(atol, rtol * opt0)
stop1.meta.tol_check = (atol, rtol, opt0) -> atol
@test Stopping._null_test(stop1, Inf) == false
#This can be determined in the initialization of the Stopping
stop3 = GenericStopping(pb, state1, tol_check = (atol, rtol, opt0) -> atol)
@test Stopping._null_test(stop3, Inf) == false

#The function _optimality_check providing the score returns Inf by default
#and must be specialized for specialized Stopping.
#If State entries have to be specified before the start!, you can use the
#function update_and_start! instead of a update! and then a start!
update_and_start!(stop3, x = zeros(2), current_time = -1.0)
@test stop3.meta.optimal == false
@test stop3.current_state.current_time == -1.0
@test stop3.meta.start_time != nothing
@test stop3.current_state.x == zeros(2)

###############################################################################
#Once the iterations begins #stop! is the main function.
#if needed an update is needed first, we can use update_and_stop!
OK = stop!(stop3) #update the Stopping and return  a boolean
@test OK == false #no reason to stop just yet!
#The stop! call check the following:
#1) meta.domainerror: check if the score is NaN
#2) meta.optimal: the score passes the _null_test
#3) meta.unbounded: check if state.x is too large
#4) meta.unbounded_pb: false by default
#5) meta.tired: check if time is exhausted
#6) meta.resources: false by default
#7) meta.iteration_limit: check the number of iterations
#8) meta.stalled: false by default
#9) meta.main_pb: false by default -> see Stopping as a subproblem tutorial
# Note that 1 and 2 are also done by start!.
#1) check unboundedness of x:
@test update_and_stop!(stop3, x = (stop3.meta.unbounded_x + 1.0) * x0 )
@test stop3.meta.unbounded == true
#5) check time
stop3.meta.start_time = 0.0 #too  force the time limit.
stop!(stop3)
@test stop3.meta.tired == true
#7) Stopping the number of iterations by the number of calls to stop!
@test stop3.meta.nb_of_stop == 3 #We called stop3 3 times already
stop3.meta.max_iter = 3
stop!(stop3)
@test stop3.meta.iteration_limit == true #as stop3.meta.nb_of_stop > 3.
#Overall we activated three flags:
@test status(stop3, list = true) == [:TimeLimit, :Unbounded, :IterationLimit]

###############################################################################
#Once we are done with an algorithm and want to reuse a stopping, we need to
#reinitialize all the entries.
reinit!(stop3)
#the status boolean are back to false
@test !stop3.meta.iteration_limit && !stop3.meta.tired && !stop3.meta.unbounded
#reinitialize also the entries updated by the start!
@test isnan(stop3.meta.start_time) && (stop3.meta.optimality0 == 1.0)
@test stop3.meta.nb_of_stop == 0 #and the counter of stop
#Note that by default reinit! does not reinitialize the current_state.
#This can be done by switching the keyword rstate to true.
#In this case, keywords are forwarded to the reinit! of current_state.
reinit!(stop3, rstate =  true, x = zeros(2))
@test isnan(stop3.current_state.current_time)
@test stop3.current_state.x == zeros(2)

end


text:
# SPDX-License-Identifier: X11
# 2020-10-18 (Timed)
# Tag (200pt)

A, V = parse.(Int, split(readline()))
B, W = parse.(Int, split(readline()))
T = parse(Int, readline())

S = V - W
D = abs(B - A)

if S <= 0
  println("NO")
elseif T * S < D
  println("NO")
else
  println("YES")
end

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
import LowRankModels: copy_estimate, copy

export copy_estimate, copy

function copy_estimate(g::GFRM)
  return GFRM(g.A,g.losses,g.r,g.k,
              g.observed_features,g.observed_examples,
              copy(g.U),copy(g.W))
end


text:
module StrFs

export StrF, @strf_str

using StaticArrays: SVector
using UnPack: @unpack

import Base: sizeof, read, write, isless, cmp, ==, typemin, repeat, promote_rule, show,
    codeunit, hash, length

"""
    StrF{S}(::String)

A string of less than `S` bytes, represented using 0-terminated UTF-8.

# Internal representation

`bytes` are stored as an `SVector{S, UInt8}`, where `S` is the maximum number of bytes. They
are not validated to be conforming UTF-8. A terminating `0x00` follows **if and only if**
the string is shorter than `S` bytes.

This follows the `str#` ("sturfs") representation of [Stata](https://www.stata.com/help.cgi?dta).

# Operations

A comparisons, conversion, and a few operations are supported, but this is primarily meant
as a storage type. For complex manipulations, it is recommended that you convert to
`String`.
"""
struct StrF{S} <: AbstractString
    bytes::SVector{S,UInt8}
end

macro strf_str(str)
    StrF(str)
end

show(io::IO, str::StrF) = show(io, String(str))

# this implementation is a modified copy from base/hashing2.jl
function hash(str::StrF, h::UInt)
    h += Base.memhash_seed
    # note: use pointer(s) here (see #6058).
    ccall(Base.memhash, UInt, (Ptr{UInt8}, Csize_t, UInt32),
          Base.cconvert(Ptr{UInt8}, str.bytes), sizeof(str), h % UInt32) + h
end

promote_rule(::Type{String}, ::Type{StrF{S}}) where S = String

promote_rule(::Type{StrF{A}}, ::Type{StrF{B}}) where {A,B} = StrF{max(A,B)}

codeunit(::StrF{S}) where S = UInt8

function sizeof(str::StrF{S}) where S
    nul = findfirst(isequal(0x00), str.bytes)
    if nul ≡ nothing
        S
    else
        nul - 1
    end
end

read(io::IO, ::Type{StrF{S}}) where S = StrF{S}(read(io, SVector{S,UInt8}))

write(io::IO, str::StrF) = write(io, str.bytes)

String(str::StrF{S}) where S = String(str.bytes[1:sizeof(str)])

function StrF{S}(str::AbstractString) where S
    @assert codeunit(str) ≡ UInt8
    len = sizeof(str)
    len ≤ S || throw(InexactError(:StrF, StrF{S}, str))
    bytes = Vector{UInt8}(undef, S)
    copyto!(bytes, codeunits(str))
    if len < S
        bytes[len + 1] = UInt8(0)
    end
    StrF{S}(SVector{S,UInt8}(bytes))
end

function StrF{A}(str::StrF{B}) where {A,B}
    if A == B
        str
    elseif A > B
        StrF{A}(vcat(str.bytes, zeros(SVector{A - B,UInt8})))
    elseif sizeof(str) ≤ A
        StrF{A}(str.bytes[1:A])
    else
        throw(InexactError(:StrF, StrF{A}, str))
    end
end

function StrF(str::AbstractString)
    @assert codeunit(str) ≡ UInt8
    S = sizeof(str)
    StrF{S}(SVector{S,UInt8}(codeunits(str)))
end

function cmp(a::StrF{A}, b::StrF{B}) where {A, B}
    for (achar, bchar) in zip(a.bytes, b.bytes)
        (achar == bchar == 0x0) && return 0
        c = cmp(achar, bchar)
        c ≠ 0 && return c
    end
    A == B && return 0
    (A < B && b.bytes[A+1] ≠ 0) && return -1
    (A > B && a.bytes[B+1] ≠ 0) && return 1
    0
end

cmp(a::StrF, b::AbstractString) = cmp(a, StrF(b)) # TODO improve

cmp(a::AbstractString, b::StrF) = -cmp(b, a)

isless(a::StrF, b::StrF) = cmp(a, b) < 0

(==)(a::StrF, b::StrF) = cmp(a, b) == 0

typemin(::StrF{S}) where S = StrF(zeros(SVector{S}))

length(str::StrF) = length(String(str)) # TODO improve

function repeat(str::StrF{S}, ::Val{n}) where {S, n}
    @unpack bytes = str
    s = sizeof(str)
    vS = n * s
    v = Vector{UInt8}(undef, vS)
    offset = 1
    for _ in 1:n
        copyto!(v, offset, bytes, 1, s)
        offset += s
    end
    if offset < vS
        v[offset] = 0
    end
    StrF{S}(SVector{S}(v))
end

function Base.iterate(str::StrF, state = (String(str), ))
    # NOTE: iteration implemented by converting to a string, and using it as the first
    # element of the state. The second element is the state for the iterator of the latter.
    y = iterate(state...)
    y ≡ nothing && return y
    first(y), (first(state), last(y))
end

Base.IteratorSize(::Type{<:StrF}) = Base.IteratorSize(String)

Base.IteratorEltype(::Type{<:StrF}) = Base.IteratorEltype(String)

Base.eltype(::Type{<:StrF}) = Base.eltype(String)

end # module


text:
"""
Thread-safe logging
"""
module Log

import Base.Threads.threadid

@inline puts(s) = ccall(:puts, Cint, (Ptr{Int8},), string(s))

function error(msg::String)
    puts("[$(threadid())] ERROR: $msg")
end

function warn(msg::String)
    puts("[$(threadid())] WARN: $msg")
end

function info(msg::String)
    puts("[$(threadid())] INFO: $msg")
end

function debug(msg::String)
    # In production mode, rather the development mode, don't log debug statements
    const is_production_run = haskey(ENV, "CELESTE_PROD") &&
                                     ENV["CELESTE_PROD"] != ""
    is_production_run || puts("[$(threadid())] DEBUG: $msg")
end

end


text:
module Sampling

export ListSampler, RejectionSampler, UniformSampler
export create_sampler

using Random
using DataStructures

using QXContexts.Contexts

# Module containing sampler objects which provide different levels of sampling features.
# Each sampler has a constructor which takes a context to perform sampling in and a set
# of keyword arguments that control the sampling behavior.
#
# Sampler(ctx; kwargs...): Initialise the sampler
#
# Each sampler is also callable with arguments that control it's execution
#
# (s::Sampler)(kwargs...): Perform sampling and return sampling results
#

"""Abstract type for samplers"""
abstract type AbstractSampler end

"""Functions to generate random bitstrings"""
random_bitstring(rng, num_qubits) = prod(rand(rng, ["0", "1"], num_qubits))
random_bitstrings(rng, num_qubits, num_samples) = [random_bitstring(rng, num_qubits) for _ in 1:num_samples]

###############################################################################
# ListSampler
###############################################################################

"""
A Sampler struct to compute the amplitudes for a list of bitstrings.
"""
struct ListSampler <: AbstractSampler
    ctx::AbstractContext
    list::Vector{String}
end

"""
    ListSampler(ctx
                ;bitstrings::Vector{String}=String[],
                rank::Integer=0,
                comm_size::Integer=1,
                kwargs...)

Constructor for a ListSampler to produce a portion of the given `bitstrings` determined by
the given `rank` and `comm_size`.
"""
function ListSampler(ctx
                     ;bitstrings::Vector{String}=String[],
                     kwargs...)
    if haskey(kwargs, :num_samples)
        n = kwargs[:num_samples]
        n = min(n, length(bitstrings))
    else
        n = length(bitstrings)
    end

    ListSampler(ctx, bitstrings[1:n])
end

"""
    (s::ListSampler)(max_amplitudes=nothing, kwargs...)

Callable for ListSampler struct. Calculates amplitudes for each bitstring in the list
"""
function (s::ListSampler)(;max_amplitudes=nothing, kwargs...)
    bs = if max_amplitudes === nothing
        s.list
    else s.list[1:min(max_amplitudes, length(s.list))] end

    amps = ctxmap(x -> compute_amplitude!(s.ctx, x; kwargs...), s.ctx, bs)
    amps = ctxgather(s.ctx, amps)
    if amps !== nothing return (bs, amps) end
end

create_sampler(ctx, sampler_params) = get_constructor(sampler_params[:method])(ctx ;sampler_params[:params]...)
get_constructor(func_name::String) = getfield(Main, Symbol(func_name*"Sampler"))

###############################################################################
# RejectionSampler
###############################################################################

"""
A Sampler struct to use rejection sampling to produce output.
"""
mutable struct RejectionSampler <: AbstractSampler
    ctx::AbstractContext
    num_qubits::Integer
    num_samples::Integer
    M::Real
    fix_M::Bool
    rng::MersenneTwister
end

"""
    function RejectionSampler(;num_qubits::Integer,
                              num_samples::Integer,
                              M::Real=0.0001,
                              fix_M::Bool=false,
                              seed::Integer=42,
                              kwargs...)

Constructor for a RejectionSampler to produce and accept a number of bitstrings.
"""
function RejectionSampler(ctx::AbstractContext;
                          num_qubits::Integer,
                          num_samples::Integer,
                          M::Real=0.0001,
                          fix_M::Bool=false,
                          seed::Integer=42,
                          kwargs...)
    # Evenly divide the number of bitstrings to be sampled amongst the subgroups of ranks.
    # num_samples = get_rank_size(num_samples, comm_size, rank)
    rng = MersenneTwister(seed) # TODO: should somehow add the rank to the seed, maybe with get_rank(ctx)?
    RejectionSampler(ctx, num_qubits, num_samples, M, fix_M, rng)
end

"""
    (s::RejectionSampler)(max_amplitudes=nothing, kwargs...)

Callable for RejectionSampler struct. Computes amplitudes for uniformly distributed bitstrings and corrects the distribution
using a rejection step.
"""
function (s::RejectionSampler)(;max_amplitudes=nothing, kwargs...)
    num_samples = max_amplitudes === nothing ? s.num_samples : max_amplitudes
    N = 2^s.num_qubits
    M = s.M

    samples = Samples()
    accepted = 0
    while accepted < num_samples
        # produce cadidate bitstrings
        bitstrings = random_bitstrings(s.rng, s.num_qubits, num_samples-accepted)

        # compute amplitudes for the bitstrings
        amps = [compute_amplitude!(s.ctx, bs; kwargs...) for bs in bitstrings]
        # bs_amp_pairs = [bs => compute_amplitude!(s.ctx, bs; kwargs...) for bs in bitstrings if !(bs in keys(samples.amplitudes))]

        # Record the computed amplitudes and update M if required
        for (bs, amp) in zip(bitstrings, amps)
            samples.amplitudes[bs] = amp
            Np = N * abs(amp)^2
            s.fix_M || (M = max(Np, M))
        end
        s.fix_M || (M = ctxreduce(max, s.ctx, M))

        # Conduct a rejection step for each bitstring to correct the distribution of samples.
        for (bs, amp) in zip(bitstrings, amps)
            Np = N * abs(amp)^2 # This is computed twice
            if rand(s.rng) < Np / M
                accepted += 1
                samples.bitstrings_counts[bs] += 1
            end
        end
    end

    ctxgather(s.ctx, samples)
    samples
end

###############################################################################
# UniformSampler
###############################################################################

"""
A Sampler struct to uniformly sample bitstrings and compute their amplitudes.
"""
mutable struct UniformSampler <: AbstractSampler
    ctx::AbstractContext
    num_qubits::Integer
    num_samples::Integer
    rng::MersenneTwister
end

"""
    UniformSampler(ctx::AbstractContext;
                    num_qubits::Integer,
                    num_samples::Integer,
                    seed::Integer=42,
                    kwargs...)

Constructor for a UniformSampler to uniformly sample bitstrings.
"""
function UniformSampler(ctx::AbstractContext;
                        num_qubits::Integer,
                        num_samples::Integer,
                        seed::Integer=42,
                        kwargs...)
    # Evenly divide the number of bitstrings to be sampled amongst the subgroups of ranks.
    # num_samples = (num_samples ÷ comm_size) + (rank < num_samples % comm_size)
    rng = MersenneTwister(seed)
    UniformSampler(ctx, num_qubits, num_samples, rng)
end

"""
    (s::UniformSampler)(max_amplitudes=nothing, kwargs...)

Callable for UniformSampler struct. Computes amplitudes for uniformly distributed bitstrings.
"""
function (s::UniformSampler)(;max_amplitudes=nothing, kwargs...)
    num_samples = max_amplitudes === nothing ? s.num_samples : max_amplitudes

    bs = random_bitstrings(s.rng, s.num_qubits, num_samples)
    amps = ctxmap(x -> compute_amplitude!(s.ctx, x; kwargs...), s.ctx, bs)
    amps = ctxgather(s.ctx, amps)

    (bs, amps)
end

###############################################################################
# Sampler Struct
###############################################################################

"""
Struct to hold the results of a simulation.
"""
struct Samples{T}
    bitstrings_counts::DefaultDict{String, <:Integer}
    amplitudes::Dict{String, T}
end

Samples() = Samples(DefaultDict{String, Int}(0), Dict{String, ComplexF32}())

Base.length(s::Samples) = sum(values(s.bitstrings_counts))
Base.unique(s::Samples) = keys(s.bitstrings_counts)

end


text:
using BenchmarkTools

function countelements_bilinear!(counts::AbstractVector{Int}, sequence)
    for idx in eachindex(counts)
        counts[idx] = count(==(idx), sequence)
    end
    return counts
end

function countelements_linear!(counts::AbstractVector{Int}, sequence)
    fill!(counts, 0)
    for item in sequence
        counts[item] += 1
    end
    return counts
end

function run_timing(m::Integer, n::Integer)
    sequence = rand(1:m, n)
    counts = Vector{Int}(undef, m)
    tq = @belapsed countelements_bilinear!($counts, $sequence) seconds=1
    tl = @belapsed countelements_linear!($counts, $sequence) seconds=1
    return tq, tl
end

function run_timing(ms, ns)
    tqs = Matrix{Float64}(undef, length(ms), length(ns))
    tls = similar(tqs)
    for (i, m) in enumerate(ms), (j, n) in enumerate(ns)
        tqs[i, j], tls[i, j] = run_timing(m, n)
    end
    return tqs, tls
end
function run_timing()
    ms, ns = 2 .^(1:10), 2 .^(10:20)
    return (ms, ns, run_timing(ms, ns)...)
end

function write_timing(filename, ms, ns, tqs, tls)
   open(filename, "w") do io
       for (i, m) in enumerate(ms)
            if i == 1
                for n in ns
                    print(io, ",", n)
                end
                println(io)
            end
            print(io, m)
            for j in eachindex(ns)
                print(io, ",(", tqs[i,j], ",", tls[i,j], ")")
            end
            println(io)
        end
    end
end

main(filename="timing_jl.csv") = write_timing(filename, run_timing()...)


text:
import StaticArrays: MArray, MVector, SVector

###################################################################################
# Mutable simplex where vertices are represented by StaticArrays.MVector instances.
###################################################################################

"""
    MutableSSimplex(pts::Vector{SVector})

Construct a mutable `SSimplex` from a vector of `SVector`s.
"""
struct MutableSSimplex{dim, T} <: AbstractMutableSimplex{dim, T}
    vertices::Vector{MVector{dim, T}}

    function MutableSSimplex(pts::Vector{MVector{dim, T}}) where {dim, T}
        if !(length(pts) == dim + 1)
            err = """ The input cannot be converted to a simplex.
            Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
            """
            throw(DomainError(pts, err))
        end

        new{dim, T}(pts)
    end
end

function Base.setindex!(s::MutableSSimplex, v::SVector{D, T}, i::Int) where {D, T}
    s[i] = v
end

function Base.setindex!(s::MutableSSimplex, v::SVector{D, T}, i::Int, j) where {D, T}
    s[i][j] = v
end

export MutableSSimplex


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1PodDNSConfigOption <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    value::Any # spec type: Union{ Nothing, String } # spec name: value

    function IoK8sApiCoreV1PodDNSConfigOption(;name=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodDNSConfigOption, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1PodDNSConfigOption, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiCoreV1PodDNSConfigOption

const _property_map_IoK8sApiCoreV1PodDNSConfigOption = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiCoreV1PodDNSConfigOption = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("value")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PodDNSConfigOption }) = collect(keys(_property_map_IoK8sApiCoreV1PodDNSConfigOption))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodDNSConfigOption }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1PodDNSConfigOption[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodDNSConfigOption }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodDNSConfigOption[property_name]

function check_required(o::IoK8sApiCoreV1PodDNSConfigOption)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodDNSConfigOption }, name::Symbol, val)
end


text:
# Author: Lukasz Bratos

# Typy zmiennopozycyjne dla ktorych bedziemy wykonywac testy
types = [Float16, Float32, Float64]

# Funkcja wyznacza iteracyjnie epsilon maszynowy dla danego typu zmiennopozycyjnego.
function macheps(type)
    x = one(type)
    while one(type) + x / 2 != one(type)
        x = x / 2
    end
    return x
end

# Funkcja wyznacza iteracyjnie wartosc eta dla danego typy zmiennopozycyjnego.
function eta(type)
    x = one(type)
    while x / 2 > 0
        x = x / 2
    end
    return x
end

# Funkcja wyznacza iteracyjnie maksymalna wartosc dla danego typu zmiennopozycyjnego.
function maxFloat(type)
    x = one(type)
    while !isinf(x * 2)
        x *= 2
    end
    y = x / 2
    while !isinf(x + y) && y >= 1.0
        x += y
        y /= 2
    end
    return x
end

for type in types
    println("Calculated macheps for $type: ", macheps(type), " vs eps($type): ", eps(type))
    println("Calculated eta for $type: ", eta(type), " vs nextfloat($type): ", nextfloat(type(0.0)))
    println("Calculated MAX for $type: ", maxFloat(type), " vs floatmax($type): ", floatmax(type))
end


text:
using Pkg
using Test
using MeshIO
using StaticArrays
using Makie
using ImageMagick

using Makie.Observables
using Makie.GeometryBasics
using Makie.PlotUtils
using Makie.FileIO
using Makie.IntervalSets
using GeometryBasics: Pyramid

# ImageIO seems broken on 1.6 ... and there doesn't
# seem to be a clean way anymore to force not to use a loader library?
filter!(x-> x !== :ImageIO, FileIO.sym2saver[:PNG])
filter!(x-> x !== :ImageIO, FileIO.sym2loader[:PNG])

include("reference_tests.jl")
include("unit_tests/runtests.jl")


text:
# discovering binary CUDA dependencies

using Pkg, Pkg.Artifacts
import Libdl

## global state

const __toolkit_version = Ref{VersionNumber}()

"""
    version()

Returns the version of the CUDA toolkit in use.
"""
toolkit_version() = @after_init(__toolkit_version[])

"""
    release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
toolkit_release() = @after_init(VersionNumber(__toolkit_version[].major, __toolkit_version[].minor))

const __nvdisasm = Ref{String}()
const __libdevice = Ref{String}()
const __libcudadevrt = Ref{String}()
const __libcupti = Ref{Union{Nothing,String}}()
const __libnvtx = Ref{Union{Nothing,String}}()
const __libcublas = Ref{String}()
const __libcusparse = Ref{String}()
const __libcusolver = Ref{String}()
const __libcufft = Ref{String}()
const __libcurand = Ref{String}()
const __libcudnn = Ref{Union{Nothing,String}}(nothing)
const __libcutensor = Ref{Union{Nothing,String}}(nothing)

nvdisasm() = @after_init(__nvdisasm[])
libdevice() = @after_init(__libdevice[])
libcudadevrt() = @after_init(__libcudadevrt[])
function libcupti()
    @after_init begin
        @assert has_cupti() "This functionality is unavailable as CUPTI is missing."
        __libcupti[]
    end
end
function libnvtx()
    @after_init begin
        @assert has_nvtx() "This functionality is unavailable as NVTX is missing."
        __libnvtx[]
    end
end

export has_cupti, has_nvtx
has_cupti() = @after_init(__libcupti[]) !== nothing
has_nvtx() = @after_init(__libnvtx[]) !== nothing

libcublas() = @after_init(__libcublas[])
libcusparse() = @after_init(__libcusparse[])
libcusolver() = @after_init(__libcusolver[])
libcufft() = @after_init(__libcufft[])
libcurand() = @after_init(__libcurand[])
function libcudnn()
    @after_init begin
        @assert has_cudnn() "This functionality is unavailabe as CUDNN is missing."
        __libcudnn[]
    end
end
function libcutensor()
    @after_init begin
        @assert has_cutensor() "This functionality is unavailabe as CUTENSOR is missing."
        __libcutensor[]
    end
end

export has_cudnn, has_cutensor
has_cudnn() = @after_init(__libcudnn[]) !== nothing
has_cutensor() = @after_init(__libcutensor[]) !== nothing

## discovery

# CUDA

# NOTE: we don't use autogenerated JLLs, because we have multiple artifacts and need to
#       decide at run time (i.e. not via package dependencies) which one to use.
const cuda_artifacts = Dict(
    v"10.2" => ()->artifact"CUDA10.2",
    v"10.1" => ()->artifact"CUDA10.1",
    v"10.0" => ()->artifact"CUDA10.0",
    v"9.2"  => ()->artifact"CUDA9.2",
    v"9.0"  => ()->artifact"CUDA9.0",
)

# utilities to look up stuff in the artifact (at known locations)
get_binary(artifact_dir, name) = joinpath(artifact_dir, "bin", Sys.iswindows() ? "$name.exe" : name)
function get_library(artifact_dir, name)
    filename = if Sys.iswindows()
        "$name.dll"
    elseif Sys.isapple()
        "lib$name.dylib"
    else
        "lib$name.so"
    end
    joinpath(artifact_dir, Sys.iswindows() ? "bin" : "lib", filename)
end
get_static_library(artifact_dir, name) = joinpath(artifact_dir, "lib", Sys.iswindows() ? "$name.lib" : "lib$name.a")
get_file(artifact_dir, path) = joinpath(artifact_dir, path)

function use_artifact_cuda()
    @debug "Trying to use artifacts..."

    # select compatible artifacts
    if haskey(ENV, "JULIA_CUDA_VERSION")
        wanted_version = VersionNumber(ENV["JULIA_CUDA_VERSION"])
        filter!(((version,artifact),) -> version == wanted_version, cuda_artifacts)
    else
        driver_version = release()
        filter!(((version,artifact),) -> version <= driver_version, cuda_artifacts)
    end

    # download and install
    artifact = nothing
    for release in sort(collect(keys(cuda_artifacts)); rev=true)
        try
            artifact = (release=release, dir=cuda_artifacts[release]())
            break
        catch
        end
    end
    if artifact == nothing
        @debug "Could not find a compatible artifact."
        return false
    end

    __nvdisasm[] = get_binary(artifact.dir, "nvdisasm")
    @assert isfile(__nvdisasm[])
    __toolkit_version[] = parse_toolkit_version(__nvdisasm[])

    # Windows libraries are tagged with the CUDA release
    long = "$(artifact.release.major)$(artifact.release.minor)"
    short = artifact.release >= v"10.1" ? string(artifact.release.major) : long

    __libcupti[] = get_library(artifact.dir, Sys.iswindows() ? "cupti64_$long" : "cupti")
    @assert isfile(__libcupti[])
    __libnvtx[] = get_library(artifact.dir, Sys.iswindows() ? "nvToolsExt64_1" : "nvToolsExt")
    @assert isfile(__libnvtx[])

    __libcudadevrt[] = get_static_library(artifact.dir, "cudadevrt")
    @assert isfile(__libcudadevrt[])
    __libdevice[] = get_file(artifact.dir, joinpath("share", "libdevice", "libdevice.10.bc"))
    @assert isfile(__libdevice[])

    for name in  ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CUDA, Symbol("__lib$name"))

        # on Windows, the library name is version dependent
        if Sys.iswindows()
            name *= if artifact.release >= v"10.1"
                "64_$(artifact.release.major)"
            else
                "64_$(artifact.release.major)$(artifact.release.minor)"
            end
        end

        handle[] = get_library(artifact.dir, name)
        Libdl.dlopen(handle[])
    end

    @debug "Using CUDA $(__toolkit_version[]) from an artifact at $(artifact.dir)"
    use_artifact_cudnn(artifact.release)
    use_artifact_cutensor(artifact.release)
    return true
end

function use_local_cuda()
    @debug "Trying to use local installation..."

    cuda_dirs = find_toolkit()

    let path = find_cuda_binary("nvdisasm", cuda_dirs)
        if path === nothing
            @debug "Could not find nvdisasm"
            return false
        end
        __nvdisasm[] = path
    end

    cuda_version = parse_toolkit_version(__nvdisasm[])
    __toolkit_version[] = cuda_version

    cupti_dirs = map(dir->joinpath(dir, "extras", "CUPTI"), cuda_dirs) |> x->filter(isdir,x)
    __libcupti[] = find_cuda_library("cupti", [cuda_dirs; cupti_dirs], [cuda_version])
    __libnvtx[] = find_cuda_library("nvtx", cuda_dirs, [v"1"])

    let path = find_libcudadevrt(cuda_dirs)
        if path === nothing
            @debug "Could not find libcudadevrt"
            return false
        end
        __libcudadevrt[] = path
    end
    let path = find_libdevice(cuda_dirs)
        if path === nothing
            @debug "Could not find libdevice"
            return false
        end
        __libdevice[] = path
    end

    for name in  ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CUDA, Symbol("__lib$name"))

        path = find_cuda_library(name, cuda_dirs, [cuda_version])
        if path === nothing
            @debug "Could not find $name"
            return false
        end
        handle[] = path
    end

    @debug "Found local CUDA $(cuda_version) at $(join(cuda_dirs, ", "))"
    use_local_cudnn(cuda_dirs)
    use_local_cutensor(cuda_dirs)
    return true
end

# CUDNN

const cudnn_artifacts = Dict(
    v"10.2" => ()->artifact"CUDNN+CUDA10.2",
    v"10.1" => ()->artifact"CUDNN+CUDA10.1",
    v"10.0" => ()->artifact"CUDNN+CUDA10.0",
    v"9.2"  => ()->artifact"CUDNN+CUDA9.2",
    v"9.0"  => ()->artifact"CUDNN+CUDA9.0",
)

function use_artifact_cudnn(release)
    artifact_dir = try
        cudnn_artifacts[release]()
    catch ex
        @debug "Could not use CUDNN from artifacts" exception=(ex, catch_backtrace())
        return false
    end

    __libcudnn[] = get_library(artifact_dir, Sys.iswindows() ? "cudnn64_7" : "cudnn")
    Libdl.dlopen(__libcudnn[])
    @debug "Using CUDNN from an artifact at $(artifact_dir)"
    return true
end

function use_local_cudnn(cuda_dirs)
    path = find_cuda_library("cudnn", cuda_dirs, [v"7"])
    path === nothing && return false

    __libcudnn[] = path
    @debug "Using local CUDNN at $(path)"
    return true
end

# CUTENSOR

const cutensor_artifacts = Dict(
    v"10.2" => ()->artifact"CUTENSOR+CUDA10.2",
    v"10.1" => ()->artifact"CUTENSOR+CUDA10.1",
)

function use_artifact_cutensor(release)
    artifact_dir = try
        cutensor_artifacts[release]()
    catch ex
        @debug "Could not use CUTENSOR from artifacts" exception=(ex, catch_backtrace())
        return false
    end

    __libcutensor[] = get_library(artifact_dir, "cutensor")
    Libdl.dlopen(__libcutensor[])
    @debug "Using CUTENSOR from an artifact at $(artifact_dir)"
    return true
end

function use_local_cutensor(cuda_dirs)
    path = find_cuda_library("cutensor", cuda_dirs, [v"1"])
    path === nothing && return false

    __libcutensor[] = path
    @debug "Using local CUTENSOR at $(path)"
    return true
end

function __init_dependencies__()
    found = false

    if parse(Bool, get(ENV, "JULIA_CUDA_USE_BINARYBUILDER", "true"))
        found = use_artifact_cuda()
    end

    if !found
        found = use_local_cuda()
    end

    return found
end


text:
# indicator of the ball of matrices with (at most) a given rank

"""
  IndBallRank(r::Int=1)

Returns the function `g = ind{X : rank(X) ⩽ r}`, for an integer parameter `r > 0`.
"""

immutable IndBallRank{I <: Integer} <: IndicatorNonconvex
  r::I
  function IndBallRank(r::I)
    if r <= 0
      error("parameter r must be a positive integer")
    else
      new(r)
    end
  end
end

IndBallRank{I <: Integer}(r::I=1) = IndBallRank{I}(r)

function (f::IndBallRank){T <: RealOrComplex}(x::AbstractArray{T,2})
  maxr = minimum(size(x))
  if maxr <= f.r return 0.0 end
  svdobj = svds(x, nsv=f.r+1)[1]
  # the tolerance in the following line should be customizable
  if svdobj[:S][end]/svdobj[:S][1] <= 1e-14
    return 0.0
  end
  return +Inf
end

function prox!{T <: Real}(f::IndBallRank, x::AbstractArray{T,2}, y::AbstractArray{T,2}, gamma::Real=1.0)
  maxr = minimum(size(x))
  if maxr <= f.r
    y[:] = x
    return 0.0
  end
  svdobj = svds(x, nsv=f.r)[1]
  for i = 1:size(x,1)
    for j = 1:size(x,2)
      y[i,j] = 0.0
      for k = 1:f.r
          if VERSION >= v"0.6.0-dev.2071"
              y[i,j] += svdobj[:U][i,k]*svdobj[:S][k]*svdobj[:Vt][k,j]
          else
              y[i,j] += svdobj[:U][i,k]*svdobj[:S][k]*svdobj[:Vt][j,k]
          end
      end
    end
  end
  return 0.0
end

function prox!{T <: Complex}(f::IndBallRank, x::AbstractArray{T,2}, y::AbstractArray{T,2}, gamma::Real=1.0)
  maxr = minimum(size(x))
  if maxr <= f.r
    y[:] = x
    return 0.0
  end
  svdobj = svds(x, nsv=f.r)[1]
  for i = 1:size(x,1)
    for j = 1:size(x,2)
      y[i,j] = 0.0 + 0.0im
      for k = 1:f.r
          if VERSION >= v"0.6.0-dev.2071"
              y[i,j] += svdobj[:U][i,k]*svdobj[:S][k]*svdobj[:Vt][k,j]
          else
              y[i,j] += svdobj[:U][i,k]*svdobj[:S][k]*conj(svdobj[:Vt][j,k])
          end
      end
    end
  end
  return 0.0
end

fun_name(f::IndBallRank) = "indicator of the set of rank-r matrices"
fun_dom(f::IndBallRank) = "AbstractArray{Real,2}, AbstractArray{Complex,2}"
fun_expr(f::IndBallRank) = "x ↦ 0 if rank(x) ⩽ r, +∞ otherwise"
fun_params(f::IndBallRank) = "r = $(f.r)"

function prox_naive{T <: RealOrComplex}(f::IndBallRank, x::AbstractArray{T,2}, gamma::Real=1.0)
  maxr = minimum(size(x))
  if maxr <= f.r
    y = x
    return y, 0.0
  end
  U, S, V = svd(x)
  M = U[:,1:f.r]*spdiagm(S[1:f.r])
  y = M*V[:,1:f.r]'
  return y, 0.0
end


text:
using TeraSort
yarncm = startcluster(true)
try
    @everywhere using TeraSort
    teragen()
    terasample()
finally
    stopcluster(yarncm)
end


text:
"""
    walkbacktrellis(x,message,h)

    Return the the stego signal with embedded message
"""
function walkbacktrellis(path,x,message,hhat,h)
    numofblocks=length(message);
    indx=length(x)
    indm=length(message)
    state=0
    nstates=2^h-1;

    y=zeros(UInt8,length(x))  # the stego signal
    @inbounds for i = numofblocks:-1:1
        hh=(length(message)-indm+1<h) ? hh=2^(length(message)-indm+1)-1 : nstates;
        state = (2*state + message[indm])&nstates
        indm-=1
        @inbounds for j = length(hhat):-1:1
            y[indx] = path[state+1,indx]
            state = state ⊻ (y[indx]*(hhat[j]&hh))
            indx-=1
        end
    end
    return(y)
end

"""
    extract(y,hhat,h)

    Extract the message from stego object y.

"""
function extract(y,hhat,h)
    numofblocks=div(length(y),length(hhat))
    state=0
    indy=1
    exm=zeros(Int,numofblocks)
    for b in 1:numofblocks
        for i=1:length(hhat)
            state=state ⊻ (y[indy]*hhat[i])
            indy+=1
        end
        exm[b]=mod(state,2)
        state>>=1
    end
    return(exm)
end

function swap(a,b)
    return(b,a)
end


text:
@testset "nonparametric" begin
    @test_throws ArgumentError Periodogram(ones(1000), -1)
    @test_throws ArgumentError Periodogram(ones(1000,2), -1)
    @test Periodogram(rand(1000),1/2) isa Periodogram
    @test BartlettPeriodogram(rand(1000),1/2) isa BartlettPeriodogram
    @test Periodogram(rand(1000,2),1/2) isa Periodogram
    @test BartlettPeriodogram(rand(1000,2),1/2) isa BartlettPeriodogram
    @test CoherancyEstimate(rand(1000,2),1/2) isa CoherancyEstimate
    @test Periodogram(rand(1000),1/2) + Periodogram(rand(1000),1/2) isa Periodogram
    @test BartlettPeriodogram(rand(1000),1/2) + BartlettPeriodogram(rand(1000),1/2) isa BartlettPeriodogram
    @test Periodogram(rand(1000,2),1/2) + Periodogram(rand(1000,2),1/2) isa Periodogram
    @test BartlettPeriodogram(rand(1000,2),1/2) + BartlettPeriodogram(rand(1000,2),1/2) isa BartlettPeriodogram
    @test CoherancyEstimate(rand(1000,2),1/2) + CoherancyEstimate(rand(1000,2),1/2) isa CoherancyEstimate
end


text:
# Completes Nt loops of loopNs!, beginning at current optimum
@inline function loopNt!(current::State, result::Result, options::Options)
    current.f = result.fopt
    current.x .= result.xopt
    @inbounds for _ in 1:options.Nt
        # loopNs and adjust search space
        loopNs!(current, result, options)
        adjust_step!(current, options)

        # Reset accepted points for search space
        current.n_accepted .= 0
    end

    return nothing
end

# Parallel version of loopNt!, communicates among processes at end of each Ns loop too update search space
@inline function paraloopNt!(current::State, result::Result, options::Options)
    current.f = result.fopt
    current.x .= result.xopt
    @inbounds for _ in 1:options.Nt
        # loopNs, collect results and adjust search space
        @everywhere workers() SimulatedAnnealing.loopNs!(current, result, options)
        collect_results!(current, result)
        adjust_step!(current, options)

        # Reset accepted points for search space
        current.n_accepted .= 0

        # Send updated search space and results to workers
        result.n_moves += options.N * options.Ns
        @everywhere workers() current, result = $current, $result
    end

    return nothing
end

# Complete Ns loops of loopN!
@inline function loopNs!(current::State, result::Result, options::Options)
    @inbounds for _ in 1:options.Ns
        loopN!(current, result, options)
    end

    return nothing
end

# Pertub each dimension of x independently and accept/reject new point
@inline function loopN!(current::State, result::Result, options::Options)
    @inbounds for i in 1:options.N
        # Perturb current x along dimension i and evaluate new point
        current.xtemp .= current.x
        perturb!(current.xtemp, current.v, i)

        # Check if new point in bounds - if not pick random point in bounds
        check_bounds!(current, options)

        # Evaluate function at new point
        current.ftemp = (2 * options.max - 1) * options.f(current.xtemp)

        # Decide whether to accept point or not
        decision(current) && update_state!(current, result, i)

        # Record if new global optimum
        current.f > result.fopt && update_result!(result, current)

        # Increment number of moves
        result.n_moves += 1

        # Check number of moves
        result.n_moves ≥ options.max_eval && (result.terminate = :maxeval)

    end

    return nothing
end

# Adjust step length based on number of accepted points
function adjust_step!(current::State, options::Options)
    current.n_accepted ./= options.Ns  # Ratio of accepted points
    current.v .= criteria.(current.v, current.n_accepted, options.c)

    return nothing
end

# Criteria aims to keep a 50% acceptance rate
function criteria(v, accepted, c)
    accepted < 0.0 || accepted > 1.0 ? error("Accepted points out of bounds $accepted") : nothing
    if accepted > 0.6
        v = v * (1 + c * (accepted - 0.6) / 0.4)
    elseif accepted < 0.4
        v = v / (1 + c * (0.4 - accepted) / 0.4)
    end

    return v
end

# Perturbs x by an random amount (-1.0 to 1.0) * step length for ith dimension
perturb!(x, v, i) = x[i] = x[i]  + (2 * rand() - 1) * v[i]

# Check new point in bounds, otherwise pick random point in bounds
function check_bounds!(current::State, options::Options)
    @inbounds for i in 1:options.N
        if current.xtemp[i] < options.lb[i] || current.xtemp[i] > options.ub[i]
            current.xtemp[i] = options.lb[i] + (options.ub[i] - options.lb[i]) * rand()
        end
    end

    return nothing
end

# Decide whether to accept new point or not
decision(current::State) = current.ftemp > current.f ? true : metropolis(current)

# Metropolis criterion: accept with probability 1 if f > fopt, or accept with probability = exp(f-fopt / T)
metropolis(current::State) = rand() < exp((current.ftemp - current.f) / current.T)

# Reduce temperature
reduce_temp!(current::State, options::Options) = current.T *= options.R

# Update current struct
function update_state!(current::State, result::Result, i)
    current.f = current.ftemp
    current.x .= current.xtemp
    current.n_accepted[i] += 1.0
    result.n_accepted += 1

    return nothing
end

function update_state!(current::State, trial::State)
    current.f = trial.f
    current.x .= trial.x
end

# Update result struct
function update_result!(result::Result, current::State)
    result.fopt = current.f
    result.xopt .= current.x
    result.n_optimums += 1

    return nothing
end

# Collect results from workers and update search space and result
function collect_results!(current::State, result::Result)
    for w in workers()
        trial = @getfrom w current

        # Update current if greater optimum found
        trial.f > current.f && update_state!(current, trial)

        # Update result if greater optimum found
        trial.f > result.fopt && update_result!(result, trial)

        # Update number of accepted point
        current.n_accepted .+= trial.n_accepted
        result.n_accepted += sum(trial.n_accepted)
    end

end

# Termination criteria for stopping optimisation
check_termination(result::Result) = result.terminate !== :none

# Check if improvement less than tolerance
check_improvement!(bests, result::Result, tol) = all((result.fopt .- bests) .< tol) && (result.terminate = :converge)

##### Printing functions
function print_intitial(current::State, options::Options, parallel)
    println("Starting values for simulated annealing:
            Initial f: $(current.f)
            Initial x: $(current.x)
            Initial step length: $(current.v)
            Lower bounds on x: $(options.lb)
            Upper bounds on x: $(options.ub)
            Initial temperature: $(current.T)
            Maximise function: $(options.max)
            N: $(options.N), Ns: $(options.Ns), Nt: $(options.Nt), tol: $(options.tol), max evaluations: $(options.max_eval)
            Parallel evaluation: $parallel")
end

function print_status(current::State, result::Result, options::Options)
    println("Result after temperature reduction:
             f: $(options.max ? result.fopt : -result.fopt)
             x: $(result.xopt)
             Step length: $(current.v),
             Moves: $(result.n_moves), accepted: $(result.n_accepted), optimums: $(result.n_optimums)
             ")
end

function print_termination(current::State, result::Result, options::Options)
    println("Simulated annealing reached termination criteria:
             Criteria: $(result.terminate)
             Optimum f found: $(result.fopt)
             Optimum x found: $(result.xopt)
             Final step length: $(current.v)
             Final temperature: $(current.T)
             Number of moves: $(result.n_moves), accepted: $(result.n_accepted), optimums: $(result.n_optimums)
             ")

end

##### Test function
# Multidimensional Rosenbrock, optimum at repeat([1], N)
f(x) = begin
    sleep(0.001)
    [100 * (x[i + 1] - x[i]^2)^2 + (1 - x[i])^2 for i ∈ 1:length(x) - 1] |> sum
end
# Initialise structs to optimise n dimension rosenbrock
function init(n)
    x₀ = rand(n)
    options = Options(func=f, N=n)
    current = State(f=f(x₀), x=x₀, v=options.ub - options.lb)
    result = Result(f(x₀), x₀)

    return current, result, options
end


text:
using Test, Joseki

@testset "Utilities" begin
    j1 = Dict("a" => 1, "b" => 2)
    j2 = Dict("a" => 1, "c" => 2)
    required_keys = ["a", "b"]
    @test has_all_required_keys(required_keys, j1)
    @test !has_all_required_keys(required_keys, j2)

    req = sample_request()
    j = Joseki.body_as_dict(req)
    @test j == Dict{String, Any}("x" => 5, "y" => 3)
end