text:
#=
 TODOs:

    - get quadratic objective functions

=#

__precompile__()
module LinQuadOptInterface

using MathOptInterface
using MathOptInterface.Utilities

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

# functions
const Linear = MOI.ScalarAffineFunction{Float64}
const Quad   = MOI.ScalarQuadraticFunction{Float64}
const SinVar = MOI.SingleVariable
const VecVar = MOI.VectorOfVariables
const VecLin = MOI.VectorAffineFunction{Float64}
# sets
const LE = MOI.LessThan{Float64}
const GE = MOI.GreaterThan{Float64}
const EQ = MOI.EqualTo{Float64}
const IV = MOI.Interval{Float64}
const SOS1 = MOI.SOS1{Float64}
const SOS2 = MOI.SOS2{Float64}
# constraint references
const CI{F,S} = MOI.ConstraintIndex{F,S}
const LCI{S} = CI{Linear,S}
const VLCI{S} = CI{VecLin,S}
const QCI{S} = CI{Quad,S}
const SVCI{S} = CI{SinVar,S}
const VVCI{S} = CI{VecVar,S}
# variable reference
const VarInd = MOI.VariableIndex

const LinSets = Union{LE, GE, EQ, IV}
const VecLinSets = Union{MOI.Nonnegatives, MOI.Nonpositives, MOI.Zeros}

struct ConstraintMapping
    # rows in constraint matrix
    less_than::Dict{LCI{LE}, Int}
    greater_than::Dict{LCI{GE}, Int}
    equal_to::Dict{LCI{EQ}, Int}
    interval::Dict{LCI{IV}, Int}

    # vectors of rows in constraint matrix
    nonnegatives::Dict{VLCI{MOI.Nonnegatives}, Vector{Int}}
    nonpositives::Dict{VLCI{MOI.Nonpositives}, Vector{Int}}
    zeros::Dict{VLCI{MOI.Zeros}, Vector{Int}}

    # rows in quadratic constraint matrix
    q_less_than::Dict{QCI{LE}, Int}
    q_greater_than::Dict{QCI{GE}, Int}
    q_equal_to::Dict{QCI{EQ}, Int}

    # references to variable
    upper_bound::Dict{SVCI{LE}, VarInd}
    lower_bound::Dict{SVCI{GE}, VarInd}
    fixed_bound::Dict{SVCI{EQ}, VarInd}
    interval_bound::Dict{SVCI{IV}, VarInd}

    # vectors of rows in constraint matrix
    vv_nonnegatives::Dict{VVCI{MOI.Nonnegatives}, Vector{Int}}
    vv_nonpositives::Dict{VVCI{MOI.Nonpositives}, Vector{Int}}
    vv_zeros::Dict{VVCI{MOI.Zeros}, Vector{Int}}

    integer::Dict{SVCI{MOI.Integer}, VarInd}
    #=
     for some reason CPLEX doesn't respect bounds on a binary variable, so we
     should store the previous bounds so that if we delete the binary constraint
     we can revert to the old bounds
    =#
    binary::Dict{SVCI{MOI.ZeroOne}, Tuple{VarInd, Float64, Float64}}
    sos1::Dict{VVCI{SOS1}, Int}
    sos2::Dict{VVCI{SOS2}, Int}
    semicontinuous::Dict{SVCI{MOI.Semicontinuous{Float64}}, VarInd}
    semiinteger::Dict{SVCI{MOI.Semiinteger{Float64}}, VarInd}

end
ConstraintMapping() = ConstraintMapping(
    Dict{LCI{LE}, Int}(),
    Dict{LCI{GE}, Int}(),
    Dict{LCI{EQ}, Int}(),
    Dict{LCI{IV}, Int}(),
    Dict{VLCI{MOI.Nonnegatives}, Vector{Int}}(),
    Dict{VLCI{MOI.Nonpositives}, Vector{Int}}(),
    Dict{VLCI{MOI.Zeros}, Vector{Int}}(),
    Dict{QCI{LE}, Int}(),
    Dict{QCI{GE}, Int}(),
    Dict{QCI{EQ}, Int}(),
    Dict{SVCI{LE}, VarInd}(),
    Dict{SVCI{GE}, VarInd}(),
    Dict{SVCI{EQ}, VarInd}(),
    Dict{SVCI{IV}, VarInd}(),
    Dict{VVCI{MOI.Nonnegatives}, Vector{Int}}(),
    Dict{VVCI{MOI.Nonpositives}, Vector{Int}}(),
    Dict{VVCI{MOI.Zeros}, Vector{Int}}(),
    Dict{SVCI{MOI.Integer}, VarInd}(),
    Dict{SVCI{MOI.ZeroOne}, Tuple{VarInd, Float64, Float64}}(),
    Dict{VVCI{SOS1}, Int}(),
    Dict{VVCI{SOS2}, Int}(),
    Dict{SVCI{MOI.Semicontinuous{Float64}}, VarInd}(),
    Dict{SVCI{MOI.Semiinteger{Float64}}, VarInd}()
)

"""
    shift_references_after_delete_affine!(m, row)

This function updates all of the references in `m`
after we have deleted row `row` in the affine constraint matrix.
"""
function shift_references_after_delete_affine!(m, row)
    for scalar_affine in [
            cmap(m).less_than,
            cmap(m).greater_than,
            cmap(m).equal_to,
            cmap(m).interval
        ]
        for (key, val) in scalar_affine
            if val > row
                scalar_affine[key] -= 1
            end
        end
    end

    for vector_affine in [
            cmap(m).nonnegatives,
            cmap(m).nonpositives,
            cmap(m).zeros,
            cmap(m).vv_nonnegatives,
            cmap(m).vv_nonpositives,
            cmap(m).vv_zeros
        ]
        for (key, vals) in vector_affine
            for (i, val) in enumerate(vals)
                if val > row
                    vector_affine[key][i] -= 1
                end
            end
        end
    end
end

"""
    shift_references_after_delete_quadratic!(m, row)

This function updates all of the references in `m`
after we have deleted row `row` in the quadratic constraint matrix.
"""
function shift_references_after_delete_quadratic!(m, row)
    for scalar_quadratic in [
            cmap(m).q_less_than,
            cmap(m).q_greater_than,
            cmap(m).q_equal_to
        ]
        for (key, val) in scalar_quadratic
            if val > row
                scalar_quadratic[key] -= 1
            end
        end
    end
end

function Base.isempty(map::ConstraintMapping)

    ret = true
    ret = ret && isempty(map.less_than)
    ret = ret && isempty(map.greater_than)
    ret = ret && isempty(map.equal_to)
    ret = ret && isempty(map.interval)
    ret = ret && isempty(map.nonnegatives)
    ret = ret && isempty(map.nonpositives)
    ret = ret && isempty(map.zeros)
    ret = ret && isempty(map.q_greater_than)
    ret = ret && isempty(map.q_greater_than)
    ret = ret && isempty(map.q_equal_to)
    ret = ret && isempty(map.upper_bound)
    ret = ret && isempty(map.lower_bound)
    ret = ret && isempty(map.fixed_bound)
    ret = ret && isempty(map.interval_bound)
    ret = ret && isempty(map.vv_nonnegatives)
    ret = ret && isempty(map.vv_nonpositives)
    ret = ret && isempty(map.vv_zeros)
    ret = ret && isempty(map.integer)
    ret = ret && isempty(map.binary)
    ret = ret && isempty(map.sos1)
    ret = ret && isempty(map.sos2)
    ret = ret && isempty(map.semiinteger)
    ret = ret && isempty(map.semicontinuous)

    return ret
end

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

@enum(ObjectiveType,
      SingleVariableObjective,
      AffineObjective,
      QuadraticObjective)

# Abstract + macro
abstract type LinQuadOptimizer <: MOI.AbstractOptimizer end
@def LinQuadOptimizerBase begin

    inner#::LinQuadOptInterface.LinQuadOptimizer

    name::String

    obj_type::LinQuadOptInterface.ObjectiveType
    single_obj_var::Nullable{MOI.VariableIndex}
    obj_sense::MOI.OptimizationSense

    last_variable_reference::UInt64
    variable_mapping::Dict{MOI.VariableIndex, Int}
    variable_names::Dict{MOI.VariableIndex, String}
    variable_names_rev::Dict{String, MOI.VariableIndex}
    variable_references::Vector{MOI.VariableIndex}

    variable_primal_solution::Vector{Float64}
    variable_dual_solution::Vector{Float64}

    last_constraint_reference::UInt64
    constraint_mapping::LinQuadOptInterface.ConstraintMapping

    constraint_constant::Vector{Float64}
    constraint_primal_solution::Vector{Float64}
    constraint_dual_solution::Vector{Float64}

    qconstraint_primal_solution::Vector{Float64}
    qconstraint_dual_solution::Vector{Float64}

    # TODO(odow): temp hack for constraint names
    constraint_names::Dict{Any, String}
    constraint_names_rev::Dict{String, Any}

    objective_constant::Float64

    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    dual_status::MOI.ResultStatusCode
    primal_result_count::Int
    dual_result_count::Int

    solvetime::Float64
end

function MOI.isempty(m::LinQuadOptimizer)

    ret = true
    ret = ret && m.name == ""
    ret = ret && m.obj_type == AffineObjective
    ret = ret && isnull(m.single_obj_var)
    ret = ret && m.obj_sense == MOI.MinSense
    ret = ret && m.last_variable_reference == 0
    ret = ret && isempty(m.variable_mapping)
    ret = ret && isempty(m.variable_names)
    ret = ret && isempty(m.variable_names_rev)
    ret = ret && isempty(m.variable_references)
    ret = ret && isempty(m.variable_primal_solution)
    ret = ret && isempty(m.variable_dual_solution)
    ret = ret && m.last_constraint_reference == 0
    ret = ret && isempty(m.constraint_mapping)
    ret = ret && isempty(m.constraint_constant)
    ret = ret && isempty(m.constraint_primal_solution)
    ret = ret && isempty(m.constraint_dual_solution)
    ret = ret && isempty(m.qconstraint_primal_solution)
    ret = ret && isempty(m.qconstraint_dual_solution)
    ret = ret && isempty(m.constraint_names)
    ret = ret && isempty(m.constraint_names_rev)
    ret = ret && m.objective_constant == 0.0
    ret = ret && m.termination_status == MOI.OtherError
    ret = ret && m.primal_status == MOI.UnknownResultStatus
    ret = ret && m.dual_status == MOI.UnknownResultStatus
    ret = ret && m.primal_result_count == 0
    ret = ret && m.dual_result_count == 0
    ret = ret && m.solvetime == 0.0

    return ret
end
function MOI.empty!(m::M, env = nothing) where M<:LinQuadOptimizer
    m.name = ""
    m.inner = LinearQuadraticModel(M,env)

    m.obj_type = AffineObjective
    m.single_obj_var = nothing
    # we assume the default is minimization
    m.obj_sense = MOI.MinSense

    m.last_variable_reference = 0
    m.variable_mapping = Dict{MathOptInterface.VariableIndex, Int}()
    m.variable_names = Dict{MathOptInterface.VariableIndex, String}()
    m.variable_names_rev = Dict{String, MathOptInterface.VariableIndex}()
    m.variable_references = MathOptInterface.VariableIndex[]

    m.variable_primal_solution = Float64[]
    m.variable_dual_solution = Float64[]

    m.last_constraint_reference = 0
    m.constraint_mapping = LinQuadOptInterface.ConstraintMapping()

    m.constraint_constant = Float64[]
    m.constraint_primal_solution = Float64[]
    m.constraint_dual_solution = Float64[]

    m.qconstraint_primal_solution = Float64[]
    m.qconstraint_dual_solution = Float64[]

    m.constraint_names = Dict{Any, String}()
    m.constraint_names_rev = Dict{String, Any}()

    m.objective_constant = 0.0

    m.termination_status = MathOptInterface.OtherError
    m.primal_status = MathOptInterface.UnknownResultStatus
    m.dual_status = MathOptInterface.UnknownResultStatus
    m.primal_result_count = 0
    m.dual_result_count = 0

    m.solvetime = 0.0

    nothing
end

function MOI.get(m::LinQuadOptimizer, ::MOI.Name)
    m.name
end
MOI.canget(m::LinQuadOptimizer, ::MOI.Name) = true
function MOI.set!(m::LinQuadOptimizer, ::MOI.Name, name::String)
    m.name = name
end
MOI.canset(m::LinQuadOptimizer, ::MOI.Name) = true

function MOI.supportsconstraint(m::LinQuadOptimizer, ft::Type{F}, st::Type{S}) where F <: MOI.AbstractFunction where S <: MOI.AbstractSet
    (ft,st) in supported_constraints(m)
end
function MOI.supports(m::LinQuadOptimizer, ::MOI.ObjectiveFunction{F}) where F <: MOI.AbstractFunction
    F in supported_objectives(m)
end

# a useful helper function
function deleteref!(dict::Dict, i::Int, ref)
    for (key, val) in dict
        if val > i
            dict[key] -= 1
        end
    end
    delete!(dict, ref)
end

include("variables.jl")
include("constraints.jl")
include("objective.jl")
include("solve.jl")
include("copy.jl")

include("solver_interface.jl")

include("mockoptimizer.jl")

end


text:
export plotcost_pernbsim, plotcost_persquare, plotcost_suboptimal

using Plots

include("markovDecision.jl")
include("simulate.jl")

function simulatedcost(board, dice, circular, start_box, nb_iter)
    allsim = Vector{Int64}()
    for i=1:nb_iter
    	push!(allsim, simulate_k(board, dice, circular, start_box))
    end
    return mean(allsim)
end

"""
    plotcost_pernbsim(board)

Plot both the simulated and expected cost for multiple number of simulations
of a game starting at square 1.
"""
function plotcost_pernbsim(board)
    fig1 = plot()
    # define board
    circular = true
    # expected cost
    (expec, dice) = markovdecision(board, circular)
    expected_cost_from_1 = expec[1]

    # cost plot
    start_box = 1
    powers = [i for i=1:7]
    nb_iter = [10^i for i in powers]
    y = [simulatedcost(board, dice, circular, start_box, x) for x in nb_iter]

    plot!(powers, y, marker="o",
        label=["Simulated cost"],
        xlabel=("Log 10 number of simulations"),
        ylabel=("Cost"))
    hline!([expected_cost_from_1]; label="Expected cost")

    # error plot
    fig2 = plot()
    error = [abs(expected_cost_from_1 - simulated_cost) for simulated_cost in y]/expected_cost_from_1
    plot!(powers, error, marker="o",
        label=["Error"],
        xlabel=("Log 10 number of simulations"),
        ylabel=("Relative error"))

    return fig1, fig2
end

"""
    plotcost_persquare(board)

Plot both the simulated and expected cost for starting the game at every square
(one square is one x-axis point).
This allows to check whether each expected cost corresponds to the simulated one,
when the game is repeated multiple times.
"""
function plotcost_persquare(board)
    fig = plot()
    # define board
    circular = true
    # expected cost
    (expec, dice) = markovdecision(board, circular)

    nb_iter = 100
    squares = [i for i=1:15]
    y = [simulatedcost(board, dice, circular, square, nb_iter) for square in squares]

    scatter!(squares, y,
        color=("blue"),
        label=["Simulated cost"]
        )

    scatter!(squares, expec,
        color=("red"),
        label=["Expected cost"],
        xlabel=("Square number"),
        ylabel=("Cost")
        )

    figname = "../../img/board_magic/cost_per_square_" * string(nb_iter) * "_iter_circ.pdf"
    return fig, figname
end

"""
    plotcost_suboptimal(board)

Plot the simulated cost for different suboptimal dice statgies
of a game starting at square 1.
"""
function plotcost_suboptimal(board)
    fig = plot()
    # define board
    circular = true
    # expected cost
    (expec, dice) = markovdecision(board, circular)
    expected_cost_from_1 = expec[1]

    # cost plot
    start_box = 1
    powers = [i for i=1:3]
    nb_iter = [10^i for i in powers]
    y_sim = [simulatedcost(board, dice, circular, start_box, x) for x in nb_iter]
    y_subopt1 = [simulatedcost(board, ones(Int64, 15), circular, start_box, x) for x in nb_iter]
    y_subopt2 = [simulatedcost(board, 2*ones(Int64, 15), circular, start_box, x) for x in nb_iter]
    y_subopt3 = [simulatedcost(board, 3*ones(Int64, 15), circular, start_box, x) for x in nb_iter]
    dice_subopt4 = ones(Int64, 15)
    dice_subopt4[1:9] = [t == 0 ? 3 : (t == 1 ? 1 : 2) for t in board[2:10]]
    dice_subopt4[11:13] = [t == 0 ? 3 : (t == 1 ? 1 : 2) for t in board[12:14]]
    y_subopt4 = [simulatedcost(board, dice_subopt4, circular, start_box, x) for x in nb_iter]

    plot!(powers, y_sim, marker="o",
        label=["Optimal cost"],
        xlabel=("Log 10 number of simulations"),
        ylabel=("Cost"))
    plot!(powers, y_subopt1, marker="o",
        label=["Cost for only security dice"])
    plot!(powers, y_subopt2, marker="o",
        label=["Cost for only normal dice"])
    plot!(powers, y_subopt3, marker="o",
        label=["Cost for only risky dice"])
    plot!(powers, y_subopt4, marker="o",
        label=["Cost for custom strategy"])

    return fig
end

board_unif_low = [0,1,0,0,2,0,1,0,2,0,0,1,0,2,0]
board_right_high = [0,0,0,0,0,2,0,0,1,1,2,2,0,0,0]
board_left_high = [0,1,2,1,2,2,1,0,0,0,0,0,0,0,0]
board_magic = [0,3,0,3,0,0,0,0,3,3,0,0,0,0,0]
# #
# plot_cost, plot_errors = plotcost_pernbsim(board_magic)
# savefig(plot_cost, "../../img/board_magic/cost_iterations_log_circ.pdf")
# savefig(plot_errors, "../../img/board_magic/error_iterations_log_circ.pdf")

# #
# p, figname = plotcost_persquare(board_magic)
# savefig(p, figname)

# #
# plot_subopt = plotcost_suboptimal(board_magic)
# savefig(plot_subopt, "../../img/board_magic/cost_subopt_log_circ.pdf")

ext: jl

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
struct Model{F, L}
    f::F
    layers::L
end

@functor Model (layers,)

(m::Model)(x) = m.f(x, m.layers)


text:
σ{T <: Real}(z::T) = 1.0/(1.0 + exp(-z))

∇σ{T <: Real}(z::T) = 1.0 / (2.0 + exp(z) + exp(-z))

max_stars_repo_name: helgith/ExampleNeuralNetworks

max_issues_repo_name: helgith/ExampleNeuralNetworks

max_forks_repo_name: helgith/ExampleNeuralNetworks


text:
frag_shader = shader_file("rectangle.frag")

@testset "Shader cache" begin
    spec = ShaderSpecification(frag_shader, GLSL)
    cache = Lava.ShaderCache(device)
    Shader(Lava.ShaderCache(device), spec) # trigger JIT compilation
    t = @elapsed Shader(cache, spec)
    @test t > 0.01
    t = @elapsed Shader(cache, spec)
    @test t < 1e-5
end


text:
"""
    getindex!(cache,a::AbstractArray,i...)

Returns the item of the array `a` associated with index `i`
by (possibly) using the scratch data passed in the `cache` object.

It defaults to

    getindex!(cache,a::AbstractArray,i...) = a[i...]

# Examples

Iterating over an array using the `getindex!` function

```jldoctest
using Gridap.Arrays

a = collect(10:15)

cache = array_cache(a)
for i in eachindex(a)
  ai = getindex!(cache,a,i)
  println("\$i -> \$ai")
end

# output
1 -> 10
2 -> 11
3 -> 12
4 -> 13
5 -> 14
6 -> 15
```

"""
getindex!(cache,a::AbstractArray,i...) = a[i...]

"""
    array_cache(a::AbstractArray)

Returns a cache object to be used in the [`getindex!`](@ref) function.
It defaults to

    array_cache(a::T) where T = nothing

for types `T` such that `uses_hash(T) == Val(false)`, and

    function array_cache(a::T) where T
      hash = Dict{UInt,Any}()
      array_cache(hash,a)
    end

for types `T` such that `uses_hash(T) == Val(true)`, see the [`uses_hash`](@ref) function. In the later case, the
type `T` should implement the following signature:

    array_cache(hash::Dict,a::AbstractArray)

where we pass a dictionary (i.e., a hash table) in the first argument. This hash table can be used to test
if the object `a` has already built a cache and re-use it as follows

    id = objectid(a)
    if haskey(hash,id)
      cache = hash[id] # Reuse cache
    else
      cache = ... # Build a new cache depending on your needs
      hash[id] = cache # Register the cache in the hash table
    end

This mechanism is needed, e.g., to re-use intermediate results in complex lazy operation trees.
In multi-threading computations, a different hash table per thread has to be used in order
to avoid race conditions.
"""
function array_cache(a::AbstractArray)
  _default_array_cache(a,uses_hash(a))
end

function array_cache(hash,a::T) where T
  if uses_hash(T) == Val{true}()
    error("array_cache(::Dict,::$T) not defined")
  end
  array_cache(a)
end

function _default_array_cache(a,::Val{false})
  nothing
end

function _default_array_cache(a,::Val{true})
  hash = Dict{UInt,Any}()
  array_cache(hash,a)
end

"""
    uses_hash(::Type{<:AbstractArray})

This function is used to specify if the type `T` uses the
hash-based mechanism to reuse caches.  It should return
either `Val(true)` or `Val(false)`. It defaults to

    uses_hash(::Type{<:AbstractArray}) = Val(false)

Once this function is defined for the type `T` it can also
be called on instances of `T`.

"""
uses_hash(::Type{<:AbstractArray}) = Val(false)

uses_hash(::T) where T = uses_hash(T)

"""
$(TYPEDSIGNATURES)

Returns an arbitrary instance of `eltype(a)`. The default returned value is the first entry
in the array if `length(a)>0` and `testvalue(eltype(a))` if `length(a)==0`
See the [`testvalue`](@ref) function.

# Examples

```jldoctest
using Gridap.Arrays

a = collect(3:10)
ai = testitem(a)

b = Int[]
bi = testitem(b)

(ai, bi)

# output
(3, 0)

```

"""
function testitem(a::AbstractArray{T}) where T
  if length(a) >0
    first(a)
  else
    testvalue(T)
  end::T
end

function testitem(a::Fill)
  a.value
end

"""
    get_array(a::AbstractArray)

Returns `a`.
"""
function get_array(a::AbstractArray)
  a
end

"""
"""
function get_arrays(a,b...)
  (get_array(a),get_arrays(b...)...)
end

function get_arrays(a)
  (get_array(a),)
end

# Test the interface

"""
    test_array(
      a::AbstractArray{T,N}, b::AbstractArray{S,N},cmp=(==)) where {T,S,N}

Checks if the entries in `a` and `b` are equal using the comparison function `cmp`.
It also stresses the new methods added to the `AbstractArray` interface.
"""
function test_array(
  a::AbstractArray{T,N}, b::AbstractArray{S,N},cmp=(==)) where {T,S,N}
  @test cmp(a,b)
  cache = array_cache(a)
  t = true
  for i in eachindex(a)
    bi = b[i]
    ai = getindex!(cache,a,i)
    t = t && cmp(bi,ai)
  end
  @test t
  t = true
  for i in eachindex(a)
    ai = getindex!(cache,a,i)
    t = t && (typeof(ai) <: eltype(a))
    t = t && (typeof(ai) <: T)
  end
  @test t
  @test IndexStyle(a) == IndexStyle(b)
  @test isa(testitem(a),eltype(a))
  if length(a) > 0
    @test testitem(a) == first(a)
  end
  true
end

# Some API

"""
    testitems(b::AbstractArray...) -> Tuple

Returns a tuple with the result of `testitem` applied to each of the
arrays in `b`.

# Examples

```jldoctest
using Gridap.Arrays

a = collect(3:10)
b = Int[]
c = Float64[]
d = ones(10)

testitems(a,b,c,d)

# output
(3, 0, 0.0, 1.0)

```
"""
function testitems(a::AbstractArray,b::AbstractArray...)
  va = testitem(a)
  vb = testitems(b...)
  (va,vb...)
end

function testitems(a::AbstractArray)
  va = testitem(a)
  (va,)
end

"""
    array_caches(a::AbstractArray...) -> Tuple

Returns a tuple with the cache of each array in `a`.
"""
function array_caches(a::AbstractArray,b::AbstractArray...)
  hash = Dict{UInt,Any}()
  array_caches(hash,a,b...)
end

function array_caches(hash::Dict,a::AbstractArray,b::AbstractArray...)
  ca = array_cache(hash,a)
  cb = array_caches(hash,b...)
  (ca,cb...)
end

function array_caches(hash::Dict,a::AbstractArray)
  ca = array_cache(hash,a)
  (ca,)
end

array_caches() = ()

"""
    getitems!(c::Tuple,a::Tuple,i...) -> Tuple

Extracts the `i`-th entry of all arrays in the tuple `a` using the caches in the tuple
`c`. The results is a tuple containing each one of the extracted entries.

# Example

Iterating over three different arrays simultaneously using `getitems!`

```jldoctest
using Gridap.Arrays

a = collect(0:5)
b = collect(10:15)
c = collect(20:25)

caches = array_caches(a,b,c)
for i in eachindex(a)
   s = getitems!(caches,(a,b,c),i)
   println("\$i -> \$s")
end

# output
1 -> (0, 10, 20)
2 -> (1, 11, 21)
3 -> (2, 12, 22)
4 -> (3, 13, 23)
5 -> (4, 14, 24)
6 -> (5, 15, 25)
```

"""
@inline function getitems!(cf::Tuple,a::Tuple{Vararg{<:AbstractArray}},i...)
  _getitems!(cf,i,a...)
end

getitems!(::Tuple{},::Tuple{},i) = ()

@inline function _getitems!(c,i,a,b...)
  ca,cb = _split(c...)
  ai = getindex!(ca,a,i...)
  bi = getitems!(cb,b,i...)
  (ai,bi...)
end

@inline function _getitems!(c,i,a)
  ca, = c
  ai = getindex!(ca,a,i...)
  (ai,)
end

# Hack to fix type-instability (use generated function?)
@inline function _getitems!(c,i,a1,a2)
  ca1,ca2 = c
  a1i = getindex!(ca1,a1,i...)
  a2i = getindex!(ca2,a2,i...)
  (a1i,a2i)
end

# Hack to fix type-instability (use generated function?)
@inline function _getitems!(c,i,a1,a2,a3)
  ca1,ca2,ca3 = c
  a1i = getindex!(ca1,a1,i...)
  a2i = getindex!(ca2,a2,i...)
  a3i = getindex!(ca3,a3,i...)
  (a1i,a2i,a3i)
end

# Hack to fix type-instability
@inline function _getitems!(c,i,a1,a2,a3,a4)
  ca1,ca2,ca3,ca4 = c
  a1i = getindex!(ca1,a1,i...)
  a2i = getindex!(ca2,a2,i...)
  a3i = getindex!(ca3,a3,i...)
  a4i = getindex!(ca4,a4,i...)
  (a1i,a2i,a3i,a4i)
end

@inline function _getitems!(c,i,a1,a2,a3,a4,a5)
  ca1,ca2,ca3,ca4,ca5 = c
  a1i = getindex!(ca1,a1,i...)
  a2i = getindex!(ca2,a2,i...)
  a3i = getindex!(ca3,a3,i...)
  a4i = getindex!(ca4,a4,i...)
  a5i = getindex!(ca5,a5,i...)
  (a1i,a2i,a3i,a4i,a5i)
end

@inline function _getitems!(c,i,a1,a2,a3,a4,a5,a6)
  ca1,ca2,ca3,ca4,ca5,ca6 = c
  a1i = getindex!(ca1,a1,i...)
  a2i = getindex!(ca2,a2,i...)
  a3i = getindex!(ca3,a3,i...)
  a4i = getindex!(ca4,a4,i...)
  a5i = getindex!(ca5,a5,i...)
  a6i = getindex!(ca6,a6,i...)
  (a1i,a2i,a3i,a4i,a5i,a6i)
end

@inline function _getitems!(c,i,a1,a2,a3,a4,a5,a6,a7)
  ca1,ca2,ca3,ca4,ca5,ca6,ca7 = c
  a1i = getindex!(ca1,a1,i...)
  a2i = getindex!(ca2,a2,i...)
  a3i = getindex!(ca3,a3,i...)
  a4i = getindex!(ca4,a4,i...)
  a5i = getindex!(ca5,a5,i...)
  a6i = getindex!(ca6,a6,i...)
  a7i = getindex!(ca7,a7,i...)
  (a1i,a2i,a3i,a4i,a5i,a6i,a7i)
end

"""
"""
@inline function getitems(a::Tuple{Vararg{<:AbstractArray}},i...)
  _getitems(i,a...)
end

@inline function _getitems(i,a,b...)
  ai = a[i...]
  bi = getitems(b,i...)
  (ai,bi...)
end

@inline function _getitems(i,a)
  ai = a[i...]
  (ai,)
end

"""
"""
function add_to_array!(a::AbstractArray{Ta,N},b::AbstractArray{Tb,N},combine=+) where {Ta,Tb,N}
  @assert size(a) == size(b) "Arrays sizes mismatch"
  @inbounds for i in eachindex(a)
    a[i] = combine(a[i],b[i])
  end
end

function add_to_array!(a::AbstractArray,b::Number,combine=+)
  @inbounds for i in eachindex(a)
    a[i] = combine(a[i],b)
  end
end


text:
import .PyPlot: PyPlot

function Base.convert(t::Type{PngImage}, figure::PyPlot.Figure)
    pb = PipeBuffer()
    show(pb, MIME("image/png"), figure)
    return PngImage(pb)
end

preprocess(name, plot::PyPlot.Figure, data) = preprocess(name, convert(PngImage, plot), data)
preprocess(name, plots::AbstractArray{<:PyPlot.Figure}, data) = begin
    for (i, plot)=enumerate(plots)
        preprocess(name*"/$i", plot, data)
    end
    return data
end


text:
const asm_syntax = Ref(:att)

function cthulhu_llvm(io::IO, mi, optimize, debuginfo, params)
    dump = InteractiveUtils._dump_function_linfo(
        mi, params.world, #=native=# false,
        #=wrapper=# false, #=strip_ir_metadata=# true,
        #=dump_module=# false, #=syntax=# asm_syntax[],
        optimize, debuginfo ? :source : :none)
    print(io, dump)
end

function cthulhu_native(io::IO, mi, optimize, debuginfo, params)
    dump = InteractiveUtils._dump_function_linfo(
        mi, params.world, #=native=# true,
        #=wrapper=# false, #=strip_ir_metadata=# true,
        #=dump_module=# false, #=syntax=# asm_syntax[],
        optimize, debuginfo ? :source : :none)
    print(io, dump)
end

cthulhu_warntype(args...) = cthulhu_warntype(stdout, args...)
function cthulhu_warntype(io::IO, src, rettype, debuginfo)
    if VERSION < v"1.1.0-DEV.762"
    elseif VERSION < v"1.2.0-DEV.229"
        lineprinter = Base.IRShow.debuginfo[debuginfo]
    else
        debuginfo = Base.IRShow.debuginfo(debuginfo)
        lineprinter = Base.IRShow.__debuginfo[debuginfo]
    end

    lambda_io::IOContext = stdout
    if src.slotnames !== nothing
        lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES =>  Base.sourceinfo_slotnames(src))
    end
    print(io, "Body")
    InteractiveUtils.warntype_type_printer(io, rettype, true)
    println(io)
    if VERSION < v"1.1.0-DEV.762"
        Base.IRShow.show_ir(lambda_io, src, InteractiveUtils.warntype_type_printer)
    else
        Base.IRShow.show_ir(lambda_io, src, lineprinter(src), InteractiveUtils.warntype_type_printer)
    end
    return nothing
end


text:
function get_iodc(obj::QZSSOrbitalElements)
    return jcall(obj, "getIODC", jint, ())
end

function get_iode(obj::QZSSOrbitalElements)
    return jcall(obj, "getIODE", jint, ())
end

function get_tgd(obj::QZSSOrbitalElements)
    return jcall(obj, "getTGD", jdouble, ())
end


text:
module ThreadPinning

# stdlibs
using Base.Threads: @threads, nthreads
using Libdl
using LinearAlgebra
using Random
using DelimitedFiles

# system information
Base.@kwdef struct SysInfo
    nsockets::Int = 1
    nnuma::Int = 1
    hyperthreading::Bool = false
    cpuids::Vector{Int} = collect(0:(Sys.CPU_THREADS - 1))
    cpuids_sockets::Vector{Vector{Int}} = [collect(0:(Sys.CPU_THREADS - 1))]
    cpuids_numa::Vector{Vector{Int}} = [collect(0:(Sys.CPU_THREADS - 1))]
    ishyperthread::Vector{Bool} = fill(false, Sys.CPU_THREADS)
end
function Base.show(io::IO, sysinfo::SysInfo)
    return print(io, "SysInfo()")
end
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, sysinfo::SysInfo)
    summary(io, sysinfo)
    println(io)
    fnames = fieldnames(SysInfo)
    for fname in fnames[1:(end - 1)]
        println(io, "├ $fname: ", getfield(sysinfo, fname))
    end
    print(io, "└ $(fnames[end]): ", getfield(sysinfo, fnames[end]))
    return nothing
end

# constants
const SYSINFO_ATTEMPT = Ref{Bool}(false) # have we yet attempted to gather the sysinfo
const SYSINFO_SUCCESS = Ref{Bool}(false) # have we succeeded in gathering the sysinfo yet
const SYSINFO = Ref{SysInfo}(SysInfo())

# includes
include("utility.jl")
include("lscpu_examples.jl")
include("libs/libc.jl")
include("libs/libuv.jl")
include("libs/libpthread.jl")
include("blas.jl")
include("querying.jl")
include("pinning.jl")
include("threadinfo.jl")
include("Core2CoreLatency/Core2CoreLatency.jl")
using .Core2CoreLatency
include("latency.jl")
export getcpuid, getcpuids, pinthread, pinthreads, threadinfo, @tspawnat
export sysinfo,
    nsockets,
    nnuma,
    hyperthreading_is_enabled,
    ishyperthread,
    cpuids_per_socket,
    cpuids_per_numa,
    cpuids_all
export threadinfo

end


text:
#=
utilities for benchmark scripts
=#

function setup_benchmark_dataframe()
    perf = DataFrames.DataFrame(
        example = String[],
        inst_set = String[],
        inst_num = Int[],
        inst_data = Tuple[],
        paj_ext = Bool[],
        paj_iter = Bool[],
        status = String[],
        solve_time = Float64[],
        nodes = Int[],
        iters = Int[],
        num_subp = Int[],
        num_cuts = Int[],
        primal_obj = Float64[],
        rel_gap = Float64[],
        setup_time = Float64[],
        optimize_time = Float64[],
        total_time = Float64[],
    )
    return perf
end

function write_perf(
    perf::DataFrames.DataFrame,
    results_path::Union{String, Nothing},
    new_perf::NamedTuple,
)
    push!(perf, new_perf, cols = :subset)
    if !isnothing(results_path)
        CSV.write(
            results_path,
            perf[end:end, :],
            transform = (col, val) -> something(val, missing),
            append = true,
        )
    end
    return
end

const limit_statuses = ["TIME_LIMIT", "ITERATION_LIMIT"]

function run_instance_set(
    inst_subset::Vector{<:Tuple},
    ex_type::Type{<:ExampleInstance},
    info_perf::NamedTuple,
    default_options::NamedTuple,
    perf::DataFrames.DataFrame,
    results_path::Union{String, Nothing},
    skip_limit::Bool,
)
    @testset "inst $inst_num: $(inst[1])" for (inst_num, inst) in enumerate(inst_subset)
        println("inst $inst_num: $(inst[1]) ...")

        total_time = @elapsed run_perf =
            run_instance(ex_type, inst..., default_options = default_options)
        @printf("%8.2e seconds\n", total_time)

        p = (; info_perf..., run_perf..., total_time, inst_num, :inst_data => inst[1])
        write_perf(perf, results_path, p)

        if skip_limit && run_perf.status in limit_statuses
            println("skipping remaining instances due to limit status\n")
            break
        end
        println()
    end
    return
end

function run_examples(
    examples::Vector{String},
    inst_sets::Vector{String},
    default_options::NamedTuple,
    results_path::Union{String, Nothing},
    skip_limit::Bool, # skip larger instances after hitting a solver limit
)
    # setup dataframe
    perf = setup_benchmark_dataframe()
    isnothing(results_path) || CSV.write(results_path, perf)

    # run examples
    @testset "$ex" for ex in examples
        # load instances
        (ex_type, ex_insts) = include(joinpath(@__DIR__, ex, "instances.jl"))

        # run instance sets
        for inst_set in inst_sets
            haskey(ex_insts, inst_set) || continue

            @testset "$inst_set" begin
                inst_subset = ex_insts[inst_set]
                isempty(inst_subset) && continue

                info_perf = (; inst_set, :example => ex)
                str = "$ex $inst_set"
                println("\nstarting $str tests\n")
                @testset "$str" begin
                    run_instance_set(
                        inst_subset,
                        ex_type,
                        info_perf,
                        default_options,
                        perf,
                        results_path,
                        skip_limit,
                    )
                end
            end
        end
    end

    println()
    DataFrames.show(perf, allrows = true, allcols = true)
    println("\n")
    return perf
end


text:
module ScatterNNlib

using CUDA
using DataStructures: DefaultDict
using FillArrays: Fill
using Statistics: mean
using StaticArrays: StaticArray
using Zygote
using ZygoteRules

export
    # scatter
    scatter_add!,
    scatter_sub!,
    scatter_max!,
    scatter_min!,
    scatter_mul!,
    scatter_div!,
    scatter_mean!,
    scatter!,

    # gather
    gather

const IntOrTuple = Union{Integer,Tuple}
const ops = [:add, :sub, :mul, :div, :max, :min, :mean]
const name2op = Dict(:add => :+, :sub => :-, :mul => :*, :div => :/)

include("scatter/densearray.jl")
include("scatter/staticarray.jl")
include("interfaces.jl")

include("grad.jl")

include("gather.jl")

include("utils.jl")

if CUDA.functional()
    using CUDA: @cuda
    const MAX_THREADS = 1024

    include("cuda/cuarray.jl")
    include("cuda/grad.jl")
    include("cuda/gather.jl")
else
    @warn "CUDA unavailable, not loading CUDA support"
end

end


text:
using Test
using Algames
using Altro
using BenchmarkTools
using ForwardDiff
using LinearAlgebra
using Parameters
using Printf
using Random
using RecipesBase
using RobotDynamics
using Rotations
using SparseArrays
using StaticArrays
using TrajectoryOptimization

# Dynamics
include("dynamics/double_integrator.jl")
include("dynamics/unicycle.jl")
include("dynamics/bicycle.jl")
include("dynamics/quadrotor.jl")

# Struct
include("struct/problem_size.jl")
include("struct/statistics.jl")
include("struct/primal_dual_traj.jl")
include("struct/violations.jl")
include("struct/regularizer.jl")
include("struct/options.jl")

# Core
include("core/stamp.jl")
include("core/newton_core.jl")

# Constraints
include("constraints/state_bound_constraint.jl")
include("constraints/control_bound_constraint.jl")
include("constraints/wall_constraint.jl")
include("constraints/cylinder_constraint.jl")
include("constraints/game_constraints.jl")
include("constraints/constraints_methods.jl")
include("constraints/constraint_derivatives.jl")
include("constraints/velocity_constraint.jl")

# Objective
include("objective/objective.jl")

# Problem
include("problem/problem.jl")
include("problem/local_quantities.jl")
include("problem/global_quantities.jl")
include("problem/solver_methods.jl")

# Equilibrium Subspace
include("active_set/active_set_stamp.jl")
include("active_set/active_set_core.jl")
include("active_set/active_set_methods.jl")

# Plots
include("plots/solver_plots.jl")


text:
abstract type StochasticDiffEqAlgorithm <: AbstractSDEAlgorithm end
abstract type StochasticDiffEqAdaptiveAlgorithm <: StochasticDiffEqAlgorithm end
abstract type StochasticDiffEqCompositeAlgorithm <: StochasticDiffEqAlgorithm end

abstract type StochasticDiffEqRODEAlgorithm <: AbstractRODEAlgorithm end
abstract type StochasticDiffEqRODEAdaptiveAlgorithm <: StochasticDiffEqRODEAlgorithm end
abstract type StochasticDiffEqRODECompositeAlgorithm <: StochasticDiffEqRODEAlgorithm end

abstract type StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} <: StochasticDiffEqAdaptiveAlgorithm end
abstract type StochasticDiffEqNewtonAlgorithm{CS,AD,Controller} <: StochasticDiffEqAlgorithm end

abstract type IteratedIntegralApprox end
struct IICommutative <:  IteratedIntegralApprox end
struct IIWiktorsson <:  IteratedIntegralApprox end

################################################################################

# Basics

struct EM{split} <: StochasticDiffEqAlgorithm end
EM(split=true) = EM{split}()

struct SplitEM <: StochasticDiffEqAlgorithm end
struct EulerHeun <: StochasticDiffEqAlgorithm end

struct LambaEM{split} <: StochasticDiffEqAdaptiveAlgorithm end
LambaEM(split=true) = LambaEM{split}()

struct LambaEulerHeun <: StochasticDiffEqAdaptiveAlgorithm end

struct SimplifiedEM <: StochasticDiffEqAlgorithm end

struct RKMil{interpretation} <: StochasticDiffEqAdaptiveAlgorithm end
RKMil(;interpretation=:Ito) = RKMil{interpretation}()

struct RKMilCommute{interpretation} <: StochasticDiffEqAdaptiveAlgorithm end
RKMilCommute(;interpretation=:Ito) = RKMilCommute{interpretation}()

struct RKMil_General{T<:IteratedIntegralApprox} <: StochasticDiffEqAdaptiveAlgorithm
  interpretation::Symbol
  ii_approx::T
  c::Int
end
function RKMil_General(;interpretation=:Ito, ii_approx=IIWiktorsson(), c = 1)
  RKMil_General(interpretation, ii_approx, c)
end

struct WangLi3SMil_A <: StochasticDiffEqAlgorithm end
struct WangLi3SMil_B <: StochasticDiffEqAlgorithm end
struct WangLi3SMil_C <: StochasticDiffEqAlgorithm end
struct WangLi3SMil_D <: StochasticDiffEqAlgorithm end
struct WangLi3SMil_E <: StochasticDiffEqAlgorithm end
struct WangLi3SMil_F <: StochasticDiffEqAlgorithm end

#SROCK methods
struct SROCK1{interpretation,E} <: StochasticDiffEqAlgorithm
  eigen_est::E
end
SROCK1(;interpretation=:Ito,eigen_est=nothing) = SROCK1{interpretation,typeof(eigen_est)}(eigen_est)

# Weak Order 2
for Alg in [:SROCK2, :KomBurSROCK2, :SROCKC2]
  @eval begin
    struct $Alg{E} <: StochasticDiffEqAlgorithm
      eigen_est::E
    end
    $Alg(;eigen_est=nothing) = $Alg(eigen_est)
  end
end

# ROCK stabilization for EM
struct SROCKEM{E} <: StochasticDiffEqAlgorithm
  strong_order_1::Bool
  eigen_est::E
end
SROCKEM(;strong_order_1=true,eigen_est=nothing) = SROCKEM(strong_order_1,eigen_est)

struct SKSROCK{E} <: StochasticDiffEqAlgorithm
  post_processing::Bool
  eigen_est::E
end
SKSROCK(;post_processing=false,eigen_est=nothing) = SKSROCK(post_processing,eigen_est)

struct TangXiaoSROCK2{E} <: StochasticDiffEqAlgorithm
  version_num::Int
  eigen_est::E
end
TangXiaoSROCK2(;version_num=5,eigen_est=nothing) = TangXiaoSROCK2(version_num,eigen_est)
###############################################################################

# Predictor Corrector
struct PCEuler{T<:Real, F} <: StochasticDiffEqAlgorithm
  theta::T
  eta::T
  ggprime::F
end

"""
    PCEuler(ggprime; theta=1/2, eta=1/2)

Predictor Corrector Euler

# Arguments
- `ggprime::Function`:
  For scalar problems, `ggprime` ``= b\\partial_x(b)``
  For multi-dimensional problems
  `bbprime_k` ``= \\sum_{j=1...M, i=1...D} b^(j)_i \\partial_i b^(j)_k``
  where ``b^(j)`` correspond to the noise vector due to the j'th noise channel.
  If problem is in place - a in place ggprime should be supplied - and
  vice versa for not in place speicification of problem.
- `theta::Real`:
  Degree of implicitness in the drift term. Set to 0.5 by default.
- `eta::Real`:
  Degree of implicitness in the diffusion term. Set to 0.5 by default.

Reference: Stochastics and Dynamics, Vol. 8, No. 3 (2008) 561–581
Note that the original paper has a typo in the definition of ggprime...
"""
PCEuler(ggprime; theta=1/2, eta=1/2) = PCEuler(theta,eta,ggprime)

################################################################################

# Rossler

struct SRA{TabType} <: StochasticDiffEqAdaptiveAlgorithm
  tableau::TabType
end
SRA(;tableau=constructSRA1()) = SRA(tableau)

struct SRI{TabType} <: StochasticDiffEqAdaptiveAlgorithm
  tableau::TabType
  error_terms::Int
end
SRI(;tableau=constructSRIW1(),error_terms=4) = SRI(tableau,error_terms)

struct SRIW1 <: StochasticDiffEqAdaptiveAlgorithm end
struct SRIW2 <: StochasticDiffEqAdaptiveAlgorithm end
struct SOSRI <: StochasticDiffEqAdaptiveAlgorithm end
struct SOSRI2 <: StochasticDiffEqAdaptiveAlgorithm end
struct SRA1 <: StochasticDiffEqAdaptiveAlgorithm end
struct SRA2 <: StochasticDiffEqAdaptiveAlgorithm end
struct SRA3 <: StochasticDiffEqAdaptiveAlgorithm end
struct SOSRA <: StochasticDiffEqAdaptiveAlgorithm end
struct SOSRA2 <: StochasticDiffEqAdaptiveAlgorithm end

################################################################################

# IIF

struct IIF1M{F} <: StochasticDiffEqAlgorithm
  nlsolve::F
end
IIF1M(;nlsolve=NLSOLVEJL_SETUP()) = IIF1M{typeof(nlsolve)}(nlsolve)

struct IIF2M{F} <: StochasticDiffEqAlgorithm
  nlsolve::F
end
IIF2M(;nlsolve=NLSOLVEJL_SETUP()) = IIF2M{typeof(nlsolve)}(nlsolve)

struct IIF1Mil{F} <: StochasticDiffEqAlgorithm
  nlsolve::F
end
IIF1Mil(;nlsolve=NLSOLVEJL_SETUP()) = IIF1Mil{typeof(nlsolve)}(nlsolve)

################################################################################

# SDIRK

struct ImplicitEM{CS,AD,F,F2,S,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  nlsolve::F2
  diff_type::S
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ImplicitEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                          extrapolant=:constant,
                          theta = 1/2,symplectic=false,
                          new_jac_conv_bound = 1e-3,
                          controller = :Predictive) =
                          ImplicitEM{chunk_size,autodiff,
                          typeof(linsolve),typeof(nlsolve),typeof(diff_type),
                          typeof(new_jac_conv_bound),controller}(
                          linsolve,nlsolve,diff_type,
                          symplectic ? 1/2 : theta,
                          extrapolant,new_jac_conv_bound,symplectic)

struct ImplicitEulerHeun{CS,AD,F,S,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  diff_type::S
  nlsolve::N
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ImplicitEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                          extrapolant=:constant,
                          theta = 1/2,symplectic = false,
                          new_jac_conv_bound = 1e-3,
                          controller = :Predictive) =
                          ImplicitEulerHeun{chunk_size,autodiff,
                          typeof(linsolve),typeof(diff_type),
                          typeof(nlsolve),
                          typeof(new_jac_conv_bound),controller}(
                          linsolve,diff_type,nlsolve,
                          symplectic ? 1/2 : theta,
                          extrapolant,
                          new_jac_conv_bound,symplectic)

struct ImplicitRKMil{CS,AD,F,S,N,T2,Controller,interpretation} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  diff_type::S
  nlsolve::N
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ImplicitRKMil(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                          linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                          extrapolant=:constant,
                          theta = 1/2,symplectic = false,
                          new_jac_conv_bound = 1e-3,
                          controller = :Predictive,interpretation=:Ito) =
                          ImplicitRKMil{chunk_size,autodiff,
                          typeof(linsolve),typeof(diff_type),
                          typeof(nlsolve),typeof(new_jac_conv_bound),
                          controller,interpretation}(
                          linsolve,diff_type,nlsolve,
                          symplectic ? 1/2 : theta,
                          extrapolant,
                          new_jac_conv_bound,symplectic)

struct ISSEM{CS,AD,F,S,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  diff_type::S
  nlsolve::N
  theta::T2
  extrapolant::Symbol
  new_jac_conv_bound::T2
  symplectic::Bool
end
ISSEM(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                       linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                       extrapolant=:constant,
                       theta = 1,symplectic=false,
                       new_jac_conv_bound = 1e-3,
                       controller = :Predictive) =
                       ISSEM{chunk_size,autodiff,
                       typeof(linsolve),typeof(diff_type),
                       typeof(nlsolve),
                       typeof(new_jac_conv_bound),controller}(
                       linsolve,diff_type,nlsolve,
                       symplectic ? 1/2 : theta,
                       extrapolant,
                       new_jac_conv_bound,symplectic)

struct ISSEulerHeun{CS,AD,F,S,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
 linsolve::F
 diff_type::S
 nlsolve::N
 theta::T2
 extrapolant::Symbol
 new_jac_conv_bound::T2
 symplectic::Bool
end
ISSEulerHeun(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                      linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                      extrapolant=:constant,
                      theta = 1,symplectic=false,
                      new_jac_conv_bound = 1e-3,
                      controller = :Predictive) =
                      ISSEulerHeun{chunk_size,autodiff,
                      typeof(linsolve),typeof(diff_type),
                      typeof(nlsolve),typeof(new_jac_conv_bound),controller}(
                      linsolve,diff_type,nlsolve,
                      symplectic ? 1/2 : theta,
                      extrapolant,
                      new_jac_conv_bound,symplectic)

struct SKenCarp{CS,AD,F,FDT,N,T2,Controller} <: StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}
  linsolve::F
  diff_type::FDT
  nlsolve::N
  smooth_est::Bool
  extrapolant::Symbol
  new_jac_conv_bound::T2
  ode_error_est::Bool
end

SKenCarp(;chunk_size=0,autodiff=true,diff_type=Val{:central},
                   linsolve=DEFAULT_LINSOLVE,nlsolve=NLNewton(),
                   smooth_est=true,extrapolant=:min_correct,
                   new_jac_conv_bound = 1e-3,controller = :Predictive,
                   ode_error_est = true) =
 SKenCarp{chunk_size,autodiff,typeof(linsolve),typeof(diff_type),
        typeof(nlsolve),typeof(new_jac_conv_bound),controller}(
        linsolve,diff_type,nlsolve,smooth_est,extrapolant,new_jac_conv_bound,
        ode_error_est)

################################################################################

# Etc.

struct StochasticCompositeAlgorithm{T,F} <: StochasticDiffEqCompositeAlgorithm
  algs::T
  choice_function::F
end

struct RandomEM <: StochasticDiffEqRODEAlgorithm end

const SplitSDEAlgorithms = Union{IIF1M,IIF2M,IIF1Mil,SKenCarp,SplitEM}


text:
@testset "CCall to SBML works and SBML returns a version" begin
    @test SBML.Version() isa VersionNumber
end


text:
# __precompile__()

module Fussy

  using PmapProgressMeter
  using DataStructures
  using Distributions
  using ProgressMeter
  using Parameters
  using SymEngine
  using Polynomials
  using Unitful
  using QuadGK
  using Roots
  using Optim

  include("abstracts.jl")
  include("defaults.jl")

  include("utils/index.jl")

  include("structs/index.jl")
  include("methods/index.jl")

  include("decks/index.jl")
  include("solutions/index.jl")

  export Sampling
  export Reactor
  export Study
  export Scan

  function __init__()
    global const rho_sym = symbols(:rho)
    global const eta_CD_sym = symbols(:eta_CD)

    global const symbol_scaling = Dict(
      zip(
        scaling_keys,
        map(cur_key -> symbols("alpha_$cur_key"), scaling_keys)
      )
    )
  end

end


text:
# Some simple tests to verify the sparse filtering code

using SubmodularMaximization
using HistogramFilters
using SparseArrays

g = Grid(10, 10)

thresh = 1e-3

sensor = RangingSensor()

robot_state = (2,2)
target_state = (5,5)

target_filter = SparseFilter(g, target_state, threshold=thresh)

@assert nnz(get_values(target_filter)) == 1

@assert get_states(target_filter) == [target_state]

process_update!(target_filter, transition_matrix(g))

@assert nnz(get_values(target_filter)) == 5
@assert length(get_states(g, target_filter)) == 5

range_observation = generate_observation(g, sensor, robot_state,
                                         target_state)
measurement_update!(target_filter, robot_state,
                    get_states(g, target_filter), sensor, g,
                    range_observation)

@assert nnz(get_values(target_filter)) == 5
@assert length(get_states(g, target_filter)) == 5

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
function VanderMonde(x::Array{T}, deg::Int64) where {T <: Real}
    """
        Get the vandermond matrix from a column vector, only support
        real numbers for now.
    """

    @assert ndims(x) == 2 "expect x to be a column vector of size (n, 1)"*
    string("But its dimension is: ", ndims(x))
    @assert size(x, 2) == 1 "Expect the second dimension to have size 1"*
    string(" but the actual size for x is: ", size(x))

    n = size(x, 1)
    V = zeros(n, deg) # fill with the same type.
    V[:, 1] = x
    for Col ∈ 2:deg
        V[:, Col] = x.^Col
    end
    return V
end

function PrintTitle(mesg::String)
    """

        Print a beautifully formated title to the console. stdout.
    """
    width = displaysize(stdout)[2]
    println(repeat("=" ,width))
    print("| ")
    println(mesg)
    println(repeat("=" ,width))

end

function Warn(mesg::String, newline=true)
    printstyled("$mesg $(newline ? '\n' : ' ')", color=:yellow)
end


text:
export add_gauss, add_gauss!, add_gauss_chn, add_gauss_chn!

 # function for a raw value
function f_awg(σ, μ, σ_c, μ_c)
    f(x::Real) = randn() * σ + μ

    function f(x::Complex{T}) where T
        return ((randn() * real(σ_c) + real(μ_c)
                + randn() * 1im * imag(σ_c) + 1im * imag(μ_c)))
    end
    return f
end

@inline function comb_add_gauss(x, n)
    return x + n
end

add_gauss(X, σ=0.1, μ=0; clip=false) = add_gauss!(copy(X), σ, μ, clip=clip)
add_gauss!(X, σ=0.1, μ=0; clip=false) = apply_noise!(comb_add_gauss, f_awg(σ, μ, complex_copy(σ), complex_copy(μ)), X, clip)
"""
    add_gauss(X; clip=false[, σ=0.1, μ=0.0])

Returns the array `X` with gauss noise (standard deviation `σ` and mean `μ`)
added.
`σ` and `μ` are optional arguments representing standard deviation and mean of gauss.
If keyword argument `clip` is provided the values are clipped to be in [0, 1].
If `X` is a RGB{Normed} or Gray{Normed} image, then the values will be automatically clipped and the keyword
`clip` is meaningless.

If `X<:Complex`, `μ` and `σ` are applied to the imaginary in the same way as for the real part.
If you want to have different behaviour for real and imaginary part, simply
choose `μ` or `σ` complex.
"""
add_gauss

add_gauss_chn(X, σ=0.1, μ=0; clip=false) = add_gauss_chn!(copy(X), σ, μ, clip=clip)
add_gauss_chn!(X, σ=0.1, μ=0; clip=false) = apply_noise_chn!(comb_add_gauss, f_awg(σ, μ, complex_copy(σ), complex_copy(μ)), X, clip)

"""
    add_gauss_chn(X; clip=false[, σ=0.1, μ=0.0])

Returns the RGB image `X` with gauss noise (standard deviation `σ` and mean `μ`)
added pixelwise. However, every channel of one pixel receives the same amount of noise.
The noise therefore acts roughly as intensity - but not color - changing noise.
If keyword argument `clip` is provided the values are clipped to be in [0, 1].
`σ` and `μ` are optional arguments representing standard deviation and mean of gauss.

If `X<:Complex`, `μ` and `σ` are applied to the imaginary in the same way as for the real part.
If you want to have different behaviour for real and imaginary part, simply
choose `μ` or `σ` complex.
"""
add_gauss_chn

max_issues_repo_name: roflmaostc/Noise

max_forks_repo_name: roflmaostc/Noise


text:
"""
In both digital filter design and spectral estimation, the choice of a windowing
function can play an important role in determining the quality of overall results.
The main role of the window is to damp out the effects of the Gibbs phenomenon that
results from truncation of an infinite series.

# Available Window Functions

    | Function       | Window Type                                    |
    | -------------- | ---------------------------------------------- |
    | barthann       | Modified Bartlett-Hann                         |
    | bartlett       | Bartlett                                       |
    | blackman       | Blackman                                       |
    | blackmanharris | Minimum four-term Blackman-Harris              |
    | bohman         | Bohman                                         |
    | flattop        | Flat top weighted                              |
    | hamming        | Hamming                                        |
    | hanning        | Hann (Hanning)                                 |
    | nuttall        | Nuttall-defined minimum 4-term Blackman-Harris |
    | parzen         | Parzen (de la Vallée Poussin)                  |
    | rectangular    | Rectangular                                    |
    | triangular     | Triangular                                     |
    | tukey          | Tukey (tapered cosine)                         |

"""
module WindowFunctions

export barthann
export bartlett
export blackman
export blackmanharris
export bohman
export flattop
export hamming
export hanning
export nuttall
export parzen
export rectangular
export triangular

"""
    barthann(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 35.9dB
barthann window has a mainlobe at the origin and asymptotically decaying
sidelobes on both sides. It is a linear combination of weighted Bartlett
and Hanning windows with near sidelobes lower than both Bartlett and Hanning
and with far sidelobes lower than both Bartlett and Hamming windows. The
mainlobe width of the modified Bartlett-Hann window is not increased relative
to either Bartlett or Hann window mainlobes.
"""
function barthann(N::Int)
    n = (0:N-1)/(N-1)
    return @. 0.62 - 0.48*abs(n-0.5) + 0.38*cos(2*pi*n-pi)
end

"""
    bartlett(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 26.5dB
Bartlett window is very similar to a triangular window as returned by the triang
function. However, the Bartlett window always has zeros at the first and last samples,
while the triangular window is nonzero at those points.
"""
function bartlett(L::Int)
    L==1 && return [1.0]
    N = L - 1
    w = zeros(L)
    for i = 1:L
        n = i - 1
        w[i] = 0<=n<=N/2 ? 2n/N : 2-2n/N
    end
    return w
end

"""
    blackman(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 58.1dB
It is recommended to use blackman window to detect two
signals with similar frequency but different amplitude.
"""
function blackman(N::Int)
    n = 2 * pi .* (0:(N-1))/(N-1)
	return @. 0.42 - 0.5 * cos(n) + 0.08 * cos(2n)
end

"""
    blackmanharris(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 92dB
Minimum four-term Blackman-Harris window.
"""
function blackmanharris(N::Int)
    n = 2 * pi .* (0:N-1)/(N-1)
    a₀ = 0.35875
    a₁ = 0.48829
    a₂ = 0.14128
    a₃ = 0.01168
    return @. a₀ - a₁*cos(n) + a₂*cos(2n) - a₃*cos(3n)
end

"""
    bohman(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 46dB
Bohman window is the convolution of two half-duration cosine lobes.
In the time domain, it is the product of a triangular window and a
single cycle of a cosine with a term added to set the first derivative
to zero at the boundary. Bohman windows fall off as 1/ω⁴
"""
function bohman(N::Int)
    n = (0:N-1)/(N-1)
    n = abs.((2n) .- 1)
    return @. (1-n)*cos(pi*n) + 1/pi*sin(pi*n)
end

"""
    flattop(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 93.6dB
flat-top window is summations of cosines with very small passband fluctuations.
"""
function flattop(N::Int)
    n = 2 * pi .* (0:N-1)/(N-1)
    a₀ = 0.215578950
    a₁ = 0.416631580
    a₂ = 0.277263158
    a₃ = 0.083578947
    a₄ = 0.006947368
    return @. a₀ - a₁*cos(n) + a₂*cos(2n) - a₃*cos(3n) + a₄*cos(4n)
end

"""
    hamming(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 43.2dB
Similar to hanning but has lower Sidelobe.
"""
function hamming(N::Int)
	return @. 0.54 - 0.46 * cos( 2*pi * (0:(N-1))/(N-1) )
end

"""
    hanning(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 31.5dB

If the signal is random or unknown, or has multiple frequency components, and
the test focuses on frequency points rather than energy, it is recommended to
select hanning window.
"""
function hanning(N::Int)
	return @. 0.5 - 0.5 * cos( 2*pi * (0:(N-1))/(N-1) )
end

"""
    nuttall(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 93.8dB
The window is minimum in the sense that its maximum sidelobes are minimized.
The coefficients for this window differ from the Blackman-Harris window coefficients
computed with blackmanharris and produce slightly lower sidelobes.
"""
function nuttall(N::Int)
    n = 2 * pi .* (0:N-1)/(N-1)
    a₀ = 0.3635819
    a₁ = 0.4891775
    a₂ = 0.1365995
    a₃ = 0.0106411
    return @. a₀ - a₁*cos(n) + a₂*cos(2n) - a₃*cos(3n)
end

"""
    parzen(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 53.1dB
Parzen windows are piecewise-cubic approximations of Gaussian windows.
Parzen window sidelobes fall off as 1/ω⁴
"""
function parzen(N::Int)
    w = zeros(N)
    for i = 1:N
        m = i - 1
        n = abs(m - (N-1)/2)
        w[i] = n<(N-1)/4 ? 1 - 6(2n/N)^2 + 6(2n/N)^3 : 2(1-2n/N)^3
    end
    return w
end

"""
    rectangular(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 13.3dB
"""
function rectangular(N::Int)
    return ones(N)
end

"""
    triangular(N::Int) -> AbstractArray
    Sidelobe peak attenuation = 26.5dB
Compared with a rectangular window, the mainlobe is about twice wide, but the
sidelobes are smaller and have no negative sidelobes.
"""
function triangular(N::Int)
    N%2==0 && return triangulareven(N)
    N%2==1 && return triangularodd(N)
end

function triangularodd(N::Int)
    w = zeros(N)
    for n = 1:N
        w[n] = 1<=n<=(N+1)/2 ? 2n/(N+1) : 2-2n/(N+1)
    end
    return w
end

function triangulareven(N::Int)
    w = zeros(N)
    for n = 1:N
        w[n] = 1<=n<=N/2 ? (2n-1)/N : 2-(2n-1)/N
    end
    return w
end

end # WindowFunctions


text:
function checkSameLength(u1::FloatVec, u2::FloatVec)
    if !(size(u1, 1) == size(u2, 1))
        error("u1 and u2 vector must be of same length")
    end
    return true
end

@doc doc"""
get family and parameter of copula instance
"""->
function getFamAndParams(cop::ParamPC_Cpp)
    return (float(getVineCppId(cop)), params(cop))
end

## this functions exists so that copula types at some point can be
## changed to follow the interface of Distributions.jl more closely.
## This means that each copula type will have multiple fields, one for
## each parameter, and function params(cop::ParamCop) will generally
## return a tuple. This allows additional flexibility, since a single
## parameter then can be more than just a single Float value, but it
## could also be a covariance matrix or a PCModification.
function params(cop::ParamPC_Cpp)
    return cop.params
end

## triangular indexing
##--------------------

@doc doc"""
transform single index of array to double index of upper triangular
matrix
""" ->
function triang2sub(nVars::Int64, ind::Int64)
    # indices are running column-wise
    elemsInRow = [(nVars-1):-1:1]
    lastInColumn = cumsum(elemsInRow)
    rowInd = findfirst(ind .<= lastInColumn)
    if rowInd == 1
        colInd = ind + 1
    else
        nInRow = ind - lastInColumn[rowInd-1]
        colInd = rowInd + nInRow
    end
    return (rowInd, colInd)
end

function triang2sub(nVars::Int64, inds::Array{Int, 1})
    nInds = length(inds)
    rowInds = zeros(Int, nInds)
    colInds = zeros(Int, nInds)
    for ii=1:nInds
        rowInd, colInd = triang2sub(nVars, inds[ii])
        rowInds[ii] = rowInd
        colInds[ii] = colInd
    end
    return rowInds, colInds
end

function sub2triang(nVars::Int64, rowInd::Int64, colInd::Int64)
    if !(rowInd < colInd <= nVars)
        error("column index must be larger than row index")
    end
    elemsInRow = [(nVars-1):-1:1]
    fullRows = sum(elemsInRow[1:(rowInd - 1)])
    return fullRows + (colInd - rowInd)
end


text:
@testset "empty" begin
    canvas = BrailleCanvas(40, 10, origin_x = 0., origin_y = 0., width = 1., height = 1.)
    p = @inferred Plot(canvas)
    @test_throws ArgumentError Plot(canvas, margin=-1)
    test_ref("references/plot/empty.txt", @show_col(p))
    canvas = BrailleCanvas(30, 5, origin_x = 0., origin_y = 0., width = 1., height = 1.)
    p = @inferred Plot(canvas)
    test_ref("references/plot/empty_small.txt", @show_col(p))
end

canvas = BrailleCanvas(40, 10, origin_x = 0., origin_y = 0., width = 1., height = 1.)
lines!(canvas, 0., 0., 1., 1., :blue)
lines!(canvas, .2, .7, 1., 0., :red)
lines!(canvas, 0., 2., .5, 0., :green)
points!(canvas, .5, .9)

@testset "canvas / xylabel" begin
    p = @inferred Plot(canvas)
    test_ref("references/plot/canvas_only.txt", @print_col(p))
    test_ref("references/plot/canvas_only.txt", @show_col(p))

    for border in (:solid, :corners, :barplot, :bold, :ascii, :none, :dashed, :dotted)
        local p = @inferred Plot(canvas, border = border)
        test_ref("references/plot/border_$(string(border)).txt", @show_col(p))
    end

    p = @inferred Plot(canvas, xlabel = "x", ylabel= "y")
    test_ref("references/plot/xylabel.txt", @show_col(p))

    p = @inferred Plot(canvas, xlabel = "x", ylabel= "y", compact = true)
    test_ref("references/plot/xylabel_compact.txt", @show_col(p))
end

@testset "padding" begin
    p = @inferred Plot(canvas, title = "testtitle", xlabel = "x", ylabel= "y", margin = 4, padding = 5)
    @test @inferred(label!(p, :l, ":l auto 1")) === p
    @test @inferred(label!(p, :r, ":r auto 1", :red)) === p
    test_ref("references/plot/padding.txt", @show_col(p))

    p = @inferred Plot(
        canvas, title = "testtitle", xlabel = "x", ylabel= "y", margin = 4, padding = 5, labels = false
    )
    @test @inferred(label!(p, :l, ":l auto 1")) === p
    @test @inferred(label!(p, :r, ":r auto 1", :red)) === p
    test_ref("references/plot/padding_nolabels.txt", @show_col(p))
end

@testset "title / labels / full deco" begin
    p = @inferred Plot(canvas, title = "testtitle")
    test_ref("references/plot/title.txt", @show_col(p))
    @test @inferred(label!(p, :l, ":l auto 1")) === p
    @test @inferred(label!(p, :r, ":r auto 1", :red)) === p
    @test @inferred(label!(p, :l, ":l auto 2", :green)) === p
    @test @inferred(label!(p, :r, ":r auto 2", :blue)) === p
    test_ref("references/plot/title_auto.txt", @show_col(p))

    for i in 1:10
        @test @inferred(label!(p, :l, i, "$i", color=:green)) === p
        @test @inferred(label!(p, :r, i, "$i", :yellow)) === p
    end
    @test_throws ArgumentError label!(p, :bl, 5, ":l  5", :red)
    @test_throws ArgumentError label!(p, :d, ":l auto", :red)
    @test @inferred(label!(p, :l, 5, ":l  5", :red)) === p
    @test @inferred(label!(p, :r, 5, "5  :r", color=:blue)) === p
    @test @inferred(label!(p, :t, ":t", color=:yellow)) === p
    @test @inferred(label!(p, :tl, ":tl", :cyan)) === p
    @test @inferred(label!(p, :tr, ":tr")) === p
    label!(p, :bl, ":bl", :blue)
    label!(p, :b, ":b", :green)
    label!(p, :br, ":br", :white)

    @test @inferred(lines!(p, 0., 1., 1., 0., color=(0,0,255))) === p
    @test @inferred(pixel!(p, 10, 1, color = (255,255,0))) === p
    @test @inferred(points!(p, 0.05, .75, color = (0,255,0))) === p

    ttl = "title!(plot, text)"
    @test @inferred(title!(p, ttl)) === p
    @test @inferred(title(p)) == ttl
    xlab = "xlabel!(plot, text)"
    @test @inferred(xlabel!(p, xlab)) === p
    @test @inferred(xlabel(p)) == xlab
    ylab = "ylabel!(plot, text)"
    @test @inferred(ylabel!(p, ylab)) === p
    @test @inferred(ylabel(p)) == ylab
    zlab = "zlabel!(plot, text)"
    @test @inferred(zlabel!(p, zlab)) === p
    @test @inferred(zlabel(p)) == zlab

    test_ref("references/plot/full_deco.txt", @show_col(p))
    tmp = tempname()
    savefig(p, tmp; color=true)
    test_ref("references/plot/full_deco.txt", read(tmp, String))

    test_ref("references/plot/full_deco_nocolor.txt", @show_nocol(p))
    tmp = tempname()
    savefig(p, tmp; color=false)
    test_ref("references/plot/full_deco_nocolor.txt", read(tmp, String))
end

@testset "annotations" begin
    p = @inferred Plot(canvas)
    @test @inferred(annotate!(p, 0.5, 0.5, "Origin")) === p
    @test @inferred(annotate!(p, 0.5, 1.0, "North")) === p
    @test @inferred(annotate!(p, 1.0, 1.0, "North East"; halign = :right)) === p
    @test @inferred(annotate!(p, 1.0, 0.5, "East"; halign = :right)) === p
    @test @inferred(annotate!(p, 1.0, 0.0, "South East"; halign = :right)) === p
    @test @inferred(annotate!(p, 0.5, 0.0, "South")) === p
    @test @inferred(annotate!(p, 0.0, 0.0, "South West"; halign = :left)) === p
    @test @inferred(annotate!(p, 0.0, 0.5, "West"; halign = :left)) === p
    @test @inferred(annotate!(p, 0.0, 1.0, "North West"; halign = :left)) === p
    test_ref("references/plot/annotations_BrailleCanvas.txt", @show_col(p))

    @test @inferred(annotate!(p, 0.0, 0.0, 'c')) === p  # Char

    for sym in (:AsciiCanvas, :DotCanvas, :BlockCanvas)
        p = lineplot([-1, 1], [-1, 1], canvas = getproperty(UnicodePlots, sym))
        annotate!(p, +0, +0, "Origin")
        annotate!(p, +0, +1, "North")
        annotate!(p, +1, +1, "North East"; halign = :right)
        annotate!(p, +1, +0, "East"; halign = :right)
        annotate!(p, +1, -1, "South East"; halign = :right)
        annotate!(p, +0, -1, "South")
        annotate!(p, -1, -1, "South West"; halign = :left)
        annotate!(p, -1, +0, "West"; halign = :left)
        annotate!(p, -1, +1, "North West"; halign = :left)
        test_ref("references/plot/annotations_$sym.txt", @show_col(p))
    end

    p = lineplot([-1, 1], [-1, 1])
    for h in (:right, :left, :center, :hcenter), v in (:top, :bottom, :center, :vcenter)
        annotate!(p, 0, 0, "Origin"; halign = h, valign = v)
    end
end


text:
# math functionality

## trigonometric

@inline cos(x::Float64) = ccall("extern __nv_cos", llvmcall, Cdouble, (Cdouble,), x)
@inline cos(x::Float32) = ccall("extern __nv_cosf", llvmcall, Cfloat, (Cfloat,), x)
@inline cos_fast(x::Float32) = ccall("extern __nv_fast_cosf", llvmcall, Cfloat, (Cfloat,), x)

@inline cospi(x::Float64) = ccall("extern __nv_cospi", llvmcall, Cdouble, (Cdouble,), x)
@inline cospi(x::Float32) = ccall("extern __nv_cospif", llvmcall, Cfloat, (Cfloat,), x)

@inline sin(x::Float64) = ccall("extern __nv_sin", llvmcall, Cdouble, (Cdouble,), x)
@inline sin(x::Float32) = ccall("extern __nv_sinf", llvmcall, Cfloat, (Cfloat,), x)
@inline sin_fast(x::Float32) = ccall("extern __nv_fast_sinf", llvmcall, Cfloat, (Cfloat,), x)

@inline sinpi(x::Float64) = ccall("extern __nv_sinpi", llvmcall, Cdouble, (Cdouble,), x)
@inline sinpi(x::Float32) = ccall("extern __nv_sinpif", llvmcall, Cfloat, (Cfloat,), x)

@inline tan(x::Float64) = ccall("extern __nv_tan", llvmcall, Cdouble, (Cdouble,), x)
@inline tan(x::Float32) = ccall("extern __nv_tanf", llvmcall, Cfloat, (Cfloat,), x)
@inline tan_fast(x::Float32) = ccall("extern __nv_fast_tanf", llvmcall, Cfloat, (Cfloat,), x)

## inverse trigonometric

@inline acos(x::Float64) = ccall("extern __nv_acos", llvmcall, Cdouble, (Cdouble,), x)
@inline acos(x::Float32) = ccall("extern __nv_acosf", llvmcall, Cfloat, (Cfloat,), x)

@inline asin(x::Float64) = ccall("extern __nv_asin", llvmcall, Cdouble, (Cdouble,), x)
@inline asin(x::Float32) = ccall("extern __nv_asinf", llvmcall, Cfloat, (Cfloat,), x)

@inline atan(x::Float64) = ccall("extern __nv_atan", llvmcall, Cdouble, (Cdouble,), x)
@inline atan(x::Float32) = ccall("extern __nv_atanf", llvmcall, Cfloat, (Cfloat,), x)

# ! CUDAnative.atan2 is equivalent to Base.atan
@inline atan2(x::Float64, y::Float64) = ccall("extern __nv_atan2", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline atan2(x::Float32, y::Float32) = ccall("extern __nv_atan2f", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline angle(x::ComplexF64) = atan2(x.im, x.re)
@inline angle(x::ComplexF32) = atan2(x.im, x.re)
@inline angle(x::Float64) = signbit(x,) * 3.141592653589793
@inline angle(x::Float32) = signbit(x,) * 3.1415927f0

## hyperbolic

@inline cosh(x::Float64) = ccall("extern __nv_cosh", llvmcall, Cdouble, (Cdouble,), x)
@inline cosh(x::Float32) = ccall("extern __nv_coshf", llvmcall, Cfloat, (Cfloat,), x)

@inline sinh(x::Float64) = ccall("extern __nv_sinh", llvmcall, Cdouble, (Cdouble,), x)
@inline sinh(x::Float32) = ccall("extern __nv_sinhf", llvmcall, Cfloat, (Cfloat,), x)

@inline tanh(x::Float64) = ccall("extern __nv_tanh", llvmcall, Cdouble, (Cdouble,), x)
@inline tanh(x::Float32) = ccall("extern __nv_tanhf", llvmcall, Cfloat, (Cfloat,), x)

## inverse hyperbolic

@inline acosh(x::Float64) = ccall("extern __nv_acosh", llvmcall, Cdouble, (Cdouble,), x)
@inline acosh(x::Float32) = ccall("extern __nv_acoshf", llvmcall, Cfloat, (Cfloat,), x)

@inline asinh(x::Float64) = ccall("extern __nv_asinh", llvmcall, Cdouble, (Cdouble,), x)
@inline asinh(x::Float32) = ccall("extern __nv_asinhf", llvmcall, Cfloat, (Cfloat,), x)

@inline atanh(x::Float64) = ccall("extern __nv_atanh", llvmcall, Cdouble, (Cdouble,), x)
@inline atanh(x::Float32) = ccall("extern __nv_atanhf", llvmcall, Cfloat, (Cfloat,), x)

## logarithmic

@inline log(x::Float64) = ccall("extern __nv_log", llvmcall, Cdouble, (Cdouble,), x)
@inline log(x::Float32) = ccall("extern __nv_logf", llvmcall, Cfloat, (Cfloat,), x)
@inline log_fast(x::Float32) = ccall("extern __nv_fast_logf", llvmcall, Cfloat, (Cfloat,), x)

@inline log(x::ComplexF64) = log(abs(x,)) + im * angle(x,)
@inline log(x::ComplexF32) = log(abs(x,)) + im * angle(x,)
@inline log_fast(x::ComplexF32) = log_fast(abs(x,)) + im * angle(x,)

@inline log10(x::Float64) = ccall("extern __nv_log10", llvmcall, Cdouble, (Cdouble,), x)
@inline log10(x::Float32) = ccall("extern __nv_log10f", llvmcall, Cfloat, (Cfloat,), x)
@inline log10_fast(x::Float32) = ccall("extern __nv_fast_log10f", llvmcall, Cfloat, (Cfloat,), x)

@inline log1p(x::Float64) = ccall("extern __nv_log1p", llvmcall, Cdouble, (Cdouble,), x)
@inline log1p(x::Float32) = ccall("extern __nv_log1pf", llvmcall, Cfloat, (Cfloat,), x)

@inline log2(x::Float64) = ccall("extern __nv_log2", llvmcall, Cdouble, (Cdouble,), x)
@inline log2(x::Float32) = ccall("extern __nv_log2f", llvmcall, Cfloat, (Cfloat,), x)
@inline log2_fast(x::Float32) = ccall("extern __nv_fast_log2f", llvmcall, Cfloat, (Cfloat,), x)

@inline logb(x::Float64) = ccall("extern __nv_logb", llvmcall, Cdouble, (Cdouble,), x)
@inline logb(x::Float32) = ccall("extern __nv_logbf", llvmcall, Cfloat, (Cfloat,), x)

@inline ilogb(x::Float64) = ccall("extern __nv_ilogb", llvmcall, Int32, (Cdouble,), x)
@inline ilogb(x::Float32) = ccall("extern __nv_ilogbf", llvmcall, Int32, (Cfloat,), x)

## exponential

@inline exp(x::Float64) = ccall("extern __nv_exp", llvmcall, Cdouble, (Cdouble,), x)
@inline exp(x::Float32) = ccall("extern __nv_expf", llvmcall, Cfloat, (Cfloat,), x)
@inline exp_fast(x::Float32) = ccall("extern __nv_fast_expf", llvmcall, Cfloat, (Cfloat,), x)

@inline exp2(x::Float64) = ccall("extern __nv_exp2", llvmcall, Cdouble, (Cdouble,), x)
@inline exp2(x::Float32) = ccall("extern __nv_exp2f", llvmcall, Cfloat, (Cfloat,), x)

@inline exp10(x::Float64) = ccall("extern __nv_exp10", llvmcall, Cdouble, (Cdouble,), x)
@inline exp10(x::Float32) = ccall("extern __nv_exp10f", llvmcall, Cfloat, (Cfloat,), x)
@inline exp10_fast(x::Float32) = ccall("extern __nv_fast_exp10f", llvmcall, Cfloat, (Cfloat,), x)

@inline expm1(x::Float64) = ccall("extern __nv_expm1", llvmcall, Cdouble, (Cdouble,), x)
@inline expm1(x::Float32) = ccall("extern __nv_expm1f", llvmcall, Cfloat, (Cfloat,), x)

@inline ldexp(x::Float64, y::Int32) = ccall("extern __nv_ldexp", llvmcall, Cdouble, (Cdouble, Int32), x, y)
@inline ldexp(x::Float32, y::Int32) = ccall("extern __nv_ldexpf", llvmcall, Cfloat, (Cfloat, Int32), x, y)

@inline exp(x::Complex{Float64}) = exp(x.re) * (cos(x.im) + 1.0im * sin(x.im))
@inline exp(x::Complex{Float32}) = exp(x.re) * (cos(x.im) + 1.0f0im * sin(x.im))
@inline exp_fast(x::Complex{Float32}) = exp_fast(x.re) * (cos_fast(x.im) + 1.0f0im * sin_fast(x.im))

## error

@inline erf(x::Float64) = ccall("extern __nv_erf", llvmcall, Cdouble, (Cdouble,), x)
@inline erf(x::Float32) = ccall("extern __nv_erff", llvmcall, Cfloat, (Cfloat,), x)

@inline erfinv(x::Float64) = ccall("extern __nv_erfinv", llvmcall, Cdouble, (Cdouble,), x)
@inline erfinv(x::Float32) = ccall("extern __nv_erfinvf", llvmcall, Cfloat, (Cfloat,), x)

@inline erfc(x::Float64) = ccall("extern __nv_erfc", llvmcall, Cdouble, (Cdouble,), x)
@inline erfc(x::Float32) = ccall("extern __nv_erfcf", llvmcall, Cfloat, (Cfloat,), x)

@inline erfcinv(x::Float64) = ccall("extern __nv_erfcinv", llvmcall, Cdouble, (Cdouble,), x)
@inline erfcinv(x::Float32) = ccall("extern __nv_erfcinvf", llvmcall, Cfloat, (Cfloat,), x)

@inline erfcx(x::Float64) = ccall("extern __nv_erfcx", llvmcall, Cdouble, (Cdouble,), x)
@inline erfcx(x::Float32) = ccall("extern __nv_erfcxf", llvmcall, Cfloat, (Cfloat,), x)

## integer handling (bit twiddling)

@inline brev(x::Int32) =   ccall("extern __nv_brev", llvmcall, Int32, (Int32,), x)
@inline brev(x::Int64) =   ccall("extern __nv_brevll", llvmcall, Int64, (Int64,), x)

@inline clz(x::Int32) =   ccall("extern __nv_clz", llvmcall, Int32, (Int32,), x)
@inline clz(x::Int64) =   ccall("extern __nv_clzll", llvmcall, Int32, (Int64,), x)

@inline ffs(x::Int32) = ccall("extern __nv_ffs", llvmcall, Int32, (Int32,), x)
@inline ffs(x::Int64) = ccall("extern __nv_ffsll", llvmcall, Int32, (Int64,), x)

@inline byte_perm(x::Int32, y::Int32, z::Int32) = ccall("extern __nv_byte_perm", llvmcall, Int32, (Int32, Int32, Int32), x, y, z)

@inline popc(x::Int32) = ccall("extern __nv_popc", llvmcall, Int32, (Int32,), x)
@inline popc(x::Int64) = ccall("extern __nv_popcll", llvmcall, Int32, (Int64,), x)

## floating-point handling

@inline isfinite(x::Float32) = (ccall("extern __nv_finitef", llvmcall, Int32, (Cfloat,), x)) != 0
@inline isfinite(x::Float64) = (ccall("extern __nv_isfinited", llvmcall, Int32, (Cdouble,), x)) != 0

@inline isinf(x::Float64) = (ccall("extern __nv_isinfd", llvmcall, Int32, (Cdouble,), x)) != 0
@inline isinf(x::Float32) = (ccall("extern __nv_isinff", llvmcall, Int32, (Cfloat,), x)) != 0

@inline isnan(x::Float64) = (ccall("extern __nv_isnand", llvmcall, Int32, (Cdouble,), x)) != 0
@inline isnan(x::Float32) = (ccall("extern __nv_isnanf", llvmcall, Int32, (Cfloat,), x)) != 0

@inline nearbyint(x::Float64) = ccall("extern __nv_nearbyint", llvmcall, Cdouble, (Cdouble,), x)
@inline nearbyint(x::Float32) = ccall("extern __nv_nearbyintf", llvmcall, Cfloat, (Cfloat,), x)

@inline nextafter(x::Float64, y::Float64) = ccall("extern __nv_nextafter", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline nextafter(x::Float32, y::Float32) = ccall("extern __nv_nextafterf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

## sign handling

@inline signbit(x::Float64) = (ccall("extern __nv_signbitd", llvmcall, Int32, (Cdouble,), x)) != 0
@inline signbit(x::Float32) = (ccall("extern __nv_signbitf", llvmcall, Int32, (Cfloat,), x)) != 0

@inline copysign(x::Float64, y::Float64) = ccall("extern __nv_copysign", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline copysign(x::Float32, y::Float32) = ccall("extern __nv_copysignf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline abs(x::Int32) =   ccall("extern __nv_abs", llvmcall, Int32, (Int32,), x)
@inline abs(f::Float64) = ccall("extern __nv_fabs", llvmcall, Cdouble, (Cdouble,), f)
@inline abs(f::Float32) = ccall("extern __nv_fabsf", llvmcall, Cfloat, (Cfloat,), f)
@inline abs(x::Int64) =   ccall("extern __nv_llabs", llvmcall, Int64, (Int64,), x)

@inline abs(x::Complex{Float64}) = hypot(x.re, x.im)
@inline abs(x::Complex{Float32}) = hypot(x.re, x.im)

## roots and powers

@inline sqrt(x::Float64) = ccall("extern __nv_sqrt", llvmcall, Cdouble, (Cdouble,), x)
@inline sqrt(x::Float32) = ccall("extern __nv_sqrtf", llvmcall, Cfloat, (Cfloat,), x)

@inline rsqrt(x::Float64) = ccall("extern __nv_rsqrt", llvmcall, Cdouble, (Cdouble,), x)
@inline rsqrt(x::Float32) = ccall("extern __nv_rsqrtf", llvmcall, Cfloat, (Cfloat,), x)

@inline cbrt(x::Float64) = ccall("extern __nv_cbrt", llvmcall, Cdouble, (Cdouble,), x)
@inline cbrt(x::Float32) = ccall("extern __nv_cbrtf", llvmcall, Cfloat, (Cfloat,), x)

@inline rcbrt(x::Float64) = ccall("extern __nv_rcbrt", llvmcall, Cdouble, (Cdouble,), x)
@inline rcbrt(x::Float32) = ccall("extern __nv_rcbrtf", llvmcall, Cfloat, (Cfloat,), x)

@inline pow(x::Float64, y::Float64) = ccall("extern __nv_pow", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline pow(x::Float32, y::Float32) = ccall("extern __nv_powf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)
@inline pow_fast(x::Float32, y::Float32) = ccall("extern __nv_fast_powf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)
@inline pow(x::Float64, y::Int32) =   ccall("extern __nv_powi", llvmcall, Cdouble, (Cdouble, Int32), x, y)
@inline pow(x::Float32, y::Int32) =   ccall("extern __nv_powif", llvmcall, Cfloat, (Cfloat, Int32), x, y)
@inline pow(x::Union{Float32, Float64}, y::Int64) = pow(x, Int32(y,))

@inline abs2(x::Complex{Float64}) = x.re * x.re + x.im * x.im
@inline abs2(x::Complex{Float32}) = x.re * x.re + x.im * x.im

## rounding and selection

# TODO: differentiate in return type, map correctly
# @inline round(x::Float64) = ccall("extern __nv_llround", llvmcall, Int64, (Cdouble,), x)
# @inline round(x::Float32) = ccall("extern __nv_llroundf", llvmcall, Int64, (Cfloat,), x)
# @inline round(x::Float64) = ccall("extern __nv_round", llvmcall, Cdouble, (Cdouble,), x)
# @inline round(x::Float32) = ccall("extern __nv_roundf", llvmcall, Cfloat, (Cfloat,), x)

# TODO: differentiate in return type, map correctly
# @inline rint(x::Float64) = ccall("extern __nv_llrint", llvmcall, Int64, (Cdouble,), x)
# @inline rint(x::Float32) = ccall("extern __nv_llrintf", llvmcall, Int64, (Cfloat,), x)
# @inline rint(x::Float64) = ccall("extern __nv_rint", llvmcall, Cdouble, (Cdouble,), x)
# @inline rint(x::Float32) = ccall("extern __nv_rintf", llvmcall, Cfloat, (Cfloat,), x)

# TODO: would conflict with trunc usage in this module
# @inline trunc(x::Float64) = ccall("extern __nv_trunc", llvmcall, Cdouble, (Cdouble,), x)
# @inline trunc(x::Float32) = ccall("extern __nv_truncf", llvmcall, Cfloat, (Cfloat,), x)

@inline ceil(x::Float64) = ccall("extern __nv_ceil", llvmcall, Cdouble, (Cdouble,), x)
@inline ceil(x::Float32) = ccall("extern __nv_ceilf", llvmcall, Cfloat, (Cfloat,), x)

@inline floor(f::Float64) = ccall("extern __nv_floor", llvmcall, Cdouble, (Cdouble,), f)
@inline floor(f::Float32) = ccall("extern __nv_floorf", llvmcall, Cfloat, (Cfloat,), f)

@inline min(x::Int32, y::Int32) = ccall("extern __nv_min", llvmcall, Int32, (Int32, Int32), x, y)
@inline min(x::Int64, y::Int64) = ccall("extern __nv_llmin", llvmcall, Int64, (Int64, Int64), x, y)
@inline min(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umin", llvmcall, Int32, (Int32, Int32), x, y))
@inline min(x::UInt64, y::UInt64) = convert(UInt64, ccall("extern __nv_ullmin", llvmcall, Int64, (Int64, Int64), x, y))
@inline min(x::Float64, y::Float64) = ccall("extern __nv_fmin", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline min(x::Float32, y::Float32) = ccall("extern __nv_fminf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline max(x::Int32, y::Int32) = ccall("extern __nv_max", llvmcall, Int32, (Int32, Int32), x, y)
@inline max(x::Int64, y::Int64) = ccall("extern __nv_llmax", llvmcall, Int64, (Int64, Int64), x, y)
@inline max(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umax", llvmcall, Int32, (Int32, Int32), x, y))
@inline max(x::UInt64, y::UInt64) = convert(UInt64, ccall("extern __nv_ullmax", llvmcall, Int64, (Int64, Int64), x, y))
@inline max(x::Float64, y::Float64) = ccall("extern __nv_fmax", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline max(x::Float32, y::Float32) = ccall("extern __nv_fmaxf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline saturate(x::Float32) = ccall("extern __nv_saturatef", llvmcall, Cfloat, (Cfloat,), x)

## division and remainder

@inline mod(x::Float64, y::Float64) = ccall("extern __nv_fmod", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline mod(x::Float32, y::Float32) = ccall("extern __nv_fmodf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline rem(x::Float64, y::Float64) = ccall("extern __nv_remainder", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline rem(x::Float32, y::Float32) = ccall("extern __nv_remainderf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline div_fast(x::Float32, y::Float32) = ccall("extern __nv_fast_fdividef", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

## gamma function

@inline lgamma(x::Float64) = ccall("extern __nv_lgamma", llvmcall, Cdouble, (Cdouble,), x)
@inline lgamma(x::Float32) = ccall("extern __nv_lgammaf", llvmcall, Cfloat, (Cfloat,), x)

@inline tgamma(x::Float64) = ccall("extern __nv_tgamma", llvmcall, Cdouble, (Cdouble,), x)
@inline tgamma(x::Float32) = ccall("extern __nv_tgammaf", llvmcall, Cfloat, (Cfloat,), x)

## Bessel

@inline j0(x::Float64) = ccall("extern __nv_j0", llvmcall, Cdouble, (Cdouble,), x)
@inline j0(x::Float32) = ccall("extern __nv_j0f", llvmcall, Cfloat, (Cfloat,), x)

@inline j1(x::Float64) = ccall("extern __nv_j1", llvmcall, Cdouble, (Cdouble,), x)
@inline j1(x::Float32) = ccall("extern __nv_j1f", llvmcall, Cfloat, (Cfloat,), x)

@inline jn(n::Int32, x::Float64) = ccall("extern __nv_jn", llvmcall, Cdouble, (Int32, Cdouble), n, x)
@inline jn(n::Int32, x::Float32) = ccall("extern __nv_jnf", llvmcall, Cfloat, (Int32, Cfloat), n, x)

@inline y0(x::Float64) = ccall("extern __nv_y0", llvmcall, Cdouble, (Cdouble,), x)
@inline y0(x::Float32) = ccall("extern __nv_y0f", llvmcall, Cfloat, (Cfloat,), x)

@inline y1(x::Float64) = ccall("extern __nv_y1", llvmcall, Cdouble, (Cdouble,), x)
@inline y1(x::Float32) = ccall("extern __nv_y1f", llvmcall, Cfloat, (Cfloat,), x)

@inline yn(n::Int32, x::Float64) = ccall("extern __nv_yn", llvmcall, Cdouble, (Int32, Cdouble), n, x)
@inline yn(n::Int32, x::Float32) = ccall("extern __nv_ynf", llvmcall, Cfloat, (Int32, Cfloat), n, x)

## distributions

@inline normcdf(x::Float64) = ccall("extern __nv_normcdf", llvmcall, Cdouble, (Cdouble,), x)
@inline normcdf(x::Float32) = ccall("extern __nv_normcdff", llvmcall, Cfloat, (Cfloat,), x)

@inline normcdfinv(x::Float64) = ccall("extern __nv_normcdfinv", llvmcall, Cdouble, (Cdouble,), x)
@inline normcdfinv(x::Float32) = ccall("extern __nv_normcdfinvf", llvmcall, Cfloat, (Cfloat,), x)

#
# Unsorted
#

@inline hypot(x::Float64, y::Float64) = ccall("extern __nv_hypot", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline hypot(x::Float32, y::Float32) = ccall("extern __nv_hypotf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline fma(x::Float64, y::Float64, z::Float64) = ccall("extern __nv_fma", llvmcall, Cdouble, (Cdouble, Cdouble, Cdouble), x, y, z)
@inline fma(x::Float32, y::Float32, z::Float32) = ccall("extern __nv_fmaf", llvmcall, Cfloat, (Cfloat, Cfloat, Cfloat), x, y, z)

@inline sad(x::Int32, y::Int32, z::Int32) = ccall("extern __nv_sad", llvmcall, Int32, (Int32, Int32, Int32), x, y, z)
@inline sad(x::UInt32, y::UInt32, z::UInt32) = convert(UInt32, ccall("extern __nv_usad", llvmcall, Int32, (Int32, Int32, Int32), x, y, z))

@inline dim(x::Float64, y::Float64) = ccall("extern __nv_fdim", llvmcall, Cdouble, (Cdouble, Cdouble), x, y)
@inline dim(x::Float32, y::Float32) = ccall("extern __nv_fdimf", llvmcall, Cfloat, (Cfloat, Cfloat), x, y)

@inline mul24(x::Int32, y::Int32) = ccall("extern __nv_mul24", llvmcall, Int32, (Int32, Int32), x, y)
@inline mul24(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umul24", llvmcall, Int32, (Int32, Int32), x, y))

@inline mul64hi(x::Int64, y::Int64) = ccall("extern __nv_mul64hi", llvmcall, Int64, (Int64, Int64), x, y)
@inline mul64hi(x::UInt64, y::UInt64) = convert(UInt64, ccall("extern __nv_umul64hi", llvmcall, Int64, (Int64, Int64), x, y))
@inline mulhi(x::Int32, y::Int32) = ccall("extern __nv_mulhi", llvmcall, Int32, (Int32, Int32), x, y)
@inline mulhi(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_umulhi", llvmcall, Int32, (Int32, Int32), x, y))

@inline hadd(x::Int32, y::Int32) = ccall("extern __nv_hadd", llvmcall, Int32, (Int32, Int32), x, y)
@inline hadd(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_uhadd", llvmcall, Int32, (Int32, Int32), x, y))

@inline rhadd(x::Int32, y::Int32) = ccall("extern __nv_rhadd", llvmcall, Int32, (Int32, Int32), x, y)
@inline rhadd(x::UInt32, y::UInt32) = convert(UInt32, ccall("extern __nv_urhadd", llvmcall, Int32, (Int32, Int32), x, y))

@inline scalbn(x::Float64, y::Int32) = ccall("extern __nv_scalbn", llvmcall, Cdouble, (Cdouble, Int32), x, y)
@inline scalbn(x::Float32, y::Int32) = ccall("extern __nv_scalbnf", llvmcall, Cfloat, (Cfloat, Int32), x, y)


text:
function LinearAlgebra.svd(A::Array{T,2}) where T<:ForwardDiff.Dual
    m = size(A,1)
    n = size(A,2)
    @assert m >= n

    Av = map(ForwardDiff.value,A)
    Uv,sv,Vv = svd(Av)

    Ap = map(ForwardDiff.partials,A)
    UvtApVv = Uv'*Ap*Vv
    sp = diag(UvtApVv)

    # if allunique(sv)
        # If all unique, then gradient well defined
        Sv = Matrix(Diagonal(sv))
        svi = 1. ./ sv
        svi[isinf.(svi)] .= 0.
        Svi = Matrix(Diagonal(svi))

        F = zeros(eltype(Av),m,n)
        for i = 1:m
            for j = 1:n
                if i!=j
                    F[i,j] = 1. / (sv[j]^2 - sv[i]^2)
                end
            end
        end
        F[isinf.(F)] .= 0.

        VvtApUv = Vv'*Ap[reshape(1:m*n,(m,n))']*Uv

        Up = Uv*(F .* (UvtApVv*Sv + Sv*VvtApUv)) + (I - Uv*Uv')*Ap*Vv*Svi
        Vp = Vv*(F .* (Sv*UvtApVv+VvtApUv*Sv)) + (I - Vv*Vv')*Ap[reshape(1:m*n,(m,n))']*Uv*Svi
    # else
    #     # Otherwise need to least-squares
    #     dU = zeros(size(Uv,1),size(Uv,2),size(A,1),size(A,2))
    #     dV = zeros(size(Vv,1),size(Vv,2),size(A,1),size(A,2))
    #     D = Array{Float64, 2}(undef, 2, 2)
    #     u = Array{Float64, 1}(undef, 2)
    #     for i = 1:m
    #         for j = 1:n
    #             for k = 1:m
    #                 for l = 1:n
    #                     D[1,1] = sv[l]
    #                     D[1,2] = sv[k]
    #                     D[2,1] = sv[k]
    #                     D[2,2] = sv[l]
    #                     u[1] = Uv[i,k]*Vv[j,l]
    #                     u[2] =  -Uv[i,l]*Vv[j,k]
    #                     if (sv[l] == sv[k])
    #                         Ω = pinv(D) * u
    #                     else
    #                         Ω = D \ u
    #                     end
    #                     dU[k,l,i,j] = Ω[1]
    #                     dV[k,l,i,j] = Ω[2]
    #                 end
    #             end
    #         end
    #     end
    #     Up = zeros(typeof(A[1].partials),size(Uv))
    #     Vp = zeros(typeof(A[1].partials),size(Vv))
    #     for k = 1:size(Up,1)
    #         for l = 1:size(Up,2)
    #             for i = 1:m
    #                 for j = 1:n
    #                     Up[k,l] += dU[k,l,i,j] * A[i,j].partials
    #                 end
    #             end
    #         end
    #     end
    #     for k = 1:size(Vp,1)
    #         for l = 1:size(Vp,2)
    #             for i = 1:m
    #                 for j = 1:n
    #                     Vp[k,l] += dV[k,l,i,j] * A[i,j].partials
    #                 end
    #             end
    #         end
    #     end
    # end

    U = map(T,Uv,Up)
    s = map(T,sv,sp)
    V = map(T,Vv,Vp)

    return U, s, V
end


text:
tickname(i) = ticktype[i + 1]

const ticktype = ["BID_SIZE",                  #   0
                  "BID",                       #   1
                  "ASK",                       #   2
                  "ASK_SIZE",                  #   3
                  "LAST",                      #   4
                  "LAST_SIZE",                 #   5
                  "HIGH",                      #   6
                  "LOW",                       #   7
                  "VOLUME",                    #   8
                  "CLOSE",                     #   9
                  "BID_OPTION",                #  10
                  "ASK_OPTION",                #  11
                  "LAST_OPTION",               #  12
                  "MODEL_OPTION",              #  13
                  "OPEN",                      #  14
                  "LOW_13_WEEK",               #  15
                  "HIGH_13_WEEK",              #  16
                  "LOW_26_WEEK",               #  17
                  "HIGH_26_WEEK",              #  18
                  "LOW_52_WEEK",               #  19
                  "HIGH_52_WEEK",              #  20
                  "AVG_VOLUME",                #  21
                  "OPEN_INTEREST",             #  22
                  "OPTION_HISTORICAL_VOL",     #  23
                  "OPTION_IMPLIED_VOL",        #  24
                  "OPTION_BID_EXCH",           #  25
                  "OPTION_ASK_EXCH",           #  26
                  "OPTION_CALL_OPEN_INTEREST", #  27
                  "OPTION_PUT_OPEN_INTEREST",  #  28
                  "OPTION_CALL_VOLUME",        #  29
                  "OPTION_PUT_VOLUME",         #  30
                  "INDEX_FUTURE_PREMIUM",      #  31
                  "BID_EXCH",                  #  32
                  "ASK_EXCH",                  #  33
                  "AUCTION_VOLUME",            #  34
                  "AUCTION_PRICE",             #  35
                  "AUCTION_IMBALANCE",         #  36
                  "MARK_PRICE",                #  37
                  "BID_EFP_COMPUTATION",       #  38
                  "ASK_EFP_COMPUTATION",       #  39
                  "LAST_EFP_COMPUTATION",      #  40
                  "OPEN_EFP_COMPUTATION",      #  41
                  "HIGH_EFP_COMPUTATION",      #  42
                  "LOW_EFP_COMPUTATION",       #  43
                  "CLOSE_EFP_COMPUTATION",     #  44
                  "LAST_TIMESTAMP",            #  45
                  "SHORTABLE",                 #  46
                  "FUNDAMENTAL_RATIOS",        #  47
                  "RT_VOLUME",                 #  48
                  "HALTED",                    #  49
                  "BID_YIELD",                 #  50
                  "ASK_YIELD",                 #  51
                  "LAST_YIELD",                #  52
                  "CUST_OPTION_COMPUTATION",   #  53
                  "TRADE_COUNT",               #  54
                  "TRADE_RATE",                #  55
                  "VOLUME_RATE",               #  56
                  "LAST_RTH_TRADE",            #  57
                  "RT_HISTORICAL_VOL",         #  58
                  "IB_DIVIDENDS",              #  59
                  "BOND_FACTOR_MULTIPLIER",    #  60
                  "REGULATORY_IMBALANCE",      #  61
                  "NEWS_TICK",                 #  62
                  "SHORT_TERM_VOLUME_3_MIN",   #  63
                  "SHORT_TERM_VOLUME_5_MIN",   #  64
                  "SHORT_TERM_VOLUME_10_MIN",  #  65
                  "DELAYED_BID",               #  66
                  "DELAYED_ASK",               #  67
                  "DELAYED_LAST",              #  68
                  "DELAYED_BID_SIZE",          #  69
                  "DELAYED_ASK_SIZE",          #  70
                  "DELAYED_LAST_SIZE",         #  71
                  "DELAYED_HIGH",              #  72
                  "DELAYED_LOW",               #  73
                  "DELAYED_VOLUME",            #  74
                  "DELAYED_CLOSE",             #  75
                  "DELAYED_OPEN",              #  76
                  "RT_TRD_VOLUME",             #  77
                  "CREDITMAN_MARK_PRICE",      #  78
                  "CREDITMAN_SLOW_MARK_PRICE", #  79
                  "DELAYED_BID_OPTION",        #  80
                  "DELAYED_ASK_OPTION",        #  81
                  "DELAYED_LAST_OPTION",       #  82
                  "DELAYED_MODEL_OPTION",      #  83
                  "LAST_EXCH",                 #  84
                  "LAST_REG_TIME",             #  85
                  "FUTURES_OPEN_INTEREST",     #  86
                  "AVG_OPT_VOLUME",            #  87
                  "DELAYED_LAST_TIMESTAMP",    #  88
                  "SHORTABLE_SHARES"]          #  89


text:
@testset "GDALext" begin

	G = mat2grid(rand(Float32, 4,4));
	GMT.setproj!(G, "+proj=longlat")
	@test G.proj4 == "+proj=longlat"
	I = mat2img(rand(UInt8, 4,4));
	GMT.setproj!(I, G)
	@test I.proj4 == "+proj=longlat"

	bf = buffer([0 0], 1)
	@test bf[1].geom == Gdal.wkbPolygon
	c = centroid(bf)
	#@test c[1].data ≈ [0. 0.]
	bf2 = buffer([0.5 0], 1);
	polyunion(bf, bf2);
	Gdal.difference(bf, bf2);
	Gdal.symdifference(bf, bf2);
	@test Gdal.distance(bf, bf2) == 0.0
	@test Gdal.equals(bf, bf2) == false
	@test Gdal.disjoint(bf, bf2) == false
	@test Gdal.touches(bf, bf2) == false
	@test Gdal.crosses(bf, bf2) == false
	@test Gdal.within(bf, bf2) == false
	@test Gdal.contains(bf, bf2) == false
	@test Gdal.overlaps(bf, bf2) == true
	@test intersects(bf, [-1 -1;1 1]) == true

	Gdal.geomarea(bf);
	@test Gdal.geomlength([0 0; 1 1]) ≈ sqrt(2)
	Gdal.envelope(bf);
	Gdal.envelope3d(bf);
	Gdal.boundary(bf);
	Gdal.convexhull(bf);
	bf[1].geom = wkbLineString;
	Gdal.pointalongline(bf, 0.3);
	g1 = fromWKT("MULTIPOINT(0 0, 10 0, 10 10, 11 10)");
	g2 = delaunay(g1,2.0,true);
	@test toWKT(g2) == "MULTILINESTRING ((0 0,10 10),(0 0,10 0),(10 0,10 10))"
	D1 = mat2ds([0 0;10 0;10 10;11 10], geom=wkbMultiPoint);
	g2 = delaunay(D1,2.0,true, gdataset=true);		# Doesn't error but returns MULTILINESTRING EMPTY
	@test length(Gdal.simplify([0. 0; 1 1.1; 2.1 2], 0.2)[1]) == 4

	D1 = mat2ds([0 0; 10 0; 10 10; 11 10]);
	gdalwrite("lixo1.gmt", D1);
	D2 = gdalread("lixo1.gmt");
	@test D1 == D2
	I1 = mat2img(UInt8.([1 2 3; 4 5 6; 7 8 9]));
	gdalwrite("lixo.png", I1);
	I2 = gdalread("lixo.png");
	@test I1 == I2'		# Because the layout is different

	function test_method(f::Function, wkt1::AbstractString, wkt2::AbstractString, wkt3::AbstractString)
		geom1 = Gdal.fromWKT(wkt1)
		geom2 = Gdal.fromWKT(wkt2)
		result = f(geom1, geom2)
		@test Gdal.toWKT(result) == wkt3
		@test Gdal.toWKT(f(geom1, geom2)) == wkt3
	end

	function test_predicate(f::Function, wkt1, wkt2, result::Bool)
		geom1 = Gdal.fromWKT(wkt1)
		geom2 = Gdal.fromWKT(wkt2)
		@test f(geom1, geom2) == result
	end

	@testset "Intersection" begin
		#test_method(intersection, "POLYGON EMPTY", "POLYGON EMPTY", "GEOMETRYCOLLECTION EMPTY")
		test_method(intersection, "POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", "POINT (2 2)")
		#test_method(intersection, "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", "POLYGON((-1 1,-1 2,2 2,2 1,-1 1))",
					#"POLYGON ((0 1,0 2,2 2,2 1,0 1))")
		#=
		test_method(intersection,
					"MULTIPOLYGON(((0 0,5 10,10 0,0 0),(1 1,1 2,2 2,2 1,1 1),(100 100,100 102,102 102,102 100,100 100)))",
					"POLYGON((0 1,0 2,10 2,10 1,0 1))",
					"GEOMETRYCOLLECTION (POLYGON ((0.5 1.0,1 2,1 1,0.5 1.0)),POLYGON ((2 2,9 2,9.5 1.0,2 1,2 2)),LINESTRING (2 1,1 1),LINESTRING (1 2,2 2))")
		=#
		intersection([0 -1;0 1], [-1 -1;1 1])

		@testset "Intersects" begin
			test_predicate(Gdal.intersects, "POLYGON EMPTY", "POLYGON EMPTY", false)
			test_predicate(Gdal.intersects, "POLYGON((1 1,1 5,5 5,5 1,1 1))", "POINT(2 2)", true)
			test_predicate(Gdal.intersects, "POINT(2 2)", "POLYGON((1 1,1 5,5 5,5 1,1 1))", true)
			test_predicate(Gdal.intersects, "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", "POLYGON((1 1,1 2,2 2,2 1,1 1))", true)
			test_predicate(Gdal.intersects, "POLYGON((1 1,1 2,2 2,2 1,1 1))", "MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))", true)
		end
	end

	@test lonlat2xy([150.0 -27.0], "+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs") ≈
		[202273.912995055 7010024.033113679]  atol=1e-6
	@test lonlat2xy(mat2ds([150.0 -27.0])[1], "+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs")[1].data ≈
		[202273.912995055 7010024.033113679]  atol=1e-6

	@test xy2lonlat([202273.912995055 7010024.033113679], "+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs") ≈
		[150.0 -27.0] atol=1e-6
	@test xy2lonlat(mat2ds([202273.912995055 7010024.033113679])[1], "+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs")[1].data ≈
		[150.0 -27.0] atol=1e-6
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, DataFrames, CSV, CategoricalArrays
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.07"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	data = CSV.read(datadir("purchaseData.csv"), DataFrame; delim=",")
	#data = CSV.read("../data/purchaseData.csv", copycols = true)

end;

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
Text("Levels of Grade: $(levels(data.Grade))")

# ╔═╡ 5aab575a-3b25-11eb-3b51-afcaba4f247b
Text("Data points: $(nrow(data))")

# ╔═╡ 5ad4b776-3b25-11eb-26e5-3baf22b283ac
begin
	n = sum(.!(ismissing.(data.Grade)))
	Text("Non-missing data points: $n")
end

# ╔═╡ 5ad52f76-3b25-11eb-2b38-b7776e775586
begin
	data2 = dropmissing(data[:,[:Grade]],:Grade)

	gradeInQuestion = "E"
	indicatorVector = data2.Grade .== gradeInQuestion
	numSuccess = sum(indicatorVector)
	phat = numSuccess/n
	serr = sqrt(phat*(1-phat)/n)

	alpha = 0.05
	confidencePercent = 100*(1-alpha)
	zVal = quantile(Normal(),1-alpha/2)
	confInt = (phat - zVal*serr, phat + zVal*serr)
end

# ╔═╡ 5ae331a2-3b25-11eb-29bb-2da08e7c872d
Text("Out of $n non-missing observations $numSuccess are at level $gradeInQuestion.")

# ╔═╡ 5ae3bbd4-3b25-11eb-21e5-bf074b4ff547
Text("Hence a point estimate for the proportion of grades at level $gradeInQuestion is $phat.")

# ╔═╡ 5aef015a-3b25-11eb-2b29-3f4ba5c59a2c
Text("A $confidencePercent% confidence interval for the proprotion of level $gradeInQuestion is:\n$confInt.")

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.07"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═5aab575a-3b25-11eb-3b51-afcaba4f247b
# ╠═5ad4b776-3b25-11eb-26e5-3baf22b283ac
# ╠═5ad52f76-3b25-11eb-2b38-b7776e775586
# ╠═5ae331a2-3b25-11eb-29bb-2da08e7c872d
# ╠═5ae3bbd4-3b25-11eb-21e5-bf074b4ff547
# ╠═5aef015a-3b25-11eb-2b29-3f4ba5c59a2c
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using PastaQ
using ITensors
using Random

# Set the random seed so the results are the same
# each time it is run
Random.seed!(1234)

# Initialize the MPS state ψ = |0,0,0⟩
ψ = qubits(3)

# Apply the X gate on qubit 2
g = ("X", 2)
ψ = runcircuit(ψ, g)

# Show samples from P(x) = |⟨x|ψ⟩|²
println("Sample from |ψ⟩ = X₂|0,0,0⟩:")
display(getsamples(ψ, 3))
println()

# Custom gates
#
# There might be often problems that require quantum gates not included
# in the standard PastaQ gate set. If that is the case, a new gate can be
# added by overloading the gate(::GateName"...") function, using the format
# defined in gates.jl.
#

# Import gate so it can be overloaded
import PastaQ: gate

gate(::GateName"iX") = [
  0 im
  im 0
]

# Set the state back to |0,0,0⟩
resetqubits!(ψ)

# Show samples from P(x) = |⟨x|ψ⟩|²
println("Sample from |ψ⟩ = iX₁|0,0,0⟩:")
g = ("iX", 1)
ψ = runcircuit(ψ, g)
display(getsamples(ψ, 3))
println()