text:
module PairDistribution

using StatsBase
using Interpolations

export rdf, rdf1d, rdf2p, rdfinter2d

"""
    distpbc(v1, v2, box)

Calculate a distance vector wrt accoring to periodic boundary conditions, defined in `box`.
"""
function distpbc(v1, v2, box)
    rij = v2 - v1
    for (k, l) in enumerate(box)
        hl = l/2
        if rij[k] > hl
            rij[k] = rij[k] - l
        elseif rij[k] < -hl
            rij[k] = rij[k] + l
        end
    end
    return rij
end

"""
    rdfhist(func, data::Array{T,3}, bins::Int...) -> Array{N,Int}

Calculate `N`D histogrm with `bins` bins using a function `func` for evaluating
values of an empirical distribution.
"""
function rdfhist(func, data::Array{T,3}, bins::Int...; debug=0,
                 minmax=fill((zero(T), maximum(data)/2),length(bins)) ) where T
    pnum, _, fnum = size(data)
    edges = tuple([range(mn, mx, length=b+1) for ((mn,mx),b) in zip(minmax,bins)]...)
    hist = Histogram(edges)
    Δs = map(step, hist.edges)
    for t in 1:fnum # t = 1
        debug>0 && println("Frame $t")
        for i in 1:pnum # i = 1
            ri = @view data[i,2:end,t]
            for j in i+1:pnum # j = 2
                rj = @view data[j,2:end,t]
                for v in func((Δs, ri, rj))
                    push!(hist, v)
                end
            end
        end
    end
    hist
end

"""
    rdf(func, data::Array{T,3}, box::Vector{T}, bins::Int...; kwargs...)

Calculate a `N`D radial distribution function (RDF) using a function `func` for constructing
an empirical distribution for RDF.
"""
function rdf(func, data::Array{T,3}, box::Vector{T}, bins::Int...;
            volfunc = (up,lw)->(up^3 - lw^3), kwargs...) where T
    hist = rdfhist(func, data, bins...; kwargs...)
    Δs = map(step, hist.edges)
    pnum, _, fnum = size(data)
    N = fnum * pnum
    ρ = pnum / prod(box) #  Density of the particles
    c = 4 / 3 * π * ρ

    idxs = CartesianIndices(tuple((1:i for i in bins)...))
    gr = zeros(T, bins)
    for ci in idxs, i in 1:length(ci)
        if i == 1
            lw = (ci.I[i] - 1)*Δs[i]
            up = lw + Δs[i]
            dVρ = c * volfunc(up, lw)
            gr[ci] = hist.weights[ci] / N / dVρ
        end
    end
    rs = [collect(r)[1:end-1].+dr for (r,dr) in zip(hist.edges, Δs)]
    (gr,rs)
end

"""
    rdf1d(data::Array{T,3}, box::Vector{T}, bins::Int...; kwargs...)

Calculate a 1D radial distribution function (RDF) using distance between particles.
"""
function rdf1d(data::Array{T,3}, box::Vector{T}, bins::Int...; kwargs...) where T
    rdf(data, box, bins...; kwargs...) do (_, ri, rj)
        rij = PairDistribution.distpbc(rj, ri, box)
        r2 = sum(abs2, rij)
        [sqrt(r2), sqrt(r2)]
    end
end

"""
    rdf2p(data::Array{T,3}, box::Vector{T}, bins::Int...; kwargs...)

Calculate a 2D radial distribution function (RDF) using distance between particles
and azimuthal angle.
"""
function rdf2p(data::Array{T,3}, box::Vector{T}, bins::Int...; kwargs...) where T
    rdf(data, box, bins...; kwargs...) do (_, ri, rj)
        rij = PairDistribution.distpbc(rj, ri, box)
        r2 = sum(abs2, rij)
        r = sqrt(r2)
        ϕ = atan(rij[2],rij[3]) # azimuthal angle
        ϕs = atan(-rij[2],rij[3]) # azimuthal angle
        [(r,ϕ),(r,ϕs)]
    end
end

"""
    rdfinter2d(gr, rs; bins=30, maxval=0.1, dist=(0.0,Inf))

Construct interpolated 2D RDF with dimensions `bins`x`bins`.

While constructing interpolation filter all value larger then `maxval`,
and which are not in the distance range `dist`.
"""
function rdfinter2d(gr, rs; bins=30, maxval=0.1, dist=(0.0,Inf))
    z = hcat(copy(gr), gr[:,1])
    # filter distances
    minidx = findfirst(d->d>=dist[1],rs[1])
    minidx = minidx === nothing ? 1 : minidx
    maxidx = findfirst(d->d>=dist[2],rs[1])
    maxidx = maxidx === nothing ? length(rs[1]) : maxidx
    dd = rs[1][minidx:maxidx]
    # convert angles to [0,2π] range
    rr = [0; rs[2].+π]
    itp = interpolate((dd, rr), z[minidx:maxidx,:], Gridded(Linear()))
    xx = range(extrema(rr)..., length=bins) |> collect
    yy = range(extrema(dd)..., length=bins) |> collect
    zz = [
        let v = itp(y,x)
            v > maxval ? 0.0 : v
        end
        for y in yy, x in xx]
    xx, yy, zz # dim(yy) x dim(xx)
end

end

max_stars_repo_name: wildart/PairDistribution

max_issues_repo_name: wildart/PairDistribution

max_forks_repo_name: wildart/PairDistribution


text:
using Distributed
addprocs(2, exeflags="--project")

@everywhere begin
include("../../src/Constructors.jl")
include("../../src/RemoteBufferChannel.jl")
include("../../src/SelfPlay.jl")
include("../../src/ReplayBuffer.jl")
include("../../src/Learning.jl")
include("game.jl")
include("params.jl")
end

env=TicTacToe()
training_step = RemoteChannel(()->Channel{Int}(1))
num_played_games = RemoteChannel(()->Channel{Int}(1))
num_played_steps = RemoteChannel(()->Channel{Int}(1))
num_reanalysed_games = RemoteChannel(()->Channel{Int}(1))
total_samples = RemoteChannel(()->Channel{Int}(1))
remote_NNs = RemoteChannel(()->Channel{NamedTuple{(:representation, :prediction, :dynamics), Tuple{Any, Any, Any}}}(1))
remote_buffer = RemoteChannel(()->BufferChannel()) #TODO Start saving Buffer to disk at periodic intervals after a certain stage in the training process

put!(remote_NNs, (representation=init_representation(hyper), prediction=init_prediction(hyper), dynamics=init_dynamics(hyper)))
put!(training_step, 0)
put!(num_played_games, 0)
put!(num_played_steps, 0)
put!(num_reanalysed_games, 0)
put!(total_samples, 0)

sp = @spawnat :any self_play!(env,
training_step,
num_played_games,
num_played_steps,
total_samples,
remote_NNs,
remote_buffer)

learn = @spawnat :any learning!(num_played_games,
training_step,
remote_NNs,
remote_buffer)

# include("./src/MuZero.jl")
# using .MuZero
# include("./games/tictactoe/game.jl")
# ttt=TicTacToe()
# include("./games/tictactoe/params.jl")
# MuZero.self_play(hyper, ttt, MuZero.progress, MuZero.buffer)

ext: jl


text:
#####
##### Divergence operators
#####

"""
    divᶜᶜᶜ(i, j, k, grid, u, v, w)

Calculates the divergence ∇·U of a vector field U = (u, v, w),

    1/V * [δxᶜᵃᵃ(Ax * u) + δxᵃᶜᵃ(Ay * v) + δzᵃᵃᶜ(Az * w)],

which will end up at the cell centers `ccc`.
"""
@inline function divᶜᶜᶜ(i, j, k, grid, u, v, w)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_ψᵃᵃᶠ, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_ψᵃᵃᶠ, v) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_ψᵃᵃᵃ, w))
end

"""
    div_xyᶜᶜᵃ(i, j, k, grid, u, v)

Calculates the 2D divergence ∂x u + ∂y v via

    1/Azᶜᶜᵃ * [δxᶜᵃᵃ(Δy * u) + δyᵃᶜᵃ(Δx * v)]

where `Azᶜᶜᵃ` is the area of the cell centered on (Center, Center, Any) --- a tracer cell,
`Δy` is the length of the cell centered on (Face, Center, Any) in `y` (a `u` cell),
and `Δx` is the length of the cell centered on (Center, Face, Any) in `x` (a `v` cell).
`div_xyᶜᶜᵃ` ends up at the location `cca`.
"""
@inline function div_xyᶜᶜᵃ(i, j, k, grid, u, v)
    return 1 / Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_uᶠᶜᶜ, u) +
                                      δyᵃᶜᵃ(i, j, k, grid, Ay_vᶜᶠᶜ, v))
end


text:
module MPITests

using MPI
using Test

#Sysimage
sysimage=nothing
if length(ARGS)==1
   @assert isfile(ARGS[1]) "$(ARGS[1]) must be a valid Julia sysimage file"
   sysimage=ARGS[1]
end

mpidir = @__DIR__
testdir = joinpath(mpidir,"..")
repodir = joinpath(testdir,"TestApp")
function run_driver(procs,file,sysimage)
  mpiexec() do cmd
    if sysimage!=nothing
       extra_args="-J$(sysimage)"
       run(`$cmd -n $procs $(Base.julia_cmd()) $(extra_args) --project=$repodir $(joinpath(mpidir,file))`)
    else
       run(`$cmd -n $procs $(Base.julia_cmd()) --project=$repodir $(joinpath(mpidir,file))`)
    end
    @test true
  end
end

run_driver(4,"runtests_np4.jl",sysimage)
run_driver(1,"runtests_np4.jl",sysimage) # Check that the degenerated case works

end # module


text:
"""
    module IPIL8

Пакет реализует исследование численные методы решения обратных задач для нелинейных сингулярно
возмущённых уравнений типа реакция-диффузия-адвекция с данными о положении фронта реакции.

Автор исходного кода пакета:
    - [Андрей Борзунов](mailto:aborzunov@physics.msu.ru), Кафедра математики физического факультета МГУ им. Ломоносова.

"""
module IPIL8

using LinearAlgebra
using Plots
using ProgressMeter
using Dierckx
using LaTeXStrings
using Missings
using Printf

include("direct.jl")
include("utilities.jl")
include("plotting.jl")
include("data_generation.jl")

end


text:
module TestReparametrization
include("preamble.jl")

using StaticArrays: SVector

using Bifurcations.Examples.Reparametrization: reparametrize, forward, backward
using Bifurcations.Examples.Bautin

ode0 = reparametrize(Bautin.ode)
ode1 = reparametrize(Bautin.ode; seed=1, extra=SVector(-1.0))
ode2 = reparametrize(Bautin.ode; seed=2, extra=SVector(-1.0, +1.0))

@testset begin
    x1 = [0.1, 0.1]
    y1 = forward(x1, ode0.p)
    x2 = backward(y1, ode0.p)
    y2 = forward(x2, ode0.p)
    @test x1 ≈ x2
    @test y1 ≈ y2

    x1 = [0.1, 0.1, 0.1]
    y1 = forward(x1, ode1.p)
    x2 = backward(y1, ode1.p)
    y2 = forward(x2, ode1.p)
    @test x1 ≈ x2
    @test y1 ≈ y2

    x1 = [0.1, 0.1, 0.1, 0.1]
    y1 = forward(x1, ode2.p)
    x2 = backward(y1, ode2.p)
    y2 = forward(x2, ode2.p)
    @test x1 ≈ x2
    @test y1 ≈ y2
end

@testset "Random test dim=$dim" for (dim, ode) in [(2, ode0),
                                                   (3, ode1),
                                                   (4, ode2)]
    for seed in 1:10
        rng = MersenneTwister(seed)
        x1 = randn(rng, dim)
        y1 = forward(x1, ode.p)
        x2 = backward(y1, ode.p)
        y2 = forward(x2, ode.p)
        @test x1 ≈ x2
        @test y1 ≈ y2
    end
end

end  # module


text:
#### Balance Law Interface
"""
    abstract type BalanceLaw end

An abstract type representing a PDE balance law of the form:

```math
\\frac{dq}{dt} = \\nabla \\cdot F_1(q, a, t) + \\nabla \\cdot F_2(q, \\nabla g, h, a, t) + S(q, \\nabla g, a, t)
```
where:
- ``q`` is the prognostic state,
- ``a`` is the auxiliary state,
- ``g = G(q, a, t)`` is the gradient state (variables of which we compute the
  gradient),
- ``h`` is the hyperdiffusive state.

Subtypes of `BalanceLaw` should define the following interfaces:
- [`vars_state`](@ref) to define the prognostic, auxiliary and intermediate
  variables.
- [`flux_first_order!`](@ref) to compute ``F_1``
- [`flux_second_order!`](@ref) to compute ``F_2``
- [`source!`](@ref) to compute ``S``

If `vars(bl, ::GradientFlux, FT)` is non-empty, then the following should be
defined:
- [`compute_gradient_argument!`](@ref) to compute ``G``
- [`compute_gradient_flux!`](@ref) is a linear transformation of ``\\nabla g``

If `vars(bl, ::Hyperdiffusive, FT)` is non-empty, then the following should be
defined:
- [`transform_post_gradient_laplacian!`](@ref)

Additional functions:
- [`wavespeed`](@ref) if using the Rusanov numerical flux.
- [`boundary_state!`](@ref) if using non-periodic boundary conditions.
"""
abstract type BalanceLaw end

"""
    BalanceLaws.vars_state(::BL, ::AbstractStateType, FT)

Defines the state variables of a [`BalanceLaw`](@ref) subtype `BL` with floating
point type `FT`.

For each [`AbstractStateType`](@ref), this should return a `NamedTuple` type,
with element type either `FT`, an `SArray` with element type `FT` or another
`NamedTuple` satisfying the same property.

For convenience, we recommend using the [`VariableTemplates.@vars`](@ref) macro.

# Example
```julia
struct MyBalanceLaw <: BalanceLaw end

BalanceLaws.vars_state(::MyBalanceLaw, ::Prognostic, FT) =
    @vars(x::FT, y::SVector{3, FT})
BalanceLaws.vars_state(::MyBalanceLaw, ::Auxiliary, FT) =
    @vars(components::@vars(a::FT, b::FT))
```
"""
function vars_state end

# Fallback: no variables
vars_state(::BalanceLaw, ::AbstractStateType, FT) = @vars()

"""
    init_state_prognostic!(
        ::BL,
        state_prognostic::Vars,
        state_auxiliary::Vars,
        localgeo,
        args...,
    )

Sets the initial state of the prognostic variables `state_prognostic` at each
node for a [`BalanceLaw`](@ref) subtype `BL`.
"""
function init_state_prognostic! end

# TODO: make these functions consistent with init_state_prognostic!
"""
    nodal_init_state_auxiliary!(::BL, state_auxiliary, state_temporary, geom)

Sets the initial state of the auxiliary variables `state_auxiliary` at each
node for a [`BalanceLaw`](@ref) subtype `BL`.

See also [`init_state_auxiliary!`](@ref).
"""
function nodal_init_state_auxiliary!(m::BalanceLaw, aux, tmp, geom) end

"""
    init_state_auxiliary!(
        ::BL,
        statearray_auxiliary,
        geom::LocalGeometry,
    )

Sets the initial state of the auxiliary variables `state_auxiliary` at each node
for a [`BalanceLaw`](@ref) subtype `BL`. By default this calls
[`nodal_init_state_auxiliary!`](@ref).
"""
function init_state_auxiliary!(
    balance_law::BalanceLaw,
    statearray_auxiliary,
    grid,
    direction,
)
    init_state_auxiliary!(
        balance_law,
        nodal_init_state_auxiliary!,
        statearray_auxiliary,
        grid,
        direction,
    )
end

"""
    flux_first_order!(
        ::BL,
        flux::Grad,
        state_prognostic::Vars,
        state_auxiliary::Vars,
        t::Real,
        direction
    )

Sets the first-order (hyperbolic) `flux` terms for a [`BalanceLaw`](@ref) subtype `BL`.
"""
function flux_first_order! end

"""
    flux_second_order!(
        ::BL,
        flux::Grad,
        state_prognostic::Vars,
        state_gradient_flux::Vars,
        hyperdiffusive::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Sets second-order (parabolic) `flux` terms for a [`BalanceLaw`](@ref) subtype `BL`.
"""
function flux_second_order! end

"""
    source!(
        ::BL,
        source::Vars,
        state_prognostic::Vars,
        diffusive::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Compute non-conservative source terms for a [`BalanceLaw`](@ref) subtype `BL`.
"""
function source! end

"""
    compute_gradient_argument!(
        ::BL,
        transformstate::Vars,
        state_prognostic::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Transformation of state variables `state_prognostic` to variables
`transformstate` of which gradients are computed for a [`BalanceLaw`](@ref)
subtype `BL`.
"""
function compute_gradient_argument! end

function compute_gradient_argument!(
    balance_law::BalanceLaw,
    transformstate::AbstractArray,
    state_prognostic::AbstractArray,
    state_auxiliary::AbstractArray,
    t,
)
    FT = eltype(transformstate)
    compute_gradient_argument!(
        balance_law,
        Vars{vars_state(balance_law, Gradient(), FT)}(transformstate),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary),
        t,
    )
end

"""
    compute_gradient_flux!(
        ::BL,
        state_gradient_flux::Vars,
        ∇transformstate::Grad,
        state_prognostic::Vars,
        state_auxiliary::Vars,
        t::Real
    )

Transformation of gradients to the diffusive variables for a
[`BalanceLaw`](@ref) subtype `BL`. This should be a linear function of
`∇transformstate`
"""
function compute_gradient_flux! end

function compute_gradient_flux!(
    balance_law::BalanceLaw,
    state_gradient_flux::AbstractArray,
    ∇transformstate::AbstractArray,
    state_prognostic::AbstractArray,
    state_auxiliary::AbstractArray,
    t,
)
    FT = eltype(state_gradient_flux)
    compute_gradient_flux!(
        balance_law,
        Vars{vars_state(balance_law, GradientFlux(), FT)}(state_gradient_flux),
        Grad{vars_state(balance_law, Gradient(), FT)}(∇transformstate),
        Vars{vars_state(balance_law, Prognostic(), FT)}(state_prognostic),
        Vars{vars_state(balance_law, Auxiliary(), FT)}(state_auxiliary),
        t,
    )
end

"""
    transform_post_gradient_laplacian!(
        ::BL,
        Qhypervisc_div::Vars,
        ∇Δtransformstate::Grad,
        state_auxiliary::Vars,
        t::Real
    )

Transformation of Laplacian gradients to the hyperdiffusive variables for a
[`BalanceLaw`](@ref) subtype `BL`.
"""
function transform_post_gradient_laplacian! end

"""
    wavespeed(
        ::BL,
        n⁻,
        state_prognostic::Vars,
        state_auxiliary::Vars,
        t::Real,
        direction
    )

Wavespeed in the direction `n⁻` for a [`BalanceLaw`](@ref) subtype `BL`. This is
required to be defined if using a `RusanovNumericalFlux` numerical flux.
"""
function wavespeed end

"""
    update_auxiliary_state!(
        dg::DGModel,
        m::BalanceLaw,
        statearray_aux,
        t::Real,
        elems::UnitRange,
        [diffusive=false]
    )

Hook to update the auxiliary state variables before calling any other functions.

By default, this calls [`nodal_update_auxiliary_state!`](@ref) at each node.

If `diffusive=true`, then `state_gradflux` is also passed to
`nodal_update_auxiliary_state!`.
"""
function update_auxiliary_state!(
    dg,
    balance_law::BalanceLaw,
    state_prognostic,
    t,
    elems,
    diffusive = false,
)
    update_auxiliary_state!(
        nodal_update_auxiliary_state!,
        dg,
        balance_law,
        state_prognostic,
        t,
        elems;
        diffusive = diffusive,
    )
end

"""
    nodal_update_auxiliary_state!(::BL, state_prognostic, state_auxiliary, [state_gradflux,] t)

Update the auxiliary state variables at each node for a [`BalanceLaw`](@ref)
subtype `BL`. By default it does nothing.

Called by [`update_auxiliary_state!`](@ref).
"""
function nodal_update_auxiliary_state!(args...)
    nothing
end

"""
    update_auxiliary_state_gradient!(
        dg::DGModel,
        m::BalanceLaw,
        statearray_aux,
        t::Real,
        elems::UnitRange,
        [diffusive=false]
    )

Hook to update the auxiliary state variables after the gradient computation.

By default, this calls nothing.

If `diffusive=true`, then `state_gradflux` is also passed to
`nodal_update_auxiliary_state!`.
"""
function update_auxiliary_state_gradient! end

# upward integrals
"""
    integral_load_auxiliary_state!(::BL, integrand, state_prognostic, state_aux)

Specify variables `integrand` which will have their upward integrals computed.

See also [`UpwardIntegrals`](@ref)
"""
function integral_load_auxiliary_state! end

"""
    integral_set_auxiliary_state!(::BL, state_aux, integral)

Update auxiliary variables based on the upward integral `integral` defined in
[`integral_load_auxiliary_state!`](@ref).
"""
function integral_set_auxiliary_state! end

"""
    indefinite_stack_integral!

Compute indefinite integral along stack.
"""
function indefinite_stack_integral! end

# downward integrals
"""
    reverse_integral_load_auxiliary_state!(::BL, integrand, state_prognostic, state_aux)

Specify variables `integrand` which will have their downward integrals computed.

See also [`DownwardIntegrals`](@ref)
"""
function reverse_integral_load_auxiliary_state! end

"""
    reverse_integral_set_auxiliary_state!(::BL, state_aux, integral)

Update auxiliary variables based on the downward integral `integral` defined in
[`reverse_integral_load_auxiliary_state!`](@ref).
"""
function reverse_integral_set_auxiliary_state! end

"""
    reverse_indefinite_stack_integral!

Compute reverse indefinite integral along stack.
"""
function reverse_indefinite_stack_integral! end

"""
    function state_to_entropy_variables!(
        balancelaw::BalanceLaw,
        state_entropy::Vars,
        state_prognostic::Vars,
        state_auxiliary::Vars,
    )
Compute entropy variables `state_entropy` from a given `state_prognostic` and
`state_auxiliary`
"""
function state_to_entropy_variables!(
    balancelaw::BalanceLaw,
    entropy::AbstractArray{FT},
    state::AbstractArray{FT},
    aux::AbstractArray{FT},
) where {FT}
    state_to_entropy_variables!(
        balancelaw,
        Vars{vars_state(balancelaw, Entropy(), FT)}(entropy),
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux),
    )
end

"""
    function entropy_variables_to_state!(
        balancelaw::BalanceLaw,
        state_prognostic::Vars,
        state_auxiliary::Vars,
        state_entropy::Vars,
    )
Compute `state_prognostic` and relevant componennts of `state_auxiliary` from
the entropy variables `state_entropy`
"""
function entropy_variables_to_state!(
    balancelaw::BalanceLaw,
    state::AbstractArray{FT},
    aux::AbstractArray{FT},
    entropy::AbstractArray{FT},
) where {FT}
    entropy_variables_to_state!(
        balancelaw,
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux),
        Vars{vars_state(balancelaw, Entropy(), FT)}(entropy),
    )
end

"""
    function state_to_entropy(
        balancelaw::BalanceLaw,
        state_prognostic::Vars,
        state_auxiliary::Vars,
    )
Compute entropy from a given `state_prognostic` and `state_auxiliary`
"""
function state_to_entropy(
    balancelaw::BalanceLaw,
    state::AbstractArray{FT},
    aux::AbstractArray{FT},
) where {FT}
    state_to_entropy(
        balancelaw,
        Vars{vars_state(balancelaw, Prognostic(), FT)}(state),
        Vars{vars_state(balancelaw, Auxiliary(), FT)}(aux),
    )
end

# Internal methods
number_states(m::BalanceLaw, st::AbstractStateType, FT = Int) =
    varsize(vars_state(m, st, FT))

### split explicit functions
function initialize_states! end
function tendency_from_slow_to_fast! end
function cummulate_fast_solution! end
function reconcile_from_fast_to_slow! end

parameter_set(balance_law) = balance_law.param_set

"""
    sub_model(::BalanceLaw, ::Type{AbstractSubModel})

Returns a tuple of balance law properties
whose supertypes are `AbstractSubModel`.
"""
function sub_model end


text:
module BinaryPlatformTests
import ..Pkg # ensure we are using the correct Pkg

using Test, Pkg.BinaryPlatforms
import Pkg.BinaryPlatforms: platform_name

# The platform we're running on
const platform = @inferred Platform platform_key_abi()

@testset "PlatformNames" begin
    # Ensure the platform type constructors are well behaved
    @testset "Platform constructors" begin
        @test_throws ArgumentError Linux(:not_a_platform)
        @test_throws ArgumentError Linux(:x86_64; libc=:crazy_libc)
        @test_throws ArgumentError Linux(:x86_64; libc=:glibc, call_abi=:crazy_abi)
        @test_throws ArgumentError Linux(:x86_64; libc=:glibc, call_abi=:eabihf)
        @test_throws ArgumentError Linux(:arm)
        @test_throws ArgumentError Linux(:armv7l; libc=:glibc, call_abi=:kekeke)
        @test_throws ArgumentError MacOS(:i686)
        @test_throws ArgumentError MacOS(:x86_64; libc=:glibc)
        @test_throws ArgumentError MacOS(:x86_64; call_abi=:eabihf)
        @test_throws ArgumentError Windows(:armv7l)
        @test_throws ArgumentError Windows(:x86_64; libc=:glibc)
        @test_throws ArgumentError Windows(:x86_64; call_abi=:eabihf)
        @test_throws ArgumentError FreeBSD(:not_a_platform)
        @test_throws ArgumentError FreeBSD(:x86_64; libc=:crazy_libc)
        @test_throws ArgumentError FreeBSD(:x86_64; call_abi=:crazy_abi)
        @test_throws ArgumentError FreeBSD(:x86_64; call_abi=:eabihf)

        @test_throws ArgumentError CompilerABI(;libgfortran_version=v"2")
        @test_throws ArgumentError CompilerABI(;libgfortran_version=v"99")
        @test_throws ArgumentError CompilerABI(;libstdcxx_version=v"3.3")
        @test_throws ArgumentError CompilerABI(;libstdcxx_version=v"3.5")
        @test_throws ArgumentError CompilerABI(;cxxstring_abi=:wut)

        # Test copy constructor
        cabi = CompilerABI(;
            libgfortran_version=v"3",
            libstdcxx_version=v"3.4.18",
            cxxstring_abi=:cxx03,
        )
        cabi2 = CompilerABI(cabi; cxxstring_abi=:cxx11)
        @test libgfortran_version(cabi) == libgfortran_version(cabi2)
        @test libstdcxx_version(cabi) == libstdcxx_version(cabi2)
        @test cxxstring_abi(cabi) != cxxstring_abi(cabi2)

        # Explicitly test that we can pass arguments to UnknownPlatform,
        # and it doesn't do anything.
        @test UnknownPlatform(:riscv; libc=:fuschia_libc) == UnknownPlatform()
    end

    @testset "Platform properties" begin
        # Test that we can get the name of various platforms
        for T in (Linux, MacOS, Windows, FreeBSD, UnknownPlatform)
            @test endswith(string(T), platform_name(T(:x86_64)))
        end

        # Test that we can get the arch of various platforms
        @test arch(Linux(:aarch64; libc=:musl)) == :aarch64
        @test arch(Windows(:i686)) == :i686
        @test arch(FreeBSD(:amd64)) == :x86_64
        @test arch(FreeBSD(:i386)) == :i686
        @test arch(UnknownPlatform(:ppc64le)) == nothing

        # Test that our platform_dlext stuff works
        @test platform_dlext(Linux(:x86_64)) == platform_dlext(Linux(:i686))
        @test platform_dlext(Windows(:x86_64)) == platform_dlext(Windows(:i686))
        @test platform_dlext(MacOS()) != platform_dlext(Linux(:armv7l))
        @test platform_dlext(FreeBSD(:x86_64)) == platform_dlext(Linux(:x86_64))
        @test platform_dlext(UnknownPlatform()) == "unknown"
        @test platform_dlext() == platform_dlext(platform)

        @test wordsize(Linux(:i686)) == wordsize(Linux(:armv7l)) == 32
        @test wordsize(MacOS()) == wordsize(Linux(:aarch64)) == 64
        @test wordsize(FreeBSD(:x86_64)) == wordsize(Linux(:powerpc64le)) == 64
        @test wordsize(UnknownPlatform(:x86_64)) == 0

        @test call_abi(Linux(:x86_64)) == nothing
        @test call_abi(Linux(:armv6l)) == :eabihf
        @test call_abi(Linux(:armv7l; call_abi=:eabihf)) == :eabihf
        @test call_abi(UnknownPlatform(;call_abi=:eabihf)) == nothing

        @test triplet(Windows(:i686)) == "i686-w64-mingw32"
        @test triplet(Linux(:x86_64; libc=:musl)) == "x86_64-linux-musl"
        @test triplet(Linux(:armv7l; libc=:musl)) == "armv7l-linux-musleabihf"
        @test triplet(Linux(:armv6l; libc=:musl, call_abi=:eabihf)) == "armv6l-linux-musleabihf"
        @test triplet(Linux(:x86_64)) == "x86_64-linux-gnu"
        @test triplet(Linux(:armv6l)) == "armv6l-linux-gnueabihf"
        @test triplet(MacOS()) == "x86_64-apple-darwin14"
        @test triplet(FreeBSD(:x86_64)) == "x86_64-unknown-freebsd11.1"
        @test triplet(FreeBSD(:i686)) == "i686-unknown-freebsd11.1"
        @test triplet(UnknownPlatform()) == "unknown-unknown-unknown"

        @test repr(Windows(:x86_64)) == "Windows(:x86_64)"
        @test repr(Linux(:x86_64; libc=:glibc, call_abi=nothing)) == "Linux(:x86_64, libc=:glibc)"
        @test repr(Linux(:armv7l; libc=:musl, call_abi=:eabihf)) == "Linux(:armv7l, libc=:musl, call_abi=:eabihf)"
        @test repr(MacOS()) == "MacOS(:x86_64)"
        @test repr(MacOS(compiler_abi=CompilerABI(cxxstring_abi=:cxx11))) == "MacOS(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx11))"
        @test repr(CompilerABI(;libgfortran_version=v"4", libstdcxx_version=v"3.4.24", cxxstring_abi=:cxx11)) == "CompilerABI(libgfortran_version=v\"4.0.0\", libstdcxx_version=v\"3.4.24\", cxxstring_abi=:cxx11)"
    end

    @testset "Valid DL paths" begin
        # Test some valid dynamic library paths
        @test valid_dl_path("libfoo.so.1.2.3", Linux(:x86_64))
        @test valid_dl_path("libfoo.1.2.3.so", Linux(:x86_64))
        @test valid_dl_path("libfoo-1.2.3.dll", Windows(:x86_64))
        @test valid_dl_path("libfoo.1.2.3.dylib", MacOS())
        @test !valid_dl_path("libfoo.dylib", Linux(:x86_64))
        @test !valid_dl_path("libfoo.so", Windows(:x86_64))
        @test !valid_dl_path("libfoo.dll", MacOS())
        @test !valid_dl_path("libfoo.so.1.2.3.", Linux(:x86_64))
        @test !valid_dl_path("libfoo.so.1.2a.3", Linux(:x86_64))
    end

    @testset "platform_key_abi parsing" begin
        # Make sure the platform_key_abi() with explicit triplet works
        @test @inferred(Platform, platform_key_abi("x86_64-linux-gnu")) == Linux(:x86_64)
        @test platform_key_abi("x86_64-linux-musl") == Linux(:x86_64, libc=:musl)
        @test platform_key_abi("i686-unknown-linux-gnu") == Linux(:i686)
        @test platform_key_abi("x86_64-apple-darwin14") == MacOS()
        @test platform_key_abi("x86_64-apple-darwin17.0.0") == MacOS()
        @test platform_key_abi("armv7l-pc-linux-gnueabihf") == Linux(:armv7l)
        @test platform_key_abi("armv7l-linux-musleabihf") == Linux(:armv7l, libc=:musl)
        @test platform_key_abi("armv6l-linux-gnueabihf") == Linux(:armv6l)
        # Test that the short name "arm" goes to `armv7l`
        @test platform_key_abi("arm-linux-gnueabihf") == Linux(:armv7l)
        @test platform_key_abi("aarch64-unknown-linux-gnu") == Linux(:aarch64)
        @test platform_key_abi("powerpc64le-linux-gnu") == Linux(:powerpc64le)
        @test platform_key_abi("ppc64le-linux-gnu") == Linux(:powerpc64le)
        @test platform_key_abi("x86_64-w64-mingw32") == Windows(:x86_64)
        @test platform_key_abi("i686-w64-mingw32") == Windows(:i686)
        @test platform_key_abi("x86_64-unknown-freebsd11.1") == FreeBSD(:x86_64)
        @test platform_key_abi("i686-unknown-freebsd11.1") == FreeBSD(:i686)
        @test platform_key_abi("amd64-unknown-freebsd12.0") == FreeBSD(:x86_64)
        @test platform_key_abi("i386-unknown-freebsd10.3") == FreeBSD(:i686)
        @test platform_key_abi("arm64-apple-darwin20.0") == MacOS(:aarch64)

        # Test inclusion of ABI stuff, both old-style and new-style
        @test platform_key_abi("x86_64-linux-gnu-gcc7") == Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4"))
        @test platform_key_abi("x86_64-linux-gnu-gcc4-cxx11") == Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"3", cxxstring_abi=:cxx11))
        @test platform_key_abi("x86_64-linux-gnu-cxx11") == Linux(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx11))
        @test platform_key_abi("x86_64-linux-gnu-libgfortran3-cxx03") == Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"3", cxxstring_abi=:cxx03))
        @test platform_key_abi("x86_64-linux-gnu-libstdcxx26") == Linux(:x86_64, compiler_abi=CompilerABI(libstdcxx_version=v"3.4.26"))

        # Make sure some of these things are rejected
        function test_bad_platform(p_str)
            @test_logs (:warn, r"not an officially supported platform") begin
                @test platform_key_abi(p_str) == UnknownPlatform()
            end
        end

        test_bad_platform("totally FREEFORM text!!1!!!1!")
        test_bad_platform("invalid-triplet-here")
        test_bad_platform("aarch64-linux-gnueabihf")
        test_bad_platform("x86_64-w32-mingw64")
    end

    @testset "platforms_match()" begin
        # Just do a quick combinatorial sweep for completeness' sake for platform matching
        for libgfortran_version in (nothing, v"3", v"5"),
            libstdcxx_version in (nothing, v"3.4.18", v"3.4.26"),
            cxxstring_abi in (nothing, :cxx03, :cxx11)

            cabi = CompilerABI(;
                libgfortran_version=libgfortran_version,
                libstdcxx_version=libstdcxx_version,
                cxxstring_abi=cxxstring_abi,
            )
            @test platforms_match(Linux(:x86_64), Linux(:x86_64, compiler_abi=cabi))
            @test platforms_match(Linux(:x86_64, compiler_abi=cabi), Linux(:x86_64))

            # Also test auto-string-parsing
            @test platforms_match(triplet(Linux(:x86_64)), Linux(:x86_64, compiler_abi=cabi))
            @test platforms_match(Linux(:x86_64), triplet(Linux(:x86_64, compiler_abi=cabi)))
        end

        # Ensure many of these things do NOT match
        @test !platforms_match(Linux(:x86_64), Linux(:i686))
        @test !platforms_match(Linux(:x86_64), Windows(:x86_64))
        @test !platforms_match(Linux(:x86_64), MacOS())
        @test !platforms_match(Linux(:x86_64), UnknownPlatform())

        # Make some explicitly non-matching cabi's
        base_cabi = CompilerABI(;
            libgfortran_version=v"5",
            cxxstring_abi=:cxx11,
        )
        for arch in (:x86_64, :i686, :aarch64, :armv6l, :armv7l),
            cabi in (
                CompilerABI(libgfortran_version=v"3"),
                CompilerABI(cxxstring_abi=:cxx03),
                CompilerABI(libgfortran_version=v"4", cxxstring_abi=:cxx11),
                CompilerABI(libgfortran_version=v"3", cxxstring_abi=:cxx03),
            )

            @test !platforms_match(Linux(arch, compiler_abi=base_cabi), Linux(arch, compiler_abi=cabi))
        end
    end

    @testset "DL name/version parsing" begin
        # Make sure our version parsing code is working
        @test parse_dl_name_version("libgfortran.dll", Windows(:x86_64)) == ("libgfortran", nothing)
        @test parse_dl_name_version("libgfortran-3.dll", Windows(:x86_64)) == ("libgfortran", v"3")
        @test parse_dl_name_version("libgfortran-3.4.dll", Windows(:x86_64)) == ("libgfortran", v"3.4")
        @test parse_dl_name_version("libgfortran-3.4a.dll", Windows(:x86_64)) == ("libgfortran-3.4a", nothing)
        @test_throws ArgumentError parse_dl_name_version("libgfortran", Windows(:x86_64))
        @test parse_dl_name_version("libgfortran.dylib", MacOS(:x86_64)) == ("libgfortran", nothing)
        @test parse_dl_name_version("libgfortran.3.dylib", MacOS(:x86_64)) == ("libgfortran", v"3")
        @test parse_dl_name_version("libgfortran.3.4.dylib", MacOS(:x86_64)) == ("libgfortran", v"3.4")
        @test parse_dl_name_version("libgfortran.3.4a.dylib", MacOS(:x86_64)) == ("libgfortran.3.4a", nothing)
        @test_throws ArgumentError parse_dl_name_version("libgfortran", MacOS(:x86_64))
        @test parse_dl_name_version("libgfortran.so", Linux(:x86_64)) == ("libgfortran", nothing)
        @test parse_dl_name_version("libgfortran.so.3", Linux(:x86_64)) == ("libgfortran", v"3")
        @test parse_dl_name_version("libgfortran.so.3.4", Linux(:x86_64)) == ("libgfortran", v"3.4")
        @test_throws ArgumentError parse_dl_name_version("libgfortran.so.3.4a", Linux(:x86_64))
        @test_throws ArgumentError parse_dl_name_version("libgfortran", Linux(:x86_64))
    end

    @testset "Sys.is* overloading" begin
        # Test that we can indeed ask if something is linux or windows, etc...
        @test Sys.islinux(Linux(:aarch64))
        @test !Sys.islinux(Windows(:x86_64))
        @test Sys.iswindows(Windows(:i686))
        @test !Sys.iswindows(Linux(:x86_64))
        @test Sys.isapple(MacOS())
        @test !Sys.isapple(Linux(:powerpc64le))
        @test Sys.isbsd(MacOS())
        @test Sys.isbsd(FreeBSD(:x86_64))
        @test !Sys.isbsd(Linux(:powerpc64le; libc=:musl))
    end

    @testset "Compiler ABI detection" begin
        @test detect_libgfortran_version("libgfortran.so.5", Linux(:x86_64)) == v"5"
        @test detect_libgfortran_version("libgfortran.4.dylib", MacOS()) == v"4"
        @test detect_libgfortran_version("libgfortran-3.dll", Windows(:x86_64)) == v"3"
        @test_logs (:warn, r"Unable to determine libgfortran version") begin
            @test detect_libgfortran_version("blah.so", Linux(:aarch64)) == nothing
        end

        # Let's check and ensure that we can autodetect the currently-running Julia process
        @test detect_libgfortran_version() != nothing

        # We run these to get coverage, but we can't test anything, because we could be built
        # with `clang`, which wouldn't have any `libstdc++` constraints at all
        detect_libstdcxx_version()
        detect_cxxstring_abi()
    end
end

@testset "select_platform" begin
    platforms = Dict(
        # Typical binning test
        Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"3")) => "linux4",
        Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4")) => "linux7",
        Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"5")) => "linux8",

        # Ambiguity test
        Linux(:aarch64, compiler_abi=CompilerABI(libgfortran_version=v"3")) => "linux4",
        Linux(:aarch64, compiler_abi=CompilerABI(libgfortran_version=v"3", libstdcxx_version=v"3.4.18")) => "linux5",

        MacOS(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"3")) => "mac4",
        Windows(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx11)) => "win",
    )

    @test select_platform(platforms, Linux(:x86_64)) == "linux8"
    @test select_platform(platforms, Linux(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4"))) == "linux7"

    # Ambiguity test
    @test select_platform(platforms, Linux(:aarch64)) == "linux5"
    @test select_platform(platforms, Linux(:aarch64; compiler_abi=CompilerABI(libgfortran_version=v"3"))) == "linux5"
    @test select_platform(platforms, Linux(:aarch64; compiler_abi=CompilerABI(libgfortran_version=v"4"))) == nothing

    @test select_platform(platforms, MacOS(:x86_64)) == "mac4"
    @test select_platform(platforms, MacOS(:x86_64, compiler_abi=CompilerABI(libgfortran_version=v"4"))) == nothing

    @test select_platform(platforms, Windows(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx11))) == "win"
    @test select_platform(platforms, Windows(:x86_64, compiler_abi=CompilerABI(cxxstring_abi=:cxx03))) == nothing

    # Poor little guy
    @test select_platform(platforms, FreeBSD(:x86_64)) == nothing
end

end # module


text:
function retrieve(c::Client, ::Type{VoiceRegion}, g::AbstractGuild)
    return get_guild_voice_regions(c, g.id)
end
retrieve(c::Client, ::Type{VoiceRegion}) = list_voice_regions(c)


text:
function explicityEnumeration!(data::alpheusData)

    options = data.options
    input = data.input
    preprocessor = data.preprocessor
    initialiseEnumerateOutput!(data)

    for k in 2:input.J-1
        ### enumerate ###
        z = zeros(J)
        z[k] = 1

        ### initialize model ###
        if options.solver == bonmin
            model = Model(() -> AmplNLWriter.Optimizer(Bonmin_jll.amplexe));
        elseif options.solver == couenne
            model = Model(() -> AmplNLWriter.Optimizer(Couenne_jll.amplexe));
        elseif options.solver == ipopt
            model = Model(Ipopt.Optimizer)
        end

        ### set parameters ###
        for (key, value) in options.parameters
            set_optimizer_attribute(model, key, value)
        end

        ### create variables ###
        @variable(model, 0.0 <= H[1:J] <= deepest, start = 1.0);
        @variable(model, 0.0 <= V[1:J] <= g*deepest, start = 1.0);

        ### create expressions ###
        @NLexpression(model, W[j in 1:J], z[j]*input.N_turbine*(V[j]/input.Vmax^2)^(3/2)*preprocessor.K);
        @NLexpression(model, A[j in 1:J], input.B*H[j]);
        @NLexpression(model, R[j in 1:J], A[j]/(input.B+2*H[j]));
        @NLexpression(model, S[j in 1:J], V[j]*(input.n0^2+input.nt^2*z[j]+2*input.n0*input.nt*z[j])/(R[j]^(4/3)));
        @NLexpression(model,SM[j in 1:J-1],(S[j] + S[j+1])/2);

        ### create constraint ###
        if !options.upstream
            @constraint(model, H[1] <= input.Hu);
        end
        if !options.downstream
            @constraint(model, H[J] >= input.Hd);
        end

        @constraint(model, [j in 1:J], V[j] >= input.Vmin*z[j]);
        @constraint(model, [j in 1:J], V[j] <= (input.Vmax - g*deepest)*z[j] + g*deepest);

        @constraint(model, [j in 1:J], H[j] >= input.Gt*z[j]);
        @constraint(model, [j in 1:J-1], -(1-z[j])*input.M*H[j]-z[j]*deepest <= H[j+1] - H[j]);
        @constraint(model, [j in 1:J-1], H[j+1] - H[j] <= (1-z[j])*input.M*H[j]+z[j]*deepest);
        @constraint(model, [j in 1:J], V[j] <= g*H[j]);
        @NLconstraint(model,[j in 1:J], V[j]*A[j]^2 == input.Q^2);
        @NLconstraint(model, [j in 1:J-1], V[j]/(2*g) + H[j] + preprocessor.E[j] == V[j+1]/(2*g) + H[j+1] + preprocessor.E[j+1] + SM[j]*preprocessor.Δx + W[j]);

        ### create objective ###

        @NLobjective(model, Max,sum(W[j] for j in 1:J) -(input.Hu-H[1])*input.ηu*options.upstream + (H[J]-input.Hd)*input.ηd*options.downstream)

        optimize!(model)
        status = termination_status(model)
        @show status
        if status in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
            fillOutputEnumeration!(model,data,k)
        end

    end
    nothing
end

function initialiseEnumerateOutput!(data::alpheusData)
    output = data.output
    input = data.input

    output.W_enum = zeros(input.J)
    output.H_enum = zeros(input.J,input.J)
    output.U_enum = zeros(input.J,input.J)
    output.z_enum = zeros(input.J,input.J)
    nothing
end

function fillOutputEnumeration!(model::Model,data::alpheusData,k::Int)
    output = data.output

    output.W_enum[k] = value.(model[:W])[k];
    output.H_enum[k,:] = value.(model[:H]);
    V = value.(model[:V]);
    output.U_enum[k,:] = sqrt.(V)
    output.z_enum[k,k] = 1
    nothing
end

max_stars_repo_name: Thiago-NovaesB/alpheus

max_issues_repo_name: Thiago-NovaesB/alpheus

max_forks_repo_name: Thiago-NovaesB/alpheus


text:
# Modified from https://juliagpu.gitlab.io/KernelAbstractions.jl/examples/matmul/

using KernelAbstractions, Test, CUDAapi, BenchmarkTools

if CUDAapi.has_cuda_gpu()
    using CuArrays
    CuArrays.allowscalar(false)
end

# Simple kernel for matrix multiplication
@kernel function matmul_kernel!(@Const(a), @Const(b), c)
    i, j = @index(Global, NTuple)

    # creating a temporary sum variable for matrix multiplication
    tmp_sum = zero(eltype(c))
    for k = 1:size(a)[2]
        tmp_sum += a[i,k] * b[k, j]
    end
    c[i,j] = tmp_sum
end

@kernel function matmul_kernel_nothing!(@Const(a), @Const(b), c)
    i, j = @index(Global, NTuple)

    # creating a temporary sum variable for matrix multiplication
    tmp_sum = zero(eltype(c))
    for k = 1:size(a)[2]
        tmp_sum += a[i,k] * b[k, j]
    end
    c[i,j] = tmp_sum
    return nothing
end

# Simple kernel for matrix multiplication
@kernel function matmul_kernel_2!(a, b, c)
   i, j = @index(Global, NTuple)

    # creating a temporary sum variable for matrix multiplication
    c[i,j] = zero(eltype(c))
    for k = 1:size(a)[2]
        c[i,j] += a[i,k] * b[k, j]
    end

end

@kernel function matmul_kernel_3!(@Const(a), @Const(b), c)
    i, j = @index(Global, NTuple)
    # creating a temporary sum variable for matrix multiplication
    c[i,j] = zero(eltype(c))
    for k = 1:size(a)[2]
        c[i,j] += a[i,k] * b[k, j]
    end
end

# Creating a wrapper kernel for launching with error checks
function matmul!(a, b, c)
    if size(a)[2] != size(b)[1]
        println("Matrix size mismatch!")
        return nothing
    end
    if isa(a, Array)
        kernel! = matmul_kernel!(CPU(),4)
    else
        kernel! = matmul_kernel!(CUDA(),256)
    end
    kernel!(a, b, c, ndrange=size(c))
end

function closure_matmul!(matmul_kernel!; cpu_threads = 4, gpu_blocks = 256)
    # Creating a wrapper kernel for launching with error checks
    function matmul!(a, b, c)
        if size(a)[2] != size(b)[1]
            println("Matrix size mismatch!")
            return nothing
        end
        if isa(a, Array)
            kernel! = matmul_kernel!(CPU(), cpu_threads)
        else
            kernel! = matmul_kernel!(CUDA(), gpu_blocks)
        end
        kernel!(a, b, c, ndrange=size(c))
    end
    return matmul!
end

# Creating a wrapper kernel for launching with error checks
function matmul_2!(a, b, c)
    if size(a)[2] != size(b)[1]
        println("Matrix size mismatch!")
        return nothing
    end
    if isa(a, Array)
        kernel! = matmul_kernel_2!(CPU(),4)
    else
        kernel! = matmul_kernel_2!(CUDA(),256)
    end
    kernel!(a, b, c, ndrange=size(c))

end

function matmul_3!(a, b, c)
    if size(a)[2] != size(b)[1]
        println("Matrix size mismatch!")
        return nothing
    end
    if isa(a, Array)
        kernel! = matmul_kernel_3!(CPU(),4)
    else
        kernel! = matmul_kernel_3!(CUDA(),256)
    end
    kernel!(a, b, c, ndrange=size(c))
end

function matmul_nothing!(a, b, c)
    if size(a)[2] != size(b)[1]
        println("Matrix size mismatch!")
        return nothing
    end
    if isa(a, Array)
        kernel! = matmul_kernel_nothing!(CPU(), 4)
    else
        kernel! = matmul_kernel_nothing!(CUDA(),256)
    end
    kernel!(a, b, c, ndrange=size(c))
end

a = rand(256,123)
b = rand(123, 45)
c = zeros(256, 45)

# beginning CPU tests, returns event
ev = matmul!(a,b,c)
wait(ev)

@test isapprox(c, a*b)

ev = matmul_nothing!(a,b,c)
wait(ev)

@test isapprox(c, a*b)

#=
@benchmark wait(matmul!($a, $b, $c))
@benchmark wait(matmul_2!($a, $b, $c))
@benchmark wait(matmul_3!($a, $b, $c))
=#

# beginning GPU tests
if has_cuda_gpu()
    d_a = CuArray(a)
    d_b = CuArray(b)
    d_c = CuArray(c)

    ev = matmul!(d_a, d_b, d_c)
    wait(ev)

    @test isapprox(Array(d_c), d_a*d_b)
    @benchmark wait(matmul!($d_a, $d_b, $d_c))
    @benchmark wait(matmul_2!($d_a, $d_b, $d_c))
    @benchmark wait(matmul_3!($d_a, $d_b, $d_c))
end

max_stars_repo_name: sandreza/Threadicare


text:
using JSON
using MD5
using Sockets

function formatError(p, incNr)
    write(p.socket, "You deliviered the wrong format. This incident will be reported. $incNr\n")
end

function updateNote(p::Player, codetype)
    if codetype == "bombcode"
        notename = "note"
    elseif codetype == "cheesecode"
        notename = "leaflet"
    end

    if !(isfile(string("data/.", notename)))
        #println("No notes files found")
        notes = []
    else
        #println("Notes file found")
        f = open_file_try(string("data/.", notename), "r")
        notes = JSON.parse(read(f, String))
        close(f)
    end

    current_length = length(notes)

    new_notes = [p.dimension]
    if (0 <= current_length < note_max_length)
        append!(new_notes, notes)
    elseif (current_length == note_max_length)
        append!(new_notes, notes[1:note_max_length-1])

        old_path = string("data/.", codetype, "_", notes[note_max_length])
        rm(old_path)
    else
        return
    end

    f = open_file_try(string("data/.", notename), "w")
    write(f, JSON.json(new_notes))
    close(f)

end

function updateOS(p::Player)
    print_dict(p, "cryptomat_os_update_mode")
    user_input = readline(p.socket)
    while true
        if (user_input == "💣" || user_input == "🧀")
            break
        elseif (user_input == "l")
            return
        else
            print_dict("crpytomat_os_update_invalid_mode")
            return
        end
    end
    mode = user_input

    print_dict(p, "cryptomat_os_update_1")
    new_os = readline(p.socket)

    #validate input
    try
        new_os = JSON.parse(new_os)
    catch y
        formatError(p, 0)
        return
    end

    if !(isa(new_os, Array{Any, 1})
        && length(new_os) == 2
        && isa(new_os[1], String)
        && isa(new_os[2], Array{Any, 1}))
        formatError(p, 1)
        return
    end

    for a in new_os[2]
        if !isa(a, Int)
            formatError(p, 2)
            return
        end
    end

    print_dict(p, "cryptomat_os_update_accept_format")

    new_os_json = JSON.json(new_os)

    try
        proc = run(`./cryptomat/rsa_sig.py $new_os_json`)
    catch
        write(p.socket, "Not accepting this signature\n")
        return
    end

    print_dict(p, "cryptomat_os_update_accept_signature")

    os_string = new_os[1]

    write(p.socket, "Updating...\n")
    #println(os_string)
    if mode == "💣"
        codetype = "bombcode"
    elseif mode == "🧀"
        codetype = "cheesecode"
    end
    path = string("data/.", codetype ,"_", p.dimension)
    f = open_file_try(path, "w")
    write(f, os_string)
    close(f)

    updateNote(p, codetype)
    write(p.socket, "Updated\n")

end


text:
using MOCNeutronTransport
@testset "QuadraticSegment_3D" begin
    for T in [Float32, Float64]
        @testset "Constructors" begin
            # Constructor
            x⃗₁ = Point_3D(T, 0, 0, 0)
            x⃗₂ = Point_3D(T, 2, 0, 0)
            x⃗₃ = Point_3D(T, 1, 1, 0)

            q = QuadraticSegment_3D(x⃗₁, x⃗₂, x⃗₃)
            @test q.points == (x⃗₁, x⃗₂, x⃗₃)
        end

        @testset "Methods" begin
            # interpolation
            x⃗₁ = Point_3D(T, 0, 0, 0)
            x⃗₂ = Point_3D(T, 2, 0, 0)
            x⃗₃ = Point_3D(T, 1, 1, 0)

            q = QuadraticSegment_3D(x⃗₁, x⃗₂, x⃗₃)
            @test q.points == (x⃗₁, x⃗₂, x⃗₃)
            for t = LinRange{T}(0, 1, 11)
                @test q(t) ≈ Point_3D(T, 2t, -(2t)^2 + 4t)
            end

#            # arc_length
#            x⃗₁ = Point_3D(T, 0, 0, 0)
#            x⃗₂ = Point_3D(T, 2, 0, 0)
#            x⃗₃ = Point_3D(T, 1, 0, 0)
#            q = QuadraticSegment_3D(x⃗₁, x⃗₂, x⃗₃)
#
#            # straight edge
#            @test abs(arc_length(q) - T(2)) < 1.0e-6
#            # curved
#            x⃗₃ = Point_3D(T, 1, 1, 0)
#            q = QuadraticSegment_3D(x⃗₁, x⃗₂, x⃗₃)
#            @test abs(arc_length(q) - T(2.9578857151786138)) < 1.0e-6
#
#            # intersect
#            x⃗₁ = Point_3D(T, 0, 0, 0)
#            x⃗₂ = Point_3D(T, 2, 0, 0)
#            x⃗₃ = Point_3D(T, 1, 1, 0)
#            x⃗₄ = Point_3D(T, 1, 0, 0)
#            x⃗₅ = Point_3D(T, 1, 2, 0)
#
#            # 1 intersection
#            q = QuadraticSegment_3D(x⃗₁, x⃗₂, x⃗₃)
#            l = LineSegment_3D(x⃗₄, x⃗₅)
#            bool, npoints, points = intersect(l, q)
#            @test bool
#            @test npoints == 1
#            @test points[1] ≈ Point_3D(T, 1, 1, 0)
#
#            # 2 intersections
#            x⃗₄ = Point_3D(T, 0, 3//4, 0)
#            x⃗₅ = Point_3D(T, 2, 3//4, 0)
#            l = LineSegment_3D(x⃗₄, x⃗₅)
#            bool, npoints, points = l ∩ q
#            @test bool
#            @test npoints == 2
#            @test points[1] ≈ Point_3D(T, 1//2, 3//4, 0)
#            @test points[2] ≈ Point_3D(T, 3//2, 3//4, 0)
#
#            # 0 intersections
#            x⃗₄ = Point_3D(T, 0, 3, 0)
#            x⃗₅ = Point_3D(T, 2, 3, 0)
#            l = LineSegment_3D(x⃗₄, x⃗₅)
#            bool, npoints, points = intersect(l, q)
#            @test !bool
#            @test npoints == 0
#            @test points[1] ≈ Point_3D(T, 0)
#            @test points[2] ≈ Point_3D(T, 0)
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
using FFTW

function compressaudiofft(channels, threshold)
    # Split channels
    left = channels[:,1]
    right = channels[:,2]

    # Perfom transforms
    fleft = fft(left)
    fright = fft(right)

    # Compress
    tleft = threshold * maximum(abs.(fleft))
    ileft = abs.(fleft) .> tleft
    tright = threshold * maximum(abs.(fright))
    iright = abs.(fright) .> tright
    return (fleft[ileft], ileft, fright[iright], iright, length(left))
end

function decompressaudiofft(compression)
    # Unpack compression
    cleft, ileft, cright, iright, len = compression

    # Decompress left channel
    fleft = Array{Complex{Float64},1}(undef, length(ileft))
    j = 1
    for i = 1:length(fleft)
        if ileft[i]
            fleft[i] = cleft[j]
            j += 1
        end
    end

    # Decompress right channel
    fright = Array{Complex{Float64},1}(undef, length(iright))
    j = 1
    for i = 1:length(fright)
        if iright[i]
            fright[i] = cright[j]
            j += 1
        end
    end

    # Perform inverse transforms
    left = ifft(fleft)
    right = ifft(fright)

    return real.([left[1:len] right[1:len]])
end


text:
calcMidpoints(edges::AbstractVector) = Float64[0.5 * (edges[i] + edges[i+1]) for i in 1:length(edges)-1]

"Make histogram-like bins of data"
function binData(data, nbins)
  lo, hi = extrema(data)
  edges = collect(linspace(lo, hi, nbins+1))
  midpoints = calcMidpoints(edges)
  buckets = Int[max(2, min(searchsortedfirst(edges, x), length(edges)))-1 for x in data]
  counts = zeros(Int, length(midpoints))
  for b in buckets
    counts[b] += 1
  end
  edges, midpoints, buckets, counts
end

"""
A hacky replacement for a histogram when the backend doesn't support histograms directly.
Convert it into a bar chart with the appropriate x/y values.
"""
function histogramHack(; kw...)
  d = KW(kw)

  # we assume that the y kwarg is set with the data to be binned, and nbins is also defined
  edges, midpoints, buckets, counts = binData(d[:y], d[:bins])
  d[:x] = midpoints
  d[:y] = float(counts)
  d[:seriestype] = :bar
  d[:fillrange] = d[:fillrange] == nothing ? 0.0 : d[:fillrange]
  d
end

"""
A hacky replacement for a bar graph when the backend doesn't support bars directly.
Convert it into a line chart with fillrange set.
"""
function barHack(; kw...)
  d = KW(kw)
  midpoints = d[:x]
  heights = d[:y]
  fillrange = d[:fillrange] == nothing ? 0.0 : d[:fillrange]

  # estimate the edges
  dists = diff(midpoints) * 0.5
  edges = zeros(length(midpoints)+1)
  for i in 1:length(edges)
    if i == 1
      edge = midpoints[1] - dists[1]
    elseif i == length(edges)
      edge = midpoints[i-1] + dists[i-2]
    else
      edge = midpoints[i-1] + dists[i-1]
    end
    edges[i] = edge
  end

  x = Float64[]
  y = Float64[]
  for i in 1:length(heights)
    e1, e2 = edges[i:i+1]
    append!(x, [e1, e1, e2, e2])
    append!(y, [fillrange, heights[i], heights[i], fillrange])
  end

  d[:x] = x
  d[:y] = y
  d[:seriestype] = :path
  d[:fillrange] = fillrange
  d
end

"""
A hacky replacement for a sticks graph when the backend doesn't support sticks directly.
Convert it into a line chart that traces the sticks, and a scatter that sets markers at the points.
"""
function sticksHack(; kw...)
  dLine = KW(kw)
  dScatter = copy(dLine)

  # these are the line vertices
  x = Float64[]
  y = Float64[]
  fillrange = dLine[:fillrange] == nothing ? 0.0 : dLine[:fillrange]

  # calculate the vertices
  yScatter = dScatter[:y]
  for (i,xi) in enumerate(dScatter[:x])
    yi = yScatter[i]
    for j in 1:3 push!(x, xi) end
    append!(y, [fillrange, yScatter[i], fillrange])
  end

  # change the line args
  dLine[:x] = x
  dLine[:y] = y
  dLine[:seriestype] = :path
  dLine[:markershape] = :none
  dLine[:fillrange] = nothing

  # change the scatter args
  dScatter[:seriestype] = :none

  dLine, dScatter
end

function regressionXY(x, y)
  # regress
  β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)

  # make a line segment
  regx = [minimum(x), maximum(x)]
  regy = β * regx + α
  regx, regy
end

function replace_image_with_heatmap{T<:Colorant}(z::Array{T})
    @show T, size(z)
    n, m = size(z)
    # idx = 0
    colors = ColorGradient(vec(z))
    newz = reshape(linspace(0, 1, n*m), n, m)
    newz, colors
    # newz = zeros(n, m)
    # for i=1:n, j=1:m
    #     push!(colors, T(z[i,j]...))
    #     newz[i,j] = idx / (n*m-1)
    #     idx += 1
    # end
    # newz, ColorGradient(colors)
end

function imageHack(d::KW)
    is_seriestype_supported(:heatmap) || error("Neither :image or :heatmap are supported!")
    d[:seriestype] = :heatmap
    d[:z], d[:fillcolor] = replace_image_with_heatmap(d[:z].surf)
end
# ---------------------------------------------------------------

type Segments{T}
    pts::Vector{T}
end

# Segments() = Segments{Float64}(zeros(0))

Segments() = Segments(Float64)
Segments{T}(::Type{T}) = Segments(T[])
Segments(p::Int) = Segments(NTuple{2,Float64}[])

# Segments() = Segments(zeros(0))

to_nan(::Type{Float64}) = NaN
to_nan(::Type{NTuple{2,Float64}}) = (NaN, NaN)

coords(segs::Segments{Float64}) = segs.pts
coords(segs::Segments{NTuple{2,Float64}}) = Float64[p[1] for p in segs.pts], Float64[p[2] for p in segs.pts]

function Base.push!{T}(segments::Segments{T}, vs...)
    if !isempty(segments.pts)
        push!(segments.pts, to_nan(T))
    end
    for v in vs
        push!(segments.pts, convert(T,v))
    end
    segments
end

function Base.push!{T}(segments::Segments{T}, vs::AVec)
    if !isempty(segments.pts)
        push!(segments.pts, to_nan(T))
    end
    for v in vs
        push!(segments.pts, convert(T,v))
    end
    segments
end

# -----------------------------------------------------
# helper to manage NaN-separated segments

type SegmentsIterator
    args::Tuple
    n::Int
end
function iter_segments(args...)
    tup = Plots.wraptuple(args)
    n = maximum(map(length, tup))
    SegmentsIterator(tup, n)
end

# helpers to figure out if there are NaN values in a list of array types
anynan(i::Int, args::Tuple) = any(a -> !isfinite(cycle(a,i)), args)
anynan(istart::Int, iend::Int, args::Tuple) = any(i -> anynan(i, args), istart:iend)
allnan(istart::Int, iend::Int, args::Tuple) = all(i -> anynan(i, args), istart:iend)

function Base.start(itr::SegmentsIterator)
    nextidx = 1
    if anynan(1, itr.args)
        _, nextidx = next(itr, 1)
    end
    nextidx
end
Base.done(itr::SegmentsIterator, nextidx::Int) = nextidx > itr.n
function Base.next(itr::SegmentsIterator, nextidx::Int)
    i = istart = iend = nextidx

    # find the next NaN, and iend is the one before
    while i <= itr.n + 1
        if i > itr.n || anynan(i, itr.args)
            # done... array end or found NaN
            iend = i-1
            break
        end
        i += 1
    end

    # find the next non-NaN, and set nextidx
    while i <= itr.n
        if !anynan(i, itr.args)
            break
        end
        i += 1
    end

    istart:iend, i
end

# ------------------------------------------------------------------------------------

nop() = nothing
notimpl() = error("This has not been implemented yet")

Base.cycle(wrapper::InputWrapper, idx::Int) = wrapper.obj
Base.cycle(wrapper::InputWrapper, idx::AVec{Int}) = wrapper.obj

Base.cycle(v::AVec, idx::Int) = v[mod1(idx, length(v))]
Base.cycle(v::AMat, idx::Int) = size(v,1) == 1 ? v[1, mod1(idx, size(v,2))] : v[:, mod1(idx, size(v,2))]
Base.cycle(v, idx::Int)       = v

Base.cycle(v::AVec, indices::AVec{Int}) = map(i -> cycle(v,i), indices)
Base.cycle(v::AMat, indices::AVec{Int}) = map(i -> cycle(v,i), indices)
Base.cycle(v, indices::AVec{Int})       = fill(v, length(indices))

Base.cycle(grad::ColorGradient, idx::Int) = cycle(grad.colors, idx)
Base.cycle(grad::ColorGradient, indices::AVec{Int}) = cycle(grad.colors, indices)

makevec(v::AVec) = v
makevec{T}(v::T) = T[v]

"duplicate a single value, or pass the 2-tuple through"
maketuple(x::Real)                     = (x,x)
maketuple{T,S}(x::Tuple{T,S}) = x

mapFuncOrFuncs(f::Function, u::AVec)        = map(f, u)
mapFuncOrFuncs(fs::AVec{Function}, u::AVec) = [map(f, u) for f in fs]

unzip{X,Y}(xy::AVec{Tuple{X,Y}})              = [t[1] for t in xy], [t[2] for t in xy]
unzip{X,Y,Z}(xyz::AVec{Tuple{X,Y,Z}})         = [t[1] for t in xyz], [t[2] for t in xyz], [t[3] for t in xyz]
unzip{X,Y,U,V}(xyuv::AVec{Tuple{X,Y,U,V}})    = [t[1] for t in xyuv], [t[2] for t in xyuv], [t[3] for t in xyuv], [t[4] for t in xyuv]

unzip{T}(xy::AVec{FixedSizeArrays.Vec{2,T}})  = T[t[1] for t in xy], T[t[2] for t in xy]
unzip{T}(xy::FixedSizeArrays.Vec{2,T})        = T[xy[1]], T[xy[2]]

unzip{T}(xyz::AVec{FixedSizeArrays.Vec{3,T}}) = T[t[1] for t in xyz], T[t[2] for t in xyz], T[t[3] for t in xyz]
unzip{T}(xyz::FixedSizeArrays.Vec{3,T})       = T[xyz[1]], T[xyz[2]], T[xyz[3]]

unzip{T}(xyuv::AVec{FixedSizeArrays.Vec{4,T}}) = T[t[1] for t in xyuv], T[t[2] for t in xyuv], T[t[3] for t in xyuv], T[t[4] for t in xyuv]
unzip{T}(xyuv::FixedSizeArrays.Vec{4,T})       = T[xyuv[1]], T[xyuv[2]], T[xyuv[3]], T[xyuv[4]]

# given 2-element lims and a vector of data x, widen lims to account for the extrema of x
function _expand_limits(lims, x)
  try
    e1, e2 = extrema(x)
    lims[1] = min(lims[1], e1)
    lims[2] = max(lims[2], e2)
  # catch err
  #   warn(err)
  end
  nothing
end

expand_data(v, n::Integer) = [cycle(v, i) for i=1:n]

# if the type exists in a list, replace the first occurence.  otherwise add it to the end
function addOrReplace(v::AbstractVector, t::DataType, args...; kw...)
    for (i,vi) in enumerate(v)
        if isa(vi, t)
            v[i] = t(args...; kw...)
            return
        end
    end
    push!(v, t(args...; kw...))
    return
end

function replaceType(vec, val)
  filter!(x -> !isa(x, typeof(val)), vec)
  push!(vec, val)
end

function replaceAlias!(d::KW, k::Symbol, aliases::Dict{Symbol,Symbol})
  if haskey(aliases, k)
    d[aliases[k]] = pop!(d, k)
  end
end

function replaceAliases!(d::KW, aliases::Dict{Symbol,Symbol})
  ks = collect(keys(d))
  for k in ks
      replaceAlias!(d, k, aliases)
  end
end

createSegments(z) = collect(repmat(z',2,1))[2:end]

Base.first(c::Colorant) = c
Base.first(x::Symbol) = x

sortedkeys(d::Dict) = sort(collect(keys(d)))

"create an (n+1) list of the outsides of heatmap rectangles"
function heatmap_edges(v::AVec)
  vmin, vmax = extrema(v)
  extra = 0.5 * (vmax-vmin) / (length(v)-1)
  vcat(vmin-extra, 0.5 * (v[1:end-1] + v[2:end]), vmax+extra)
end

function calc_r_extrema(x, y)
    xmin, xmax = extrema(x)
    ymin, ymax = extrema(y)
    r = 0.5 * min(xmax - xmin, ymax - ymin)
    extrema(r)
end

function convert_to_polar(x, y, r_extrema = calc_r_extrema(x, y))
    rmin, rmax = r_extrema
    phi, r = x, y
    r = 0.5 * (r - rmin) / (rmax - rmin)
    n = max(length(phi), length(r))
    x = zeros(n)
    y = zeros(n)
    for i in 1:n
        x[i] = cycle(r,i) * cos(cycle(phi,i))
        y[i] = cycle(r,i) * sin(cycle(phi,i))
    end
    x, y
end

function fakedata(sz...)
  y = zeros(sz...)
  for r in 2:size(y,1)
    y[r,:] = 0.95 * vec(y[r-1,:]) + randn(size(y,2))
  end
  y
end

isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited
isatom() = isdefined(Main, :Atom) && Main.Atom.isconnected()

function is_installed(pkgstr::AbstractString)
    try
        Pkg.installed(pkgstr) === nothing ? false: true
    catch
        false
    end
end

istuple(::Tuple) = true
istuple(::Any)   = false
isvector(::AVec) = true
isvector(::Any)  = false
ismatrix(::AMat) = true
ismatrix(::Any)  = false
isscalar(::Real) = true
isscalar(::Any)  = false

is_2tuple(v) = typeof(v) <: Tuple && length(v) == 2

isvertical(d::KW) = get(d, :orientation, :vertical) in (:vertical, :v, :vert)
isvertical(series::Series) = isvertical(series.d)

ticksType{T<:Real}(ticks::AVec{T})                      = :ticks
ticksType{T<:AbstractString}(ticks::AVec{T})            = :labels
ticksType{T<:AVec,S<:AVec}(ticks::Tuple{T,S})  = :ticks_and_labels
ticksType(ticks)                                        = :invalid

limsType{T<:Real,S<:Real}(lims::Tuple{T,S})    = :limits
limsType(lims::Symbol)                                  = lims == :auto ? :auto : :invalid
limsType(lims)                                          = :invalid

# axis_Symbol(letter, postfix) = Symbol(letter * postfix)
# axis_symbols(letter, postfix...) = map(s -> axis_Symbol(letter, s), postfix)

Base.convert{T<:Real}(::Type{Vector{T}}, rng::Range{T})         = T[x for x in rng]
Base.convert{T<:Real,S<:Real}(::Type{Vector{T}}, rng::Range{S}) = T[x for x in rng]

Base.merge(a::AbstractVector, b::AbstractVector) = sort(unique(vcat(a,b)))

nanpush!(a::AbstractVector, b) = (push!(a, NaN); push!(a, b))
nanappend!(a::AbstractVector, b) = (push!(a, NaN); append!(a, b))

function nansplit(v::AVec)
    vs = Vector{eltype(v)}[]
    while true
        idx = findfirst(isnan, v)
        if idx <= 0
            # no nans
            push!(vs, v)
            break
        elseif idx > 1
            push!(vs, v[1:idx-1])
        end
        v = v[idx+1:end]
    end
    vs
end

function nanvcat(vs::AVec)
    v_out = zeros(0)
    for v in vs
        nanappend!(v_out, v)
    end
    v_out
end

# given an array of discrete values, turn it into an array of indices of the unique values
# returns the array of indices (znew) and a vector of unique values (vals)
function indices_and_unique_values(z::AbstractArray)
    vals = sort(unique(z))
    vmap = Dict([(v,i) for (i,v) in enumerate(vals)])
    newz = map(zi -> vmap[zi], z)
    newz, vals
end

# this is a helper function to determine whether we need to transpose a surface matrix.
# it depends on whether the backend matches rows to x (transpose_on_match == true) or vice versa
# for example: PyPlot sends rows to y, so transpose_on_match should be true
function transpose_z(d, z, transpose_on_match::Bool = true)
    if d[:match_dimensions] == transpose_on_match
        # z'
        permutedims(z, [2,1])
    else
        z
    end
end

function ok(x::Number, y::Number, z::Number = 0)
    isfinite(x) && isfinite(y) && isfinite(z)
end
ok(tup::Tuple) = ok(tup...)

# compute one side of a fill range from a ribbon
function make_fillrange_side(y, rib)
    frs = zeros(length(y))
    for (i, (yi, ri)) in enumerate(zip(y, cycle(rib)))
        frs[i] = yi + ri
    end
    frs
end

# turn a ribbon into a fillrange
function make_fillrange_from_ribbon(kw::KW)
    y, rib = kw[:y], kw[:ribbon]
    rib = wraptuple(rib)
    rib1, rib2 = -first(rib), last(rib)
    kw[:ribbon] = nothing
    kw[:fillrange] = make_fillrange_side(y, rib1), make_fillrange_side(y, rib2)
end

function get_sp_lims(sp::Subplot, letter::Symbol)
    axis_limits(sp[Symbol(letter, :axis)])
end
xlims(sp::Subplot) = get_sp_lims(sp, :x)
ylims(sp::Subplot) = get_sp_lims(sp, :y)
zlims(sp::Subplot) = get_sp_lims(sp, :z)
xlims(plt::Plot, sp_idx::Int = 1) = xlims(plt[sp_idx])
ylims(plt::Plot, sp_idx::Int = 1) = ylims(plt[sp_idx])
zlims(plt::Plot, sp_idx::Int = 1) = zlims(plt[sp_idx])
xlims(sp_idx::Int = 1) = xlims(current(), sp_idx)
ylims(sp_idx::Int = 1) = ylims(current(), sp_idx)
zlims(sp_idx::Int = 1) = zlims(current(), sp_idx)

# ---------------------------------------------------------------

wraptuple(x::Tuple) = x
wraptuple(x) = (x,)

trueOrAllTrue(f::Function, x::AbstractArray) = all(f, x)
trueOrAllTrue(f::Function, x) = f(x)

allLineTypes(arg)   = trueOrAllTrue(a -> get(_typeAliases, a, a) in _allTypes, arg)
allStyles(arg)      = trueOrAllTrue(a -> get(_styleAliases, a, a) in _allStyles, arg)
allShapes(arg)      = trueOrAllTrue(a -> is_marker_supported(get(_markerAliases, a, a)), arg) ||
                        trueOrAllTrue(a -> isa(a, Shape), arg)
allAlphas(arg)      = trueOrAllTrue(a -> (typeof(a) <: Real && a > 0 && a < 1) ||
                        (typeof(a) <: AbstractFloat && (a == zero(typeof(a)) || a == one(typeof(a)))), arg)
allReals(arg)       = trueOrAllTrue(a -> typeof(a) <: Real, arg)
allFunctions(arg)   = trueOrAllTrue(a -> isa(a, Function), arg)

# ---------------------------------------------------------------
# ---------------------------------------------------------------

"""
Allows temporary setting of backend and defaults for Plots. Settings apply only for the `do` block.  Example:
```
with(:gadfly, size=(400,400), type=:histogram) do
  plot(rand(10))
  plot(rand(10))
end
```
"""
function with(f::Function, args...; kw...)
    newdefs = KW(kw)

    if :canvas in args
        newdefs[:xticks] = nothing
        newdefs[:yticks] = nothing
        newdefs[:grid] = false
        newdefs[:legend] = false
    end

  # dict to store old and new keyword args for anything that changes
  olddefs = KW()
  for k in keys(newdefs)
    olddefs[k] = default(k)
  end

  # save the backend
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end
  oldbackend = CURRENT_BACKEND.sym

  for arg in args

    # change backend?
    if arg in backends()
      backend(arg)
    end

    # # TODO: generalize this strategy to allow args as much as possible
    # #       as in:  with(:gadfly, :scatter, :legend, :grid) do; ...; end
    # # TODO: can we generalize this enough to also do something similar in the plot commands??

    # k = :seriestype
    # if arg in _allTypes
    #   olddefs[k] = default(k)
    #   newdefs[k] = arg
    # elseif haskey(_typeAliases, arg)
    #   olddefs[k] = default(k)
    #   newdefs[k] = _typeAliases[arg]
    # end

    k = :legend
    if arg in (k, :leg)
      olddefs[k] = default(k)
      newdefs[k] = true
    end

    k = :grid
    if arg == k
      olddefs[k] = default(k)
      newdefs[k] = true
    end
  end

  # display(olddefs)
  # display(newdefs)

  # now set all those defaults
  default(; newdefs...)

  # call the function
  ret = f()

  # put the defaults back
  default(; olddefs...)

  # revert the backend
  if CURRENT_BACKEND.sym != oldbackend
    backend(oldbackend)
  end

  # return the result of the function
  ret
end

# ---------------------------------------------------------------
# ---------------------------------------------------------------

type DebugMode
  on::Bool
end
const _debugMode = DebugMode(false)

function debugplots(on = true)
  _debugMode.on = on
end

debugshow(x) = show(x)
debugshow(x::AbstractArray) = print(summary(x))

function dumpdict(d::KW, prefix = "", alwaysshow = false)
  _debugMode.on || alwaysshow || return
  println()
  if prefix != ""
    println(prefix, ":")
  end
  for k in sort(collect(keys(d)))
    @printf("%14s: ", k)
    debugshow(d[k])
    println()
  end
  println()
end
DD(d::KW, prefix = "") = dumpdict(d, prefix, true)

function dumpcallstack()
  error()  # well... you wanted the stacktrace, didn't you?!?
end

# ---------------------------------------------------------------
# ---------------------------------------------------------------
# used in updating an existing series

extendSeriesByOne(v::UnitRange{Int}, n::Int = 1) = isempty(v) ? (1:n) : (minimum(v):maximum(v)+n)
extendSeriesByOne(v::AVec, n::Integer = 1)       = isempty(v) ? (1:n) : vcat(v, (1:n) + maximum(v))
extendSeriesData{T}(v::Range{T}, z::Real)        = extendSeriesData(float(collect(v)), z)
extendSeriesData{T}(v::Range{T}, z::AVec)        = extendSeriesData(float(collect(v)), z)
extendSeriesData{T}(v::AVec{T}, z::Real)         = (push!(v, convert(T, z)); v)
extendSeriesData{T}(v::AVec{T}, z::AVec)         = (append!(v, convert(Vector{T}, z)); v)

# -------------------------------------------------------
# NOTE: backends should implement the following methods to get/set the x/y/z data objects

tovec(v::AbstractVector) = v
tovec(v::Void) = zeros(0)

function getxy(plt::Plot, i::Integer)
    d = plt.series_list[i].d
    tovec(d[:x]), tovec(d[:y])
end
function getxyz(plt::Plot, i::Integer)
    d = plt.series_list[i].d
    tovec(d[:x]), tovec(d[:y]), tovec(d[:z])
end

function setxy!{X,Y}(plt::Plot, xy::Tuple{X,Y}, i::Integer)
    series = plt.series_list[i]
    series.d[:x], series.d[:y] = xy
    sp = series.d[:subplot]
    reset_extrema!(sp)
    _series_updated(plt, series)
end
function setxyz!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer)
    series = plt.series_list[i]
    series.d[:x], series.d[:y], series.d[:z] = xyz
    sp = series.d[:subplot]
    reset_extrema!(sp)
    _series_updated(plt, series)
end

function setxyz!{X,Y,Z<:AbstractMatrix}(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer)
    setxyz!(plt, (xyz[1], xyz[2], Surface(xyz[3])), i)
end

# -------------------------------------------------------
# indexing notation

# Base.getindex(plt::Plot, i::Integer) = getxy(plt, i)
Base.setindex!{X,Y}(plt::Plot, xy::Tuple{X,Y}, i::Integer) = (setxy!(plt, xy, i); plt)
Base.setindex!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z}, i::Integer) = (setxyz!(plt, xyz, i); plt)

# -------------------------------------------------------

# operate on individual series

function push_x!(series::Series, xi)
    push!(series[:x], xi)
    expand_extrema!(series[:subplot][:xaxis], xi)
    return
end
function push_y!(series::Series, yi)
    push!(series[:y], yi)
    expand_extrema!(series[:subplot][:yaxis], yi)
    return
end
function push_z!(series::Series, zi)
    push!(series[:z], zi)
    expand_extrema!(series[:subplot][:zaxis], zi)
    return
end

function Base.push!(series::Series, yi)
    x = extendSeriesByOne(series[:x])
    expand_extrema!(series[:subplot][:xaxis], x[end])
    series[:x] = x
    push_y!(series, yi)
end
Base.push!(series::Series, xi, yi) = (push_x!(series,xi); push_y!(series,yi))
Base.push!(series::Series, xi, yi, zi) = (push_x!(series,xi); push_y!(series,yi); push_z!(series,zi))

# -------------------------------------------------------

function update!(series::Series; kw...)
    d = KW(kw)
    preprocessArgs!(d)
    for (k,v) in d
        if haskey(_series_defaults, k)
            series[k] = v
        else
            warn("unused key $k in series update")
        end
    end
    _series_updated(series[:subplot].plt, series)
    series
end

function update!(sp::Subplot; kw...)
    d = KW(kw)
    preprocessArgs!(d)
    for (k,v) in d
        if haskey(_subplot_defaults, k)
            sp[k] = v
        else
            warn("unused key $k in subplot update")
        end
    end
    sp
end

# -------------------------------------------------------
# push/append for one series

# push value to first series
Base.push!(plt::Plot, y::Real) = push!(plt, 1, y)
Base.push!(plt::Plot, x::Real, y::Real) = push!(plt, 1, x, y)
Base.push!(plt::Plot, x::Real, y::Real, z::Real) = push!(plt, 1, x, y, z)

# y only
function Base.push!(plt::Plot, i::Integer, y::Real)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesByOne(xdata), extendSeriesData(ydata, y)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, y::AVec)
    xdata, ydata = plt[i]
    if !isa(xdata, UnitRange{Int})
        error("Expected x is a UnitRange since you're trying to push a y value only")
    end
    plt[i] = (extendSeriesByOne(xdata, length(y)), extendSeriesData(ydata, y))
    plt
end

# x and y
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec)
    @assert length(x) == length(y)
    xdata, ydata = getxy(plt, i)
    setxy!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y)), i)
    plt
end

# x, y, and z
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real, z::Real)
    # @show i, x, y, z
    xdata, ydata, zdata = getxyz(plt, i)
    # @show xdata, ydata, zdata
    setxyz!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y), extendSeriesData(zdata, z)), i)
    plt
end
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec, z::AVec)
    @assert length(x) == length(y) == length(z)
    xdata, ydata, zdata = getxyz(plt, i)
    setxyz!(plt, (extendSeriesData(xdata, x), extendSeriesData(ydata, y), extendSeriesData(zdata, z)), i)
    plt
end

# tuples
Base.push!{X,Y}(plt::Plot, xy::Tuple{X,Y})                  = push!(plt, 1, xy...)
Base.push!{X,Y,Z}(plt::Plot, xyz::Tuple{X,Y,Z})             = push!(plt, 1, xyz...)
Base.push!{X,Y}(plt::Plot, i::Integer, xy::Tuple{X,Y})      = push!(plt, i, xy...)
Base.push!{X,Y,Z}(plt::Plot, i::Integer, xyz::Tuple{X,Y,Z}) = push!(plt, i, xyz...)

# -------------------------------------------------------
# push/append for all series

# push y[i] to the ith series
function Base.push!(plt::Plot, y::AVec)
    ny = length(y)
    for i in 1:plt.n
        push!(plt, i, y[mod1(i,ny)])
    end
    plt
end

# push y[i] to the ith series
# same x for each series
function Base.push!(plt::Plot, x::Real, y::AVec)
    push!(plt, [x], y)
end

# push (x[i], y[i]) to the ith series
function Base.push!(plt::Plot, x::AVec, y::AVec)
    nx = length(x)
    ny = length(y)
    for i in 1:plt.n
        push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)])
    end
    plt
end

# push (x[i], y[i], z[i]) to the ith series
function Base.push!(plt::Plot, x::AVec, y::AVec, z::AVec)
    nx = length(x)
    ny = length(y)
    nz = length(z)
    for i in 1:plt.n
        push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)], z[mod1(i,nz)])
    end
    plt
end

# ---------------------------------------------------------------

# Some conversion functions
# note: I borrowed these conversion constants from Compose.jl's Measure

const PX_PER_INCH   = 100
const DPI           = PX_PER_INCH
const MM_PER_INCH   = 25.4
const MM_PER_PX     = MM_PER_INCH / PX_PER_INCH

inch2px(inches::Real)   = float(inches * PX_PER_INCH)
px2inch(px::Real)       = float(px / PX_PER_INCH)
inch2mm(inches::Real)   = float(inches * MM_PER_INCH)
mm2inch(mm::Real)       = float(mm / MM_PER_INCH)
px2mm(px::Real)         = float(px * MM_PER_PX)
mm2px(mm::Real)         = float(px / MM_PER_PX)

"Smallest x in plot"
xmin(plt::Plot) = minimum([minimum(series.d[:x]) for series in plt.series_list])
"Largest x in plot"
xmax(plt::Plot) = maximum([maximum(series.d[:x]) for series in plt.series_list])

"Extrema of x-values in plot"
Base.extrema(plt::Plot) = (xmin(plt), xmax(plt))


text:
"""
    slice(results::MEstimation_results,
          what::Int64;
          grid::Vector{Float64} = Vector{Float64}(),
          at::Vector{Float64} = Vector{Float64}(),
          n_points::Int64 = 50,
          n_sd::Real = 2)

Compute 1-dimensional slices of objective functions and estimating function surfaces for the parameter `what` over a grid of points `grid`, from a [`MEstimation_results`](@ref) object.

Arguments
===
+ `results`: An [`MEstimation_results`](@ref) object.
+ `what`: the index of the parameter for which to compute a slice for

[Keyword arguments](https://docs.julialang.org/en/v1/manual/functions/#Keyword-Arguments-1)
===
+ `grid`: a `Vector{Float64}`; if supplied, the slice is computed at each element of `grid`. The default will result in the automatic calculation of the grid; see Details.
+ `at`: a `Vector{Float64}` of the sample length as `coef(results)`, specifying the parameter values at which to compute the slice. The default results in computing the slice at `coef(results)`.
+ `n_points`: an `Int64` specifying the number of grid points to generate. Applicable only if `grid` is not supplied; see Details.
+ `n_sd`: an `Int64` specifying the number of standard errors to be used for the grid generation. Applicable only if `grid` is not supplied; see Details.

Result
===
A [`Dict`](https://docs.julialang.org/en/v1/base/collections/#Dictionaries-1) with keys "grid" and "slice", holding `grid` and the values of the slice at `grid`, respectively.

Details
===
The default value of `grid` will result in the automatic calculation of a grid of `n_points` points, between `coef(results)[what] - n_sd * stderror(results)[what]` and `coef(results)[what] + n_sd * stderror(results)[what]`.
"""
function slice(results::MEstimation_results,
               what::Int64;
               grid::Vector{Float64} = Vector{Float64}(),
               at::Vector{Float64} = Vector{Float64}(),
               n_points::Int64 = 50,
               n_sd::Real = 2.0)
    data = results.data
    theta = length(at) > 0 ? at : coef(results)
    n_par = length(theta)
    model_template = results.template
    br = results.br
    if length(grid) == 0
        sd = stderror(results)[what]
        co = coef(results)[what]
        step = 2 * n_sd * sd / n_points
        grid = range(co - n_sd * sd, co + n_sd * sd, length = n_points)
    end
    if results.has_regularizer
        regularizer = results.regularizer
        if results.has_objective
            slice_function = beta -> objective_function(beta, data, model_template, br) + regularizer(beta, data)
        else
            slice_function = beta -> estimating_function(beta, data, model_template, br)[what] + regularizer(beta, data)[what]
        end
    else
        if results.has_objective
            slice_function = beta -> objective_function(beta, data, model_template, br)
        else
            slice_function = beta -> estimating_function(beta, data, model_template, br)[what]
        end
    end
    slice_m = map(function(b) th = theta[1:n_par]; th[what] = b; slice_function(th) end, grid)
    Dict("grid" => grid, "slice" => slice_m)
end


text:
using Test
using CUSTATEVEC
import CUSTATEVEC: CuStateVec, applyMatrix!, expectation, sample

@testset "CUSTATEVEC" begin
    # build a simple state and compute expectations
    n_q = 2
    @testset for elty in [ComplexF32, ComplexF64]
        H = convert(Matrix{elty}, (1/√2).*[1 1; 1 -1])
        X = convert(Matrix{elty}, [0 1; 1 0])
        Z = convert(Matrix{elty}, [1 0; 0 -1])
        sv = CuStateVec(elty, n_q)
        sv = applyMatrix!(sv, H, false, Int32[0], Int32[])
        exp, res = expectation(sv, Z, Int32[0])
        @test exp ≈ 0.0 atol=1e-6
        exp, res = expectation(sv, X, Int32[0])
        @test exp ≈ 1.0 atol=1e-6
    end
    # build a simple state with controls and compute expectations
    n_q = 2
    @testset for elty in [ComplexF32, ComplexF64]
        H = convert(Matrix{elty}, (1/√2).*[1 1; 1 -1])
        X = convert(Matrix{elty}, [0 1; 1 0])
        Z = convert(Matrix{elty}, [1 0; 0 -1])
        sv = CuStateVec(elty, n_q)
        sv = applyMatrix!(sv, H, false, Int32[0], Int32[])
        sv = applyMatrix!(sv, X, false, Int32[1], Int32[0]) # CNOT
        exp, res = expectation(sv, Z, Int32[0])
        @test exp ≈ 0.0 atol=1e-6
        exp, res = expectation(sv, X, Int32[0])
        @test exp ≈ 0.0 atol=1e-6
    end
    # build a simple state and compute samples
    n_q = 10
    @testset for elty in [ComplexF32, ComplexF64]
        H = convert(Matrix{elty}, (1/√2).*[1 1; 1 -1])
        sv = CuStateVec(elty, n_q)
        for q in 0:n_q-1
            sv = applyMatrix!(sv, H, false, Int32[q], Int32[])
        end
        samples = sample(sv, Int32.(collect(0:n_q-1)), 10)
        # check that this succeeded
        @test length(samples) == 10
    end
end


text:
module ESeriesRounding

export E3, E6, E12, E24, E48, E96, E192
export determine_E

import Base: round, print, show
using NumericIO: formatted
include("UtilityFunctions.jl")

struct ESeries
    name::Symbol
    vals::AbstractVector
end

E3  = ESeries(:E3,   [100, 220, 470])
E6  = ESeries(:E6,   [100, 150, 220, 330, 470, 680])
E12 = ESeries(:E12,  [100, 120, 150, 180, 220, 270, 330, 390, 470, 560, 680, 820])
E24 = ESeries(:E24,  [100, 110, 120, 130, 150, 160, 180, 200, 220, 240, 270, 300, 330, 360, 390, 430, 470, 510, 560, 620, 680, 750, 820, 910])
E48 = ESeries(:E48,  [100, 105, 110, 115, 121, 127, 133, 140, 147, 154, 162, 169, 178, 187, 196, 205, 215, 226, 237, 249, 261, 274, 287, 301, 316, 332, 348, 365, 383, 402, 422, 442, 464, 487, 511, 536, 562, 590, 619, 649, 681, 715, 750, 787, 825, 866, 909, 953])
E96 = ESeries(:E96,  [100, 102, 105, 107, 110, 113, 115, 118, 121, 124, 127, 130, 133, 137, 140, 143, 147, 150, 154, 158, 162, 165, 169, 174, 178, 182, 187, 191, 196, 200, 205, 210, 216, 221, 226, 232, 237, 243, 249, 255, 261, 267, 274, 280, 287, 294, 301, 309, 316, 324, 332, 340, 348, 357, 365, 374, 383, 392, 402, 412, 422, 432, 442, 453, 464, 475, 487, 499, 511, 523, 536, 549, 562, 576, 590, 604, 619, 634, 649, 665, 681, 698, 715, 732, 750, 768, 787, 806, 825, 845, 866, 887, 909, 931, 953, 976])
E192 =ESeries(:E192, [100, 101, 102, 104, 105, 106, 107, 109, 110, 111, 113, 114, 115, 117, 118, 120, 121, 123, 124, 126, 127, 129, 130, 132, 133, 135, 137, 138, 140, 142, 143, 145, 147, 149, 150, 152, 154, 156, 158, 160, 162, 164, 165, 167, 169, 172, 174, 176, 178, 180, 182, 184, 187, 189, 191, 193, 196, 198, 200, 203, 205, 208, 210, 213, 215, 218, 221, 223, 226, 229, 232, 234, 237, 240, 243, 246, 249, 252, 255, 258, 261, 264, 267, 271, 274, 277, 280, 284, 287, 291, 294, 298, 301, 305, 309, 312, 316, 320, 324, 328, 332, 336, 340, 344, 348, 352, 357, 361, 365, 370, 374, 379, 383, 388, 392, 397, 402, 407, 412, 417, 422, 427, 432, 437, 442, 448, 453, 459, 464, 470, 475, 481, 487, 493, 499, 505, 511, 517, 523, 530, 536, 542, 549, 556, 562, 569, 576, 583, 590, 597, 604, 612, 619, 626, 634, 642, 649, 657, 665, 673, 681, 690, 698, 706, 715, 723, 732, 741, 750, 759, 768, 777, 787, 796, 806, 816, 825, 835, 845, 856, 866, 876, 887, 898, 909, 920, 931, 942, 953, 965, 976, 988])

function round(series::ESeries, input::Number, format::Union{Symbol, Bool}=false)
    input_normed, OOM = norm_to_between_100_and_1000(input)
    if input_normed>maximum(series.vals)
        candidate1, candidate2 = maximum(series.vals), 1000
    else
        candidate1, candidate2 = find_under_and_over(input_normed, series.vals)
    end
    geom_mean = geometric_mean(candidate1, candidate2)
    @assert (candidate1 ≤ input_normed ≤ candidate2) "Input value not in the range of rounding candidates. Something has gone wrong"

    result = (input_normed ≥ geom_mean ? candidate2  :  candidate1)
    result = result*10^OOM

    if format == false
        return round(result, sigdigits=3)
    else
        return formatted(result, format; ndigits=3)
    end
end

function round(series::ESeries, inputs::AbstractArray, format::Union{Symbol, Bool}=false)
    return [round(series, inputs[i], format) for i in eachindex(inputs)]  .|> identity   # identity call constrains the type
end

function print(series::ESeries)
    println("Values in $(string(series.name)):")
    for i in 1:3:length(series.vals)
        print(series.vals[i])
        print("   |   ")
        print(series.vals[i+1])
        print("   |   ")
        println(series.vals[i+2])
    end
end

function show(io::IO, series::ESeries)
    println(io, "Values in $(string(series.name)):")
    for i in 1:3:length(series.vals)
        print(io, series.vals[i])
        print(io, "   |   ")
        print(io, series.vals[i+1])
        print(io, "   |   ")
        println(io, series.vals[i+2])
    end
end

"""
    `determine_E(inputs...)`

Determine all E-series that contain each of the given inputs.
The inputs are comma-seperated.
The outputs are a vector of the names of possible series as `Symbol`s.

# Examples

```julia-repl
julia> determine_E(220, 470, 680)
3-element Vector{Symbol}:
 :E6
 :E12
 :E24

julia> determine_E(220, 470, 910)
1-element Vector{Symbol}:
 :E24
```
"""
function determine_E(inputs::Number...)
    inputs = [inputs...]  # Vector → tuple
    inputs = [norm_to_between_100_and_1000(i, return_OOM=false) for i in inputs]

    candidates = [E3, E6, E12, E24, E48, E96, E192]
    myfilt = fill(true, length(candidates))
    for input in inputs, i in eachindex(candidates)
        if input ∉ candidates[i].vals
            myfilt[i] = false
        end
    end
    possible_series = getfield.(candidates[myfilt], :name)
    #println("E-series containing all inputs:")
    #println.(possible_series)
    return possible_series
end

end


text:
using ReTest, Random, AdvancedHMC, ForwardDiff, AbstractMCMC
using Statistics: mean
include("common.jl")

@testset "AbstractMCMC w/ gdemo" begin
    rng = MersenneTwister(0)

    n_samples = 5_000
    n_adapts = 5_000

    θ_init = randn(rng, 2)

    model = AdvancedHMC.DifferentiableDensityModel(ℓπ_gdemo, ForwardDiff)
    init_eps = Leapfrog(1e-3)
    κ = NUTS(init_eps)
    metric = DiagEuclideanMetric(2)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, κ.τ.integrator))

    samples = AbstractMCMC.sample(
        rng, model, κ, metric, adaptor, n_adapts + n_samples;
        nadapts = n_adapts,
        init_params = θ_init,
        progress=false,
        verbose=false
    );

    # Transform back to original space.
    # NOTE: We're not correcting for the `logabsdetjac` here since, but
    # we're only interested in the mean it doesn't matter.
    for t in samples
        t.z.θ .= invlink_gdemo(t.z.θ)
    end
    m_est = mean(samples[n_adapts + 1:end]) do t
        t.z.θ
    end

    @test m_est ≈ [49 / 24, 7 / 6] atol=RNDATOL

    # Test that using the same AbstractRNG results in the same chain
    rng1 = MersenneTwister(42)
    rng2 = MersenneTwister(42)
    samples1 = AbstractMCMC.sample(
        rng1, model, κ, metric, adaptor, 10;
        nadapts = 0,
        progress=false,
        verbose=false
    );
    samples2 = AbstractMCMC.sample(
        rng2, model, κ, metric, adaptor, 10;
        nadapts = 0,
        progress=false,
        verbose=false
    );
    @test mapreduce(*, samples1, samples2) do s1, s2
        s1.z.θ == s2.z.θ
    end # Equivalent to using all, check that all samples are equal
end


text:
# This file provides tools for design of multi-robot target tracking experiments

using JLD2
using Printf
using Random
using Base.Threads
using Base.Iterators

export target_tracking_experiment, target_tracking_instance,
  iterate_target_tracking!, iterate_target_coverage!, visualize_experiment,
  visualize_time_step, run_experiments

# Run a target tracking experiment for a given number of steps
# By default, keep all data from each step
function target_tracking_experiment(;steps,
                                    num_robots,
                                    num_targets=default_num_targets(num_robots),
                                    configs::MultiRobotTargetTrackingConfigs,
                                    solver=solve_sequential,
                                    step_callback=copy_data
                                   )
  initialization = target_tracking_instance(num_robots=num_robots,
                                            num_targets=num_targets,
                                            configs=configs)

  robot_states = initialization.robot_states
  target_states = initialization.target_states
  target_filters = initialization.target_filters

  # Run target tracking. Store data in whatever the callback returns
  map(1:steps) do _
    # Simulate one step of tracking; update everything in place
    data = iterate_target_tracking!(robot_states=robot_states,
                                    target_states=target_states,
                                    target_filters=target_filters,
                                    configs=configs,
                                    solver=solver)

    step_callback(data)
  end
end

function target_tracking_instance(;num_robots::Int64,
                                  num_targets=default_num_targets(num_robots),
                                  configs::MultiRobotTargetTrackingConfigs)
  grid = configs.grid

  target_states = map(x->random_state(grid), 1:num_targets)

  (
   robot_states = map(x->random_state(grid), 1:num_robots),
   target_states = target_states,
   target_filters = map(x->Filter(grid, x), target_states)
  )
end

function iterate_target_tracking!(;robot_states::Vector{State},
                                  target_states::Vector{State},
                                  target_filters::Vector{<:AbstractFilter},
                                  configs::MultiRobotTargetTrackingConfigs,
                                  solver=solve_sequential)
  #
  # Update robot and target states
  #

  problem = MultiRobotTargetTrackingProblem(robot_states, target_filters,
                                            configs)
  solution = solver(problem)

  for (index, trajectory) in solution.elements
    robot_states[index] = trajectory[1]
  end
  trajectories = map(last, solution.elements)

  # Update Target states
  for (ii, target_state) in enumerate(target_states)
    target_states[ii] = target_dynamics(configs.grid, target_state)
  end

  #
  # Run filter updates (corresponding to new target states) and provide robots
  # with observations
  #

  # Process update
  for filter in target_filters
    process_update!(filter, transition_matrix(configs.grid))
  end

  # Sample ranging observations
  # (output is an array of arrays of robots' observations)
  range_observations = map(robot_states) do robot
    map(target_states) do target
      generate_observation(configs.grid, configs.sensor, robot, target)
    end
  end

  # Measurement update
  for (robot, observations) in zip(robot_states, range_observations)
    for (filter, observation) in zip(target_filters, observations)
      measurement_update!(filter, robot, get_states(configs.grid, filter),
                          configs.sensor, configs.grid, observation)
    end
  end

  # Note: by default, we do not copy the filters
  (
   robot_states=robot_states,
   target_states=target_states,
   target_filters=target_filters,
   trajectories=trajectories,
   range_observations=range_observations,
   objective=solution.value
  )
end

function iterate_target_coverage!(;robot_states::Vector{State},
                                  target_states::Vector{State},
                                  configs::MultiRobotTargetCoverageConfigs,
                                  solver=solve_sequential)
  #
  # Update robot and target states
  #

  problem = MultiRobotTargetCoverageProblem(robot_states, target_states,
                                            configs)
  solution = solver(problem)

  for (index, trajectory) in solution.elements
    robot_states[index] = trajectory[1]
  end
  trajectories = map(last, solution.elements)

  # Update Target states
  for (ii, target_state) in enumerate(target_states)
    target_states[ii] = target_dynamics(configs.grid, target_state)
  end

  # Note: by default, we do not copy the filters
  (
   robot_states=robot_states,
   target_states=target_states,
   trajectories=trajectories,
   objective=solution.value
  )
end

# The default updater modifies the states and filters in place
function copy_data(x)
  (
   robot_states=Array(x.robot_states),
   target_states=Array(x.target_states),
   target_filters=map(duplicate, x.target_filters),
   trajectories=x.trajectories,
   range_observations=x.range_observations,
   objective=x.objective
  )
end

#
# Visualize experiments
#

function visualize_experiment(data::Vector,
                              configs)
  if length(data) == 0
    println("There is no experiment data to visualize")
    return
  end

  grid = configs.grid

  figure()

  plot_state_space(grid)
  xlim([0, grid.width+1])
  ylim([0, grid.height+1])

  for (ii, step) in enumerate(data)
    println("Step: ", ii)

    plots = visualize_time_step(robot_states=step.robot_states,
                                target_states=step.target_states,
                                target_filters=step.target_filters,
                                trajectories=step.trajectories)

    # End of step (take input and clear the canvas

    line = readline()
    if line == "q"
      break
    end
    foreach(x->x.remove(), plots)
  end
end

function visualize_time_step(;robot_states,
                             target_states,
                             target_filters,
                             trajectories)
  plots=[]

  for robot in robot_states
    append!(plots, plot_quadrotor(robot, color=:blue, scale=0.3))
  end

  for trajectory in trajectories
    append!(plots, plot_states(trajectory, color=:blue, linestyle=":"))
  end

  for target in target_states
    append!(plots, plot_target(target))
  end

  append!(plots, visualize_filters(target_filters))

  plots
end

# Run a number of experiments with helpful features for batch processing
# * experiments run in parallel
# * individual results get cached
# * batch gets cached
# * provides helpful output as trials run
function run_experiments(tests;
                         trial_fun::Function,
                         print_summary::Function = x->nothing,
                         experiment_name::String,
                         data_folder::String,
                         reprocess=false,
                         threaded=true
                        )
  # produce file names
  data_file = string(data_folder, "/", experiment_name, ".jld2")
  trial_folder = string(data_folder, "/", experiment_name)
  get_trial_file = x->string(trial_folder, "/", experiment_name, " ", x, ".jld2")

  results = Dict{Any,Any}(key=>nothing for key in tests)

  if !isfile(data_file) || reprocess
    if !isdir(trial_folder)
      mkdir(trial_folder)
    end

    num_tests_completed = Atomic{Int64}(0)

    # First process any data that has been saved
    remaining_tests = filter(collect(tests)) do trial_spec
      run_test = true

      trial_file = get_trial_file(trial_spec)

      # Cache data from each trial
      if !reprocess && isfile(trial_file)
        try
          print("Loading: ", trial_file, "...")

          @load trial_file trial_results
          println("Loaded")

          results[trial_spec] = trial_results

          run_test = false
        catch e
          println("Failed to load ", trial_spec)
        end
      end

      run_test
    end

    println("\n", length(remaining_tests), " tests remain\n")

    # Mutex for load/save/output
    load_save_lock = ReentrantLock()

    start = time()
    shuffle!(remaining_tests)
    trial_backtrace = nothing

    spawn_for_each(remaining_tests, threaded=threaded) do trial_spec
      id = threadid()

      lock(load_save_lock)
      println("Thread-", id, " running: ", trial_spec)
      unlock(load_save_lock)

      trial_start = time()

      # Run trial
      trial_results = trial_fun(trial_spec)

      # Save results
      try
        lock(load_save_lock)
        print(id, "-Saving: ", trial_spec, "...")

        trial_file = get_trial_file(trial_spec)
        @save trial_file trial_results

        println(id, "-Saved")
      catch e
        # Save the error for later
        trial_backtrace = catch_backtrace()

        println(id, "-Failed to save ", trial_spec)
        return
      finally
        unlock(load_save_lock)
      end

      # Store value
      results[trial_spec] = trial_results

      # Summarize completion status
      elapsed = time() - start
      trial_elapsed = time() - trial_start

      # (returns old value)
      completed = atomic_add!(num_tests_completed, 1) + 1
      completion = completed/length(remaining_tests)
      projected = elapsed / completion
      hour = 3600

      lock(load_save_lock)
      print_summary(trial_spec)
      println(" (Done, ", completed, "/", length(remaining_tests),
              @sprintf(" %0.2f", 100completion), "%",
              " Trial time: ", @sprintf("%0.0fs", trial_elapsed),
              " Elapsed: ", @sprintf("%0.2fh", elapsed / hour),
              " Total: ", @sprintf("%0.2fh", projected / hour),
              " Remaining: ", @sprintf("%0.2fh", (projected - elapsed) / hour),
              ")"
             )
      unlock(load_save_lock)
    end

    if !isnothing(trial_backtrace)
      println("Loading or producing results produce errors. Backtrace:")
      map(println, stacktrace(trial_backtrace))
      println()

      error("One or more trials failed.")
    end

    @save data_file results
  else
    @load data_file results
  end

  results
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using UlamSpiral
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@testset "position function" begin
    @test UlamSpiral.position(1) == (0,0)
    @test UlamSpiral.position(10,10) == (0,0)
    @test UlamSpiral.position(10) == (2, -1)

    @test_throws DomainError UlamSpiral.position(0)
    @test_throws DomainError UlamSpiral.position(9, 10)
end

#
# at present these tests are a bit trivial, because most of this
# is about plotting. Perhaps can use
#     https://github.com/JuliaPlots/VisualRegressionTests.jl
# to do this, but perhaps is overkill for a niche package like this.
#