text:
include("../src/global/ImportFolders.jl")

import Config
import Client
import Connection
import Event
import Helper

function onBeforeResponse(sock::Base.TCPSocket, data::Array{Any ,1}, ans::Connection.Answer)
  println("before-response\n", data, ans)
  ans.id   = UInt8(1)
  ans.data = [UInt64(666)]
end

function onAfterResponse(data::Array{Any ,1})
  println("after-response\n", data)
end

con = Client.create(ip"127.0.0.1", Config.val(:conServerPort), true)
Event.on(con.observer, Connection.EVENT_BEFORE_RESPONSE, onBeforeResponse)
Event.on(con.observer, Connection.EVENT_AFTER_REQUEST, onAfterResponse)

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"""
TODO
"""
mutable struct HyperParameters{Fl<:AbstractFloat}
    num::Int
    names::Vector{String}
    # Poistion of the hyperparameter on the constrained and uncontrained values
    position::Dict{String,Int}
    # Minimizer vector corresponds to which hyperparameter
    minimizer_hyperparameter_position::Dict{Int,String}
    constrained_values::Vector{Fl}
    unconstrained_values::Vector{Fl}
    fixed_constrained_values::Dict{String,Fl}

    function HyperParameters{Fl}(names::Vector{String}) where Fl
        num = length(names)
        positions = Dict(names .=> collect(1:num))
        minimizer_hyperparameter_position = Dict{Int,String}()
        constrained_values = fill(Fl(NaN), num)
        unconstrained_values = fill(Fl(NaN), num)
        fixed_constrained_values = Dict{String,Fl}()
        return new{Fl}(
            num,
            names,
            positions,
            minimizer_hyperparameter_position,
            constrained_values,
            unconstrained_values,
            fixed_constrained_values,
        )
    end
end

function has_hyperparameter(model::StateSpaceModel)
    return :hyperparameters in fieldnames(typeof(model))
end

function register_unconstrained_values!(
    model::StateSpaceModel, unconstrained_values::Vector{Fl}
) where Fl
    return register_unconstrained_values!(model.hyperparameters, unconstrained_values)
end
function register_unconstrained_values!(
    hyperparameters::HyperParameters, unconstrained_values::Vector{Fl}
) where Fl
    for (i, val) in enumerate(unconstrained_values)
        update_unconstrained_value!(
            hyperparameters, hyperparameters.minimizer_hyperparameter_position[i], val
        )
    end
    return hyperparameters
end

function fill_minimizer_hyperparameter_position!(model::StateSpaceModel)
    fill_minimizer_hyperparameter_position!(model.hyperparameters)
    return model
end

function fill_minimizer_hyperparameter_position!(hyperparameters::HyperParameters)
    # Reset minimizer hyperparameter_position
    hyperparameters.minimizer_hyperparameter_position = Dict{Int,String}()
    i = 1
    for name in get_names(hyperparameters)
        if is_free(name, hyperparameters)
            hyperparameters.minimizer_hyperparameter_position[i] = name
            i += 1
        end
    end
    return hyperparameters
end

function handle_optim_initial_hyperparameters(model::StateSpaceModel)
    # Register default initial hyperparameters
    initial_hyperparameters!(model)
    # Possibly fix some user specified initial hyperparameters
    fix_hyperparameters!(model)
    # Associate the position of each index of the optimizer vector with a hyperparameter
    fill_minimizer_hyperparameter_position!(model)
    # Register the unconstrained initial_parameters
    unconstrain_hyperparameters!(model)
    # Return the unconstrained values of the initial_parameters
    return get_free_unconstrained_values(model)
end

get_hyperparameters(model::StateSpaceModel) = model.hyperparameters

"""
    get_names(model::StateSpaceModel)

Get the names of the hyperparameters registered on a `StateSpaceModel`.
"""
get_names(model::StateSpaceModel) = get_names(model.hyperparameters)

get_names(hyperparameters::HyperParameters) = hyperparameters.names

"""
    number_hyperparameters(model::StateSpaceModel)

Get the number of hyperparameters registered on a `StateSpaceModel`.
"""
function number_hyperparameters(model::StateSpaceModel)
    return number_hyperparameters(model.hyperparameters)
end

number_hyperparameters(hyperparameters::HyperParameters) = hyperparameters.num

position(str::String, hyperparameters::HyperParameters) = hyperparameters.position[str]

get_position(hyperparameters::HyperParameters) = hyperparameters.position

function get_constrained_values(model::StateSpaceModel)
    return get_constrained_values(model.hyperparameters)
end

function get_constrained_values(hyperparameters::HyperParameters)
    return hyperparameters.constrained_values
end
function get_unconstrained_values(model::StateSpaceModel)
    return get_unconstrained_values(model.hyperparameters)
end

function get_unconstrained_values(hyperparameters::HyperParameters)
    return hyperparameters.unconstrained_values
end

function get_fixed_constrained_values(hyperparameters::HyperParameters)
    return hyperparameters.fixed_constrained_values
end

function get_constrained_value(model::StateSpaceModel, str::String)
    return get_constrained_value(model.hyperparameters, str)
end

function get_constrained_value(hyperparameters::HyperParameters, str::String)
    return hyperparameters.constrained_values[position(str, hyperparameters)]
end

function update_constrained_value!(model::StateSpaceModel, str::String, value)
    return update_constrained_value!(model.hyperparameters, str, value)
end

function update_constrained_value!(hyperparameters::HyperParameters, str::String, value)
    return hyperparameters.constrained_values[position(str, hyperparameters)] = value
end

function get_unconstrained_value(model::StateSpaceModel, str::String)
    return get_unconstrained_value(model.hyperparameters, str)
end

function get_unconstrained_value(hyperparameters::HyperParameters, str::String)
    return hyperparameters.unconstrained_values[position(str, hyperparameters)]
end

function update_unconstrained_value!(model::StateSpaceModel, str::String, value)
    return update_unconstrained_value!(model.hyperparameters, str, value)
end

function update_unconstrained_value!(hyperparameters::HyperParameters, str::String, value)
    return hyperparameters.unconstrained_values[position(str, hyperparameters)] = value
end

function get_free_unconstrained_values(model::StateSpaceModel)
    return get_free_unconstrained_values(model.hyperparameters)
end

function get_free_unconstrained_values(hyperparameters::HyperParameters{Fl}) where Fl
    free_unconstrained_values = Vector{Fl}(
        undef, number_free_hyperparameters(hyperparameters)
    )
    for (k, v) in hyperparameters.minimizer_hyperparameter_position
        free_unconstrained_values[k] = get_unconstrained_value(hyperparameters, v)
    end
    return free_unconstrained_values
end

function has_hyperparameter(str::String, hyperparameters::HyperParameters)
    return haskey(hyperparameters.position, str)
end

function has_hyperparameter(str::Vector{String}, hyperparameters::HyperParameters)
    for s in str
        if !has_hyperparameter(s, hyperparameters)
            return false
        end
    end
    return true
end

"""
    set_initial_hyperparameters!(model::StateSpaceModel,
                                 initial_hyperparameters::Dict{String, <:Real})

Fill a model with user inputed initial points for hyperparameter optimzation.

# Example
```jldoctest
julia> model = LocalLevel(rand(100))
LocalLevel model

julia> get_names(model)
2-element Array{String,1}:
 "sigma2_ε"
 "sigma2_η"

julia> set_initial_hyperparameters!(model, Dict("sigma2_η" => 100.0))
LocalLevel model

julia> model.hyperparameters.constrained_values
2-element Array{Float64,1}:
 NaN
 100.0
```
"""
function set_initial_hyperparameters!(
    model::StateSpaceModel, initial_hyperparameters::Dict{String,<:Real}
)
    for (k, v) in initial_hyperparameters
        # Only update if it has not been updated before.
        if isnan(get_constrained_value(model, k))
            update_constrained_value!(model, k, v)
        end
    end
    return model
end

function is_fixed(str::String, hyperparameters::HyperParameters)
    return haskey(hyperparameters.fixed_constrained_values, str)
end

function number_free_hyperparameters(model::StateSpaceModel)
    return number_free_hyperparameters(model.hyperparameters)
end

function number_free_hyperparameters(hyperparameters::HyperParameters)
    return number_hyperparameters(hyperparameters) -
           number_fixed_hyperparameters(hyperparameters)
end

function is_free(str::String, hyperparameters::HyperParameters)
    return !is_fixed(str, hyperparameters)
end

function number_fixed_hyperparameters(model::StateSpaceModel)
    return number_fixed_hyperparameters(model.hyperparameters)
end

function number_fixed_hyperparameters(hyperparameters::HyperParameters)
    return length(hyperparameters.fixed_constrained_values)
end

"""
    fix_hyperparameters!(model::StateSpaceModel, fixed_hyperparameters::Dict)

Fixes the desired hyperparameters so that they are not considered as decision variables in
the model estimation.

# Example
```jldoctest
julia> model = LocalLevel(rand(100))
LocalLevel model

julia> get_names(model)
2-element Array{String,1}:
 "sigma2_ε"
 "sigma2_η"

julia> fix_hyperparameters!(model, Dict("sigma2_ε" => 100.0))
LocalLevel model

julia> model.hyperparameters.fixed_constrained_values
Dict{String,Float64} with 1 entry:
  "sigma2_ε" => 100.0
```
"""
function fix_hyperparameters! end

function fix_hyperparameters!(model::StateSpaceModel, fixed_hyperparameters::Dict)
    fix_hyperparameters!(model.hyperparameters, fixed_hyperparameters)
    return model
end

function fix_hyperparameters!(
    hyperparameters::HyperParameters{Fl}, fixed_hyperparameters::Dict{String,Fl}
) where Fl
    # assert we have these hyperparameters
    @assert has_hyperparameter(collect(keys(fixed_hyperparameters)), hyperparameters)
    # build dict
    hyperparameters.fixed_constrained_values = fixed_hyperparameters
    return hyperparameters
end

fix_hyperparameters!(model::StateSpaceModel) = fix_hyperparameters!(model.hyperparameters)
function fix_hyperparameters!(hyperparameters::HyperParameters{Fl}) where Fl
    for (name, value) in hyperparameters.fixed_constrained_values
        update_constrained_value!(hyperparameters, name, value)
    end
    return hyperparameters
end

"""
    isfitted(model::StateSpaceModel) -> Bool

Verify if `model` is fitted, i.e., returns `false` if there is at least one `NaN` entry in
the hyperparameters.
"""
function isfitted(model::StateSpaceModel)
    c1 = any(isnan, model.hyperparameters.unconstrained_values)
    c2 = any(isnan, model.hyperparameters.constrained_values)
    if c1 || c2
        return false
    else
        return true
    end
end

# Some special constraint functions
"""
    constrain_box!(model::StateSpaceModel, str::String, lb::Fl, ub::Fl) where Fl

Map a constrained hyperparameter ``\\psi \\in [lb, ub]`` to an unconstrained hyperparameter
``\\psi_* \\in \\mathbb{R}``.

The mapping is
``\\psi = lb + \\frac{ub - lb}{1 + \\exp(-\\psi_{*})}``
"""
function constrain_box!(model::StateSpaceModel, str::String, lb::Fl, ub::Fl) where Fl
    update_constrained_value!(
        model, str, lb + ((ub - lb) / (1 + exp(-get_unconstrained_value(model, str))))
    )
    return model
end

"""
    unconstrain_box!(model::StateSpaceModel, str::String, lb::Fl, ub::Fl) where Fl

Map an unconstrained hyperparameter ``\\psi_* \\in \\mathbb{R}`` to a constrained hyperparameter
``\\psi \\in [lb, ub]``.

The mapping is
``\\psi_* = -\\ln \\frac{ub - lb}{\\psi - lb} - 1``
"""
function unconstrain_box!(model::StateSpaceModel, str::String, lb::Fl, ub::Fl) where Fl
    update_unconstrained_value!(
        model, str, -log((ub - lb) / (get_constrained_value(model, str) - lb) - 1)
    )
    return model
end

"""
    constrain_variance!(model::StateSpaceModel, str::String)

Map a constrained hyperparameter ``\\psi \\in \\mathbb{R}^+`` to an unconstrained hyperparameter
``\\psi_* \\in \\mathbb{R}``.

The mapping is
``\\psi = \\psi_*^2``
"""
function constrain_variance!(model::StateSpaceModel, str::String)
    update_constrained_value!(model, str, (get_unconstrained_value(model, str))^2)
    return model
end

"""
    unconstrain_variance!(model::StateSpaceModel, str::String)

Map an unconstrained hyperparameter ``\\psi_{*} \\in \\mathbb{R}`` to a constrained hyperparameter
``\\psi \\in \\mathbb{R}^+``.

The mapping is ``\\psi_{*} = \\sqrt{\\psi}``.
"""
function unconstrain_variance!(model::StateSpaceModel, str::String)
    update_unconstrained_value!(model, str, sqrt(get_constrained_value(model, str)))
    return model
end

"""
    constrain_identity!(model::StateSpaceModel, str::String)

Map an constrained hyperparameter ``\\psi \\in \\mathbb{R}`` to an unconstrained hyperparameter
``\\psi_* \\in \\mathbb{R}``. This function is necessary to copy values from a location to another
inside `HyperParameters`

The mapping is
``\\psi = \\psi_*``
"""
function constrain_identity!(model::StateSpaceModel, str::String)
    update_constrained_value!(model, str, get_unconstrained_value(model, str))
    return model
end

"""
    unconstrain_identity!(model::StateSpaceModel, str::String)

Map an unconstrained hyperparameter ``\\psi_* \\in \\mathbb{R}`` to a constrained hyperparameter
``\\psi \\in \\mathbb{R}``. This function is necessary to copy values from a location to another
inside `HyperParameters`

The mapping is
``\\psi_* = \\psi``
"""
function unconstrain_identity!(model::StateSpaceModel, str::String)
    update_unconstrained_value!(model, str, get_constrained_value(model, str))
    return model
end

ext: jl


text:
using Compat
using PyCall
import PyCall: PyObject, pygui, pycall, pyexists
import Base: convert, show, download, detach

const _libcloud = PyCall.PyNULL()
const _base_driver_fns = [:list_regions]
global list_regions

if isless(Base.VERSION, v"0.5.0-")
convert{T}(::Type{Nullable{T}}, x) = convert(Nullable{T}, convert(T, x))
end

function _map_types(pymod, types)
    for t in types
        n = rsplit(string(t), '.', limit=2)[end]
        pytype_mapping(pymod[n], t)
    end
end

function __init__()
    copy!(_libcloud, pyimport_conda("libcloud", "apache-libcloud"))

    global const version = try
        convert(VersionNumber, _libcloud[:__version__])
    catch
        v"0.0" # fallback
    end
end


text:
# auxiliary function to compute correlation integrals
function integrate_χ_boxes(
    integrand :: Function,
    width     :: Float64,
    χ_tol     :: NTuple{2, Float64}
    )         :: Float64

    # map width to interval [0.0, 1.0]
    widthp = (width - 2.0 + sqrt(width * width + 4.0)) / (2.0 * width)

    # determine box size
    size = abs(widthp - 0.5)

    # perform integation for top line
    I  = hcubature_v((vv, buff) -> integrand(vv, buff), Float64[       0.0, 0.5 + size], Float64[0.5 - size,        1.0], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[0.5 - size, 0.5 + size], Float64[0.5 + size,        1.0], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[0.5 + size, 0.5 + size], Float64[       1.0         1.0], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]

    # perform integation for middle line
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[       0.0, 0.5 - size], Float64[0.5 - size, 0.5 + size], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[0.5 - size, 0.5 - size], Float64[0.5 + size, 0.5 + size], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[0.5 + size, 0.5 - size], Float64[       1.0, 0.5 + size], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]

    # perform integation for bottom line
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[       0.0,        0.0], Float64[0.5 - size, 0.5 - size], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[0.5 - size,        0.0], Float64[0.5 + size, 0.5 - size], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]
    I += hcubature_v((vv, buff) -> integrand(vv, buff), Float64[0.5 + size,        0.0], Float64[       1.0, 0.5 - size], abstol = χ_tol[1], reltol = χ_tol[2], maxevals = 10^7)[1]

    return I
end

# load correlations for different symmetries
include("correlation_su2.jl")
include("correlation_u1_dm.jl")

# generate correlation dummy
function get_χ_empty(
    symmetry :: String,
    r        :: Reduced_lattice,
    m        :: Mesh
    )        :: Vector{Matrix{Float64}}

    if symmetry == "su2"
        return get_χ_su2_empty(r, m)
    elseif symmetry == "u1-dm"
        return get_χ_u1_dm_empty(r, m)
    end
end


text:
# Common
# ------------------------------------------------------------
function MOI.get(regularizer::AbstractRegularizer, param::RawRegularizationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(regularizer.parameters)))
        error("Unrecognized parameter name: $(name) for regularizer $(typeof(regularizer)).")
    end
    return getfield(regularizer.parameters, name)
end

function MOI.set(regularizer::AbstractRegularizer, param::RawRegularizationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(regularizer.parameters)))
        error("Unrecognized parameter name: $(name) for regularizer $(typeof(regularizer)).")
    end
    setfield!(regularizer.parameters, name, value)
    return nothing
end

function add_projection_targets!(regularization::AbstractRegularization, model::MOI.AbstractOptimizer)
    ξ = regularization.ξ
    for i in eachindex(ξ)
        name = add_subscript(:ξ, i)
        set = SingleDecisionSet(1, ξ[i], NoSpecifiedConstraint(), false)
        var_index, _ = MOI.add_constrained_variable(model, set)
        set_decision!(regularization.decisions, var_index, ξ[i])
        MOI.set(model, MOI.VariableName(), var_index, name)
        regularization.projection_targets[i] = var_index
    end
    return nothing
end

function decision(::AbstractLShaped, regularization::AbstractRegularization)
    return map(regularization.ξ) do ξᵢ
        ξᵢ.value
    end
end

function objective_value(::AbstractLShaped, regularization::AbstractRegularization)
    return regularization.data.Q̃
end

function gap(lshaped::AbstractLShaped, regularization::AbstractRegularization)
    @unpack θ = lshaped.data
    @unpack Q̃ = regularization.data
    return abs(θ-Q̃)/(abs(Q̃)+1e-10)
end

function process_cut!(lshaped::AbstractLShaped, cut::AbstractHyperPlane, ::AbstractRegularization)
    return nothing
end


text:
export NDMSymm

struct NDMSymm{T,T2,T3} <: MatrixNeuralNetwork
    bare_net::T
    symm_net::T
    ∇b_mat::T2
    ∇h_mat::T2
    ∇d_mat::T2
    ∇w_mat::T2
    ∇u_mat::T2

    bare_der::T3
    symm_map
end

function NDMSymm(T::Type{<:Real}, n_in, αh, αa, permutations)
    n_symm = length(permutations)
    @assert length(first(permutations)) == n_in

    symm_net = NDM(T, n_in, αh//n_in,αa//n_in)
    bare_net = NDM(T, n_in, αh//n_in*n_symm, αa//n_in*n_symm)
    set_bare_params!(bare_net, symm_net, permutations)
    ∇b_mat, ∇h_mat, ∇d_mat, ∇w_mat, ∇u_mat = construct_∇matrices(bare_net, symm_net, n_in, αh, αa, permutations)
    bare_der = grad_cache(bare_net)

    NDMSymm(bare_net, symm_net, ∇b_mat, ∇h_mat, ∇d_mat, ∇w_mat, ∇u_mat, bare_der, permutations)
end

input_type(net::NDMSymm)         = input_type(net.bare_net)
weight_type(net::NDMSymm)        = out_type(net.bare_net)
out_type(net::NDMSymm)           = out_type(net.bare_net)
input_shape(net::NDMSymm)        =  input_shape(net.bare_net)
random_input_state(net::NDMSymm) = random_input_state(net.bare_net)
is_analytic(net::NDMSymm)        = is_analytic(net.bare_net)
weights(cnet::NDMSymm)           = cnet.symm_net

cache(net::NDMSymm)              = cache(net.bare_net)
grad_cache(net::NDMSymm)         = grad_cache(net.symm_net)
weight_tuple(net::NDMSymm, args...) = weight_tuple(weights(net), args...)
update!(opt, cnet::NDMSymm, Δ, state=nothing) = (res = update!(opt, weights(cnet), weights(Δ), state);
                                            set_bare_params!(cnet.bare_net, cnet.symm_net, cnet.symm_map);
                                            res)

@inline (net::NDMSymm)(args...) = net.bare_net(args...)
function logψ_and_∇logψ!(∇lnψ_symm, net::NDMSymm, args...)
    lnψ, ∇lnψ = logψ_and_∇logψ!(net.bare_der, net.bare_net, args...)

    symmetrize_∇logψ_NDM!(∇lnψ_symm, ∇lnψ, net)
    lnψ, ∇lnψ_symm
end

function symmetrize_∇logψ_NDM!(∇lnψ_symm, ∇lnψ, net)
    mul!(∇lnψ_symm.b_μ, net.∇b_mat, ∇lnψ.b_μ)
    mul!(∇lnψ_symm.b_λ, net.∇b_mat, ∇lnψ.b_λ)
    mul!(∇lnψ_symm.h_μ, net.∇h_mat, ∇lnψ.h_μ)
    mul!(∇lnψ_symm.h_λ, net.∇h_mat, ∇lnψ.h_λ)
    mul!(∇lnψ_symm.d_λ, net.∇d_mat, ∇lnψ.d_λ)
    mul!(vec(∇lnψ_symm.w_μ), net.∇w_mat, vec(∇lnψ.w_μ))
    mul!(vec(∇lnψ_symm.w_λ), net.∇w_mat, vec(∇lnψ.w_λ))
    mul!(vec(∇lnψ_symm.u_μ), net.∇u_mat, vec(∇lnψ.u_μ))
    mul!(vec(∇lnψ_symm.u_λ), net.∇u_mat, vec(∇lnψ.u_λ))
    ∇lnψ_symm
end

function set_bare_params!(bare_net, symm_net, permutations)
    # Load a bit
    n_symm = length(permutations)
    n_in   = length(symm_net.b_μ)
    αh   = size(symm_net.w_μ,1)
    αa   = size(symm_net.u_μ,1)

    # Local bias
    symm_net.b_μ .= sum(symm_net.b_μ)/length(symm_net.b_μ)
    symm_net.b_λ .= sum(symm_net.b_λ)/length(symm_net.b_λ)
    bare_net.b_μ .=  symm_net.b_μ
    bare_net.b_λ .=  symm_net.b_λ

    # Hidden biases
    for (i, (h_μ_i, h_λ_i)) = enumerate(zip(symm_net.h_μ, symm_net.h_λ))
        bare_net.h_μ[((i-1)*n_symm+1):(i*n_symm)] .=  h_μ_i
        bare_net.h_λ[((i-1)*n_symm+1):(i*n_symm)] .=  h_λ_i
    end

    # Ancilla biases
    for (i, d_λ_i) = enumerate(symm_net.d_λ)
        bare_net.d_λ[((i-1)*n_symm+1):(i*n_symm)] .=  d_λ_i
    end

    # Hidden connections
    for f=1:αh
        for (j,perm)=enumerate(permutations)
            for (i, i_p)=enumerate(perm)
                bare_net.w_μ[j + (f-1)*n_symm, i_p] = symm_net.w_μ[f,i]
                bare_net.w_λ[j + (f-1)*n_symm, i_p] = symm_net.w_λ[f,i]
            end
        end
    end

    # Ancilla connections
    for f=1:αa
        for (j,perm)=enumerate(permutations)
            for (i, i_p)=enumerate(perm)
                bare_net.u_μ[j + (f-1)*n_symm, i_p] = symm_net.u_μ[f,i]
                bare_net.u_λ[j + (f-1)*n_symm, i_p] = symm_net.u_λ[f,i]
            end
        end
    end
    bare_net
end

function construct_∇matrices(bare_net, symm_net, n_in, αh, αa, permutations)
    # Load a bit
    n_symm = length(permutations)

    # Local bias matrix
    ∇b_mat   = ones(eltype(bare_net.u_λ), n_in, n_in)./n_in

    # Hidden bias matrix
    ∇h_mat   = zeros(eltype(bare_net.h_λ), αh, αh*n_symm)
    for s=1:αh
        ∇h_mat[s, (s-1)*n_symm+1:s*n_symm] .= 1
    end

    # Ancillary bias matrix
    ∇d_mat   = zeros(eltype(bare_net.d_λ), αa, αa*n_symm)
    for s=1:αa
        ∇d_mat[s, (s-1)*n_symm+1:s*n_symm] .= 1
    end

    ∇w_mat   = zeros(eltype(bare_net.w_λ), αh*n_in, αh*n_symm*n_in)
    lin_ind_bare  = LinearIndices(bare_net.w_λ)
    lin_ind_symm  = LinearIndices(symm_net.w_λ)
    for f=1:αh
        for (j,perm)=enumerate(permutations)
            for (i, i_p)=enumerate(perm)
                row = lin_ind_symm[f,i]
                col = lin_ind_bare[j + (f-1)*n_symm, i_p]
                ∇w_mat[row,col] = 1
            end
        end
    end

    ∇u_mat   = zeros(eltype(bare_net.u_λ), αa*n_in, αa*n_symm*n_in)
    lin_ind_bare  = LinearIndices(bare_net.u_λ)
    lin_ind_symm  = LinearIndices(symm_net.u_λ)
    for f=1:αa
        for (j,perm)=enumerate(permutations)
            for (i, i_p)=enumerate(perm)
                row = lin_ind_symm[f,i]
                col = lin_ind_bare[j + (f-1)*n_symm, i_p]
                ∇u_mat[row,col] = 1
            end
        end
    end
    ∇b_mat, ∇h_mat, ∇d_mat, ∇w_mat, ∇u_mat
end

######


text:
struct FixedOutput <: AbstractDeviceFormulation end

function _check_device_formulation(
    ::Type{D},
) where {D <: Union{AbstractDeviceFormulation, PSY.Device}}
    if !isconcretetype(D)
        throw(
            ArgumentError(
                "The device model must contain only concrete types, $(D) is an Abstract Type",
            ),
        )
    end
end

"""
Establishes the model for a particular device specified by type. Uses the keyword argument
feedforward to enable passing values between operation model at simulation time

# Arguments

-`::Type{D}`: Power System Device Type
-`::Type{B}`: Abstract Device Formulation

# Accepted Key Words

  - `feedforward::Array{<:AbstractAffectFeedforward}` : use to pass parameters between models

# Example

thermal_gens = DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
"""
mutable struct DeviceModel{D <: PSY.Device, B <: AbstractDeviceFormulation}
    feedforwards::Vector{<:AbstractAffectFeedforward}
    use_slacks::Bool
    duals::Vector{DataType}
    services::Vector{ServiceModel}
    time_series_names::Dict{Type{<:TimeSeriesParameter}, String}
    attributes::Dict{String, Any}

    function DeviceModel(
        ::Type{D},
        ::Type{B};
        feedforwards=Vector{AbstractAffectFeedforward}(),
        use_slacks=false,
        duals=Vector{DataType}(),
        time_series_names=get_default_time_series_names(D, B),
        attributes=get_default_attributes(D, B),
    ) where {D <: PSY.Device, B <: AbstractDeviceFormulation}
        _check_device_formulation(D)
        _check_device_formulation(B)
        new{D, B}(
            feedforwards,
            use_slacks,
            duals,
            Vector{ServiceModel}(),
            time_series_names,
            attributes,
        )
    end
end

get_component_type(
    ::DeviceModel{D, B},
) where {D <: PSY.Device, B <: AbstractDeviceFormulation} = D
get_formulation(
    ::DeviceModel{D, B},
) where {D <: PSY.Device, B <: AbstractDeviceFormulation} = B
get_feedforwards(m::DeviceModel) = m.feedforwards
get_services(m::DeviceModel) = m.services
get_services(::Nothing) = nothing
get_use_slacks(m::DeviceModel) = m.use_slacks
get_duals(m::DeviceModel) = m.duals
get_time_series_names(m::DeviceModel) = m.time_series_names
get_attributes(m::DeviceModel) = m.attributes
get_attribute(m::DeviceModel, key::String) = get(m.attributes, key, nothing)
DeviceModelForBranches = DeviceModel{<:PSY.Branch, <:AbstractDeviceFormulation}

function _set_model!(
    dict::Dict,
    model::DeviceModel{D, B},
) where {D <: PSY.Device, B <: AbstractDeviceFormulation}
    key = Symbol(D)
    if haskey(dict, key)
        @warn "Overwriting $(D) existing model"
    end
    dict[key] = model
    return
end

has_service_model(model::DeviceModel) = !isempty(get_services(model))


text:
#=
def logx(x, n):
    for i in range(0, n):
        x = math.log(x) if x>0 else 0
    return x
=#

function logx(x, n)             # iterative logarithm function
    for i in 0:n-1
    	x = x > 0.0 ? log(x) : 0.0
    end
    return x
end


text:
_mpi_rank() = nothing

@init @require MPIClusterManagers="e7922434-ae4b-11e9-05c5-9780451d2c66" begin

    using .MPIClusterManagers: MPI, start_main_loop, TCP_TRANSPORT_ALL, MPI_TRANSPORT_ALL

    """
    init_MPI_workers()

    Initialize MPI processes as Julia workers. Should be called from all MPI
    processes, and will only return on the master process.

    `transport` should be `"MPI"` or `"TCP"`, which is by default read from the
    environment variable `JULIA_MPI_TRANSPORT`, and otherwise defaults to `"TCP"`.

    If CUDA is loaded and functional in the Main module, additionally calls
    [`assign_GPU_workers()`](@ref)
    """
    function init_MPI_workers(;
        stdout_to_master = false,
        stderr_to_master = false,
        transport = get(ENV,"JULIA_MPI_TRANSPORT","TCP"),
        print_info = true,
    )

        if !MPI.Initialized()
            MPI.Init()
        end
        size = MPI.Comm_size(MPI.COMM_WORLD)
        rank = MPI.Comm_rank(MPI.COMM_WORLD)

        if size>1
            # workers don't return from this call:
            start_main_loop(
                Dict("TCP"=>TCP_TRANSPORT_ALL,"MPI"=>MPI_TRANSPORT_ALL)[transport],
                stdout_to_master=stdout_to_master,
                stderr_to_master=stderr_to_master
            )

            if @isdefined(CUDA) && CUDA.functional()
                assign_GPU_workers(;print_info=false)
            end

            print_info && proc_info()

        end

    end

    _mpi_rank() = MPI.Initialized() ? MPI.Comm_rank(MPI.COMM_WORLD) : nothing

end

"""
    assign_GPU_workers()

Assign each Julia worker process a unique GPU using `CUDA.device!`.
Workers may be distributed across different hosts, and each host can have
multiple GPUs.
"""
function assign_GPU_workers(;print_info=true)
    @everywhere @eval Main using Distributed, CMBLensing
    master_uuid = @isdefined(CUDA) ? CUDA.uuid(device()) : nothing
    accessible_gpus = Dict(map(workers()) do id
        @eval Main @fetchfrom $id begin
            ds = CUDA.devices()
            # put master's GPU last so we don't double up on it unless we need to
            $id => sort((CUDA.deviceid.(ds) .=> CUDA.uuid.(ds)), by=(((k,v),)->v==$master_uuid ? Inf : k))
        end
    end)
    claimed = Set()
    assignments = Dict(map(workers()) do myid
        for (gpu_id, gpu_uuid) in accessible_gpus[myid]
            if !(gpu_uuid in claimed)
                push!(claimed, gpu_uuid)
                return myid => gpu_id
            end
        end
        error("Can't assign a unique GPU to every worker, process $myid has no free GPUs left.")
    end)
    @everywhere workers() device!($assignments[myid()])
    print_info && proc_info()
end

"""
    proc_info()

Returns string showing info about available processes.
"""
function proc_info()
    @eval Main using Distributed
    lines = @eval Main map(procs()) do id
        @fetchfrom id begin
            info = ["myid = $id"]
            !isnothing(CMBLensing._mpi_rank()) && push!(info, "mpi-rank = $(CMBLensing._mpi_rank())")
            push!(info, "host = $(gethostname())")
            @isdefined(CUDA) && push!(info, "device = $(sprint(io->show(io, MIME("text/plain"), CUDA.device()))) $(split(string(CUDA.uuid(CUDA.device())),'-')[1]))")
            " ("*join(info, ", ")*")"
        end
    end
    @info join(["Processes:"; lines], "\n")
end

# a ProgressMeter that can be advanced from any remote worker

struct DistributedProgress <: ProgressMeter.AbstractProgress
    channel :: RemoteChannel{Channel{Any}}
end

function DistributedProgress(args...; kwargs...)
    pbar = Progress(args...; kwargs...)
    channel = RemoteChannel(()->Channel(), 1)
    @async begin
        while (x = take!(channel)) != nothing
            func, args, kwargs = x
            func(pbar, args...; kwargs...)
        end
        finish!(pbar)
    end
    DistributedProgress(channel)
end

ProgressMeter.next!(pbar::DistributedProgress, args...; kwargs...) = put!(pbar.channel, (ProgressMeter.next!, args, kwargs))
ProgressMeter.update!(pbar::DistributedProgress, args...; kwargs...) = put!(pbar.channel, (ProgressMeter.update!, args, kwargs))
ProgressMeter.finish!(pbar::DistributedProgress, args...; kwargs...) = (put!(pbar.channel, nothing); close(pbar.channel))


text:
# Makie seems to have trouble with SVectors, so just collect them into Vectors
function convert_arguments(T::Type{<:LineSegments},
                           V::SVector{NP, Polytope{K,P,N,PT}}) where {NP,K,P,N,PT}
    return convert_arguments(T, collect(V))
end

function convert_arguments(T::Type{<:GLMakieMesh},
                           V::SVector{NP, Polytope{K,P,N,PT}}) where {NP,K,P,N,PT}
    return convert_arguments(T, collect(V))
end

# LineSegment
function convert_arguments(LS::Type{<:LineSegments}, l::LineSegment)
    return convert_arguments(LS, collect(vertices(l)))
end

function convert_arguments(LS::Type{<:LineSegments},
                           L::Vector{LineSegment{T}}) where {T}
    points = Vector{T}(undef, 2length(L))
    for i in eachindex(L)
        points[2i-1] = L[i][1]
        points[2i  ] = L[i][2]
    end
    return convert_arguments(LS, points)
end

# QuadraticSegment
function convert_arguments(LS::Type{<:LineSegments}, q::QuadraticSegment{T}) where {T}
    points = Vector{T}(undef, 2*(plot_nonlinear_subdivisions-1))
    r = LinRange(0, 1, plot_nonlinear_subdivisions)
    for i = 1:plot_nonlinear_subdivisions-1
        points[2i-1] = q(r[i  ])
        points[2i  ] = q(r[i+1])
    end
    return convert_arguments(LS, points)
end

function convert_arguments(LS::Type{<:LineSegments}, Q::Vector{<:QuadraticSegment})
    return convert_arguments(LS, vcat([convert_arguments(LS, q)[1] for q in Q]...))
end

# General Polytopes
#
# LineSegments
function convert_arguments(LS::Type{<:LineSegments}, poly::Polytope)
    return convert_arguments(LS, edges(poly))
end

function convert_arguments(LS::Type{<:LineSegments}, P::Vector{<:Polytope})
    return convert_arguments(LS, vcat([convert_arguments(LS, p)[1] for p ∈  P]...))
end

# Mesh
# 2D Polygons can be triangulated using fan triangulation.
# All others need to be triangulated based on the number of nonlinear subdivisions
function convert_arguments(M::Type{<:GLMakieMesh}, tri::Triangle)
    vertices = [coordinates(v) for v in ridges(tri)]
    face = [1 2 3]
    return convert_arguments(M, vertices, face)
end

function convert_arguments(M::Type{<:GLMakieMesh},
                           tris::Vector{Triangle{Point{Dim,T}}}) where {Dim,T}
    points = Vector{Vec{Dim,T}}(undef, 3length(tris))
    for i in eachindex(tris)
        tri = tris[i]
        points[3i-2] = coordinates(tri[1])
        points[3i-1] = coordinates(tri[2])
        points[3i  ] = coordinates(tri[3])
    end
    faces = zeros(Int64, length(tris), 3)
    k = 1
    for i in 1:length(tris), j = 1:3
        faces[i, j] = k
        k += 1
    end
    return convert_arguments(M, points, faces)
end

function convert_arguments(M::Type{<:GLMakieMesh}, poly::Polygon{N, Point{2,T}}) where {N,T}
    return convert_arguments(M, triangulate(poly))
end

function convert_arguments(M::Type{<:GLMakieMesh}, poly::Polytope{2})
    return convert_arguments(M, triangulate(poly, Val(plot_nonlinear_subdivisions)))
end

function convert_arguments(M::Type{<:GLMakieMesh}, P::Vector{<:Polytope{2}})
    triangles = [ begin
                     if poly isa Polygon{N, Point{2,T}} where {N,T}
                         triangulate(poly)
                     elseif poly isa Triangle
                         Vec(poly)
                     else
                         triangulate(poly, Val(plot_nonlinear_subdivisions))
                     end
                  end
                  for poly in P]
    return convert_arguments(M, reduce(vcat, triangles))
end

function convert_arguments(M::Type{<:GLMakieMesh}, poly::Polytope{3})
    return convert_arguments(M, faces(poly))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
@testsuite "random" AT->begin
    @testset "rand" begin  # uniform
        for T in (Int8, Float32, Float64, Int64, Int32,
                  Complex{Float32}, Complex{Float64},
                  Complex{Int64}, Complex{Int32}), d in (10, (10,10))
            A = AT{T}(undef, d)
            B = copy(A)
            rand!(A)
            rand!(B)
            @test Array(A) != Array(B)

            rng = if AT <: AbstractGPUArray
                GPUArrays.default_rng(AT)
            else
                Random.default_rng()
            end
            Random.seed!(rng)
            Random.seed!(rng, 1)
            rand!(rng, A)
            Random.seed!(rng, 1)
            rand!(rng, B)
            @test all(Array(A) .== Array(B))
        end

        A = AT{Bool}(undef, 5)
        rand!(A)
        @test true in Array(A)
        @test false in Array(A)
    end

    @testset "randn" begin  # uniform
        for T in (Float32, Float64), d in (2, (2,2))
            A = AT{T}(undef, d)
            B = copy(A)
            randn!(A)
            randn!(B)
            @test !any(Array(A) .== Array(B))

            rng = if AT <: AbstractGPUArray
                GPUArrays.default_rng(AT)
            else
                Random.default_rng()
            end
            Random.seed!(rng)
            Random.seed!(rng, 1)
            randn!(rng, A)
            Random.seed!(rng, 1)
            randn!(rng, B)
            @test all(Array(A) .== Array(B))
        end
    end
end


text:
using Test
using Polyhedra
using JuMP
const MOIT = MOI.Test
const MOIB = MOI.Bridges

@testset "Continuous Linear problems with VRepOptimizer" begin
    optimizer = VRepOptimizer{Float64}()
    @testset "SolverName" begin
        @test MOI.get(optimizer, MOI.SolverName()) == "VRep"
    end
    cache = MOIU.UniversalFallback(Polyhedra._MOIModel{Float64}())
    cached = MOIU.CachingOptimizer(cache, optimizer)
    bridged = MOIB.full_bridge_optimizer(cached, Float64)
    config = MOIT.TestConfig(duals=false)
    MOIT.contlineartest(bridged, config,
                        # linear8a and linear12 will be solved by https://github.com/jump-dev/MathOptInterface.jl/pull/702
                        ["linear8a", "linear12", "partial_start"])
end
@testset "simplex chebyshev center with $T" for T in [Float64, Rational{BigInt}]
    h = HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0) ∩ HyperPlane([1, 1], 1)
    lib = DefaultLibrary{T}(VRepOptimizer{T})
    poly = polyhedron(h, lib)
    center, radius = chebyshevcenter(poly)
    @test center ≈ [1/2, 1/2]
    @test radius ≈ 1/2
    @test dim(poly) == 1
end


text:
function update!(dQ, Q, rka, rkb, dt, slow_δ, slow_dQ, slow_scaling)
  @inbounds @loop for i = (1:length(Q);
                           (blockIdx().x - 1) * blockDim().x + threadIdx().x)
    if slow_δ !== nothing
      dQ[i] += slow_δ * slow_dQ[i]
    end
    Q[i] += rkb * dt * dQ[i]
    dQ[i] *= rka
    if slow_scaling !== nothing
      slow_dQ[i] *= slow_scaling
    end
  end
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
include("datalinks.jl")
register(DataDep("PokesWildTypes",
    """
    Dataset: Hidden state foraging task, interleaved barrier, pokes WT.
    Author: Sarra et. al. (Champalimaud Research)
    """,
    pokespath,
    "a5ef1b92cd5593811e3a89a2304e7d7662056c4c1b5890b8522d218b2bf0ee15",
    post_fetch_method = x -> mv(x, x*".csv")
));

const PokesDatasetFolder = datadep"PokesWildTypes"

register(DataDep("StreaksWildTypes",
    """
    Dataset: Hidden state foraging task, interleaved barrier, streaks WT.
    Author: Sarra et. al. (Champalimaud Research)
    """,
    streakspath,
    "9693c801d5069c8c09ea3836b17e0003b207c5dd9a16bf339cb05cad71024026",
    post_fetch_method = x -> mv(x, x*".csv")
));

const StreaksDatasetFolder = datadep"StreaksWildTypes"


text:
module MLJTuning

## METHOD EXPORT

# defined in tuned_models.jl:
export TunedModel

# defined in strategies/:
export Explicit, Grid

# defined in learning_curves.jl:
export learning_curve!, learning_curve

## METHOD IMPORT

import MLJBase
using MLJBase
using RecipesBase
using Distributed
import ComputationalResources: CPU1, CPUProcesses,
    CPUThreads, AbstractResource
using Random

## INCLUDE FILES

include("utilities.jl")
include("tuning_strategy_interface.jl")
include("tuned_models.jl")
include("ranges.jl")
include("strategies/explicit.jl")
include("strategies/grid.jl")
include("plotrecipes.jl")
include("learning_curves.jl")

end


text:
type RuleSource
  rulename::Symbol
  args::Vector{Any}
  comments::Vector{AbstractString}
  source::Vector{AbstractString}
  function RuleSource(rulename::Symbol, args=Any[])
    new(rulename, args, Vector{AbstractString}[], Vector{AbstractString}[])
  end
end

function output_generator(io::IO, genname::Symbol, description, rules::Vector{RuleSource})
  println(io, "@generator $(genname) begin")
  println(io)
  println(io, "generates: [\"$(description)\"]")
  println(io)
  for rule in rules
    output_rule(io, rule)
    println(io)
  end
  println(io, "end")
end

function output_rule(io::IO, rule::RuleSource)
  for comment in rule.comments
    println(io, "# $(comment)")
  end
  args = join(map(a->string(a), rule.args),", ")
  print(io, "$(rule.rulename)($(args)) = ")
  for sourceline in rule.source
    println(io, sourceline)
  end
end

function build_rule_start(node::ASTNode)
  rule = RuleSource(node.rulename)
  rule.comments = node.source
  push!(rule.source, "begin")
  for warning in node.warnings
    push!(rule.source, "  warn(\"$(warning) in $(node.rulename)\")")
  end
  rule
end

function build_rule_end(rule::RuleSource, node::ASTNode)
  push!(rule.source, "end")
end

function build_rule_shortform_start(node::ASTNode)
  if !isempty(node.warnings)
    rule = build_rule_start(node)
  else
    rule = RuleSource(node.rulename)
    rule.comments = node.source
  end
  rule
end

function build_rule_shortform_end(rule::RuleSource, node::ASTNode)
  if !isempty(node.warnings)
    build_rule_end(rule, node)
  end
end

function build_called_rulename(callnode::ASTNode)
  if isempty(callnode.refs)
    callnode.rulename
  else
    @assert length(callnode.refs)==1
    refnode = callnode.refs[1]
    refnode.rulename
  end
end

function build_called_child_rulename(node::ASTNode, func::Symbol)
    funcchildren = filter(child -> child.func == func, node.children)
    @assert length(funcchildren) == 1
    build_called_rulename(funcchildren[1])
end

escape_rule_name(name) = replace(name, r"[^a-zA-Z0-9]", "_")

function assign_rulenames(node::ASTNode, rulenameprefix="")
  rulename = escape_rule_name(rulenameprefix * (isempty(rulenameprefix) ? "" : " ") * describe_ast_node(node))
  if isempty(rulenameprefix)
    # special case: empty prefix at root node, so name this the "start" rule as it will be the entry point to the generator
    node.rulename = Symbol("start")
  else
    node.rulename = Symbol(rulename)
  end
  for enumchild in enumerate(node.children)
    assign_rulenames(enumchild[2], rulename * " $(enumchild[1])") # note we use "canonical" rulename rather than actual rule name used for this node as prefix to children
  end
end


text:
@reexport module ProgressiveHedgingSolvers

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using Parameters
using JuMP
using StochasticPrograms
using StochasticPrograms: _WS
using StochasticPrograms: AbstractScenarioProblems, ScenarioProblems, DScenarioProblems
using StochasticPrograms: LQSolver, getsolution, getobjval, getredcosts, getduals, status, QPSolver, get_solver, loadLP
using StochasticPrograms: Execution, Serial, Synchronous, Asynchronous
using StochasticPrograms: PenaltyTerm, Quadratic, Linearized, InfNorm, ManhattanNorm, initialize_penaltyterm!, update_penaltyterm!, solve_penalized!
using MathProgBase
using ProgressMeter

import Base: show, put!, wait, isready, take!, fetch
import StochasticPrograms: StructuredModel, internal_solver, optimize_structured!, fill_solution!, solverstr

const MPB = MathProgBase

export
    ProgressiveHedgingSolver,
    ProgressiveHedging,
    Fixed,
    Adaptive,
    Serial,
    Synchronous,
    Asynchronous,
    Quadratic,
    Linearized,
    InfNorm,
    ManhattanNorm,
    Crash,
    StructuredModel,
    optimsolver,
    optimize_structured!,
    fill_solution!

# Include files
include("types/types.jl")
include("penalties/penalization.jl")
include("execution/execution.jl")
include("solver.jl")
include("spinterface.jl")

end # module


text:
struct DummyCallback end

@testset "simulation" begin
    roadway = gen_straight_roadway(1, 500.0)

    models = Dict{Int, DriverModel}()
    models[1] = IntelligentDriverModel(k_spd = 1.0, v_des = 10.0)
    models[2] = IntelligentDriverModel(k_spd = 1.0, v_des = 5.0)

    veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 0.0), roadway, 5.)
    veh1 = Vehicle(veh_state, VehicleDef(), 1)
    veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 70.0), roadway, 5.)
    veh2 = Vehicle(veh_state, VehicleDef(), 2)

    scene = Scene()
    push!(scene, veh1)
    push!(scene, veh2)

    n_steps = 40
    dt = 0.1
    rec = SceneRecord(n_steps, dt)
    simulate!(rec, scene, roadway, models, n_steps)

    reset_hidden_states!(models)

      # initializing vehicles too close
    veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 0.0), roadway, 10.)
    veh1 = Vehicle(veh_state, VehicleDef(), 1)
    veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 5.0), roadway, 5.)
    veh2 = Vehicle(veh_state, VehicleDef(), 2)

    scene = Scene()
    push!(scene, veh1)
    push!(scene, veh2)

    rec = SceneRecord(n_steps, dt)
    simulate!(rec, scene, roadway, models, 10, (CollisionCallback(),))

    @test_throws ErrorException simulate!(rec, scene, roadway, models, 10, (DummyCallback(),))

    # collision right from start
    veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 0.0), roadway, 10.)
    veh1 = Vehicle(veh_state, VehicleDef(), 1)
    veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 1.0), roadway, 5.)
    veh2 = Vehicle(veh_state, VehicleDef(), 2)

    scene = Scene()
    push!(scene, veh1)
    push!(scene, veh2)

    rec = SceneRecord(n_steps, dt)
    simulate!(rec, scene, roadway, models, 10, (CollisionCallback(),))
end

@testset "trajdata simulation" begin
  roadway = get_test_roadway()
  trajdata = get_test_trajdata(roadway)

  veh_state = VehicleState(Frenet(roadway[LaneTag(1,1)], 6.0), roadway, 10.)
  ego = Vehicle(veh_state, VehicleDef(), 2)
  model = ProportionalSpeedTracker()
  dt = get_timestep(trajdata)
  rec = SceneRecord(3, dt)
  simulate!(rec, model, ego.id, trajdata, roadway, 1, 2)
  @test findfirst(ego.id, rec[0]) != nothing
end


text:
# ---
# title: 611. Valid Triangle Number
# id: problem611
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/valid-triangle-number/description/>
# hidden: true
# ---
#
# Given an array consists of non-negative integers, your task is to count the
# number of triplets chosen from the array that can make triangles if we take
# them as side lengths of a triangle.
#
# **Example 1:**
#
#
#
#     Input: [2,2,3,4]
#     Output: 3
#     Explanation:
#     Valid combinations are:
#     2,3,4 (using the first 2)
#     2,3,4 (using the second 2)
#     2,2,3
#
#
# **Note:**
#
#   1. The length of the given array won't exceed 1000.
#   2. The integers in the given array are in the range of [0, 1000].
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
for dtype in ("", "Dist")
    for stype in ("", "Sparse")
        for (elty, ext) in ((:Float32, :s),
                            (:Float64, :d),
                            (:ComplexF32, :c),
                            (:ComplexF64, :z))
            mat = Symbol(dtype, stype, "Matrix")
            @eval begin
                function solve!(A::$mat{$elty}, B::$mat{$elty})
                    ElError(ccall(($(string("ElLinearSolve", dtype, stype, "_", ext)), libEl), Cuint,
                        (Ptr{Cvoid}, Ptr{Cvoid}),
                        A.obj, B.obj))
                    return B
                end
            end
        end
    end
end


text:
"""
Abstract type that allows users to pass in their preferred Algorithm
"""
abstract AbstractMultirouteFlowAlgorithm

"""
Forces the multiroute_flow function to use the Kishimoto algorithm.
"""
type KishimotoAlgorithm <: AbstractMultirouteFlowAlgorithm
end

"""
Forces the multiroute_flow function to use the Extended Multiroute Flow algorithm.
"""
type ExtendedMultirouteFlowAlgorithm <: AbstractMultirouteFlowAlgorithm
end

# Methods when the number of routes is more than the connectivity
# 1) When using Boykov-Kolmogorov as a flow subroutine
# 2) Other flow algorithm
function empty_flow{T<:Real}(
  capacity_matrix::AbstractArray{T, 2},     # edge flow capacities
  flow_algorithm::BoykovKolmogorovAlgorithm # keyword argument for algorithm
  )
  n = size(capacity_matrix, 1)
  return zero(T), zeros(T, n, n), zeros(T, n)
end
# 2) Other flow algorithm
function empty_flow{T<:Real}(
  capacity_matrix::AbstractArray{T, 2},     # edge flow capacities
  flow_algorithm::AbstractFlowAlgorithm     # keyword argument for algorithm
  )
  n = size(capacity_matrix, 1)
  return zero(T), zeros(T, n, n)
end

# Method for Kishimoto algorithm
function multiroute_flow{T<:Real}(
  flow_graph::DiGraph,                   # the input graph
  source::Int,                           # the source vertex
  target::Int,                           # the target vertex
  capacity_matrix::AbstractArray{T, 2},  # edge flow capacities
  flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
  mrf_algorithm::KishimotoAlgorithm,     # keyword argument for algorithm
  routes::Int                            # keyword argument for routes
  )
  return kishimoto(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
end

## Methods for Extended Multiroute Flow Algorithm
#1 When the breaking points are not already known
function multiroute_flow{T<:Real, R<:Real}(
  flow_graph::DiGraph,                            # the input graph
  source::Int,                                    # the source vertex
  target::Int,                                    # the target vertex
  capacity_matrix::AbstractArray{T, 2},           # edge flow capacities
  flow_algorithm::AbstractFlowAlgorithm,          # keyword argument for algorithm
  mrf_algorithm::ExtendedMultirouteFlowAlgorithm, # keyword argument for algorithm
  routes::R                                       # keyword argument for routes
  )
  return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
end
#2 When the breaking points are already known
#2-a Output: flow value (paired with the associated restriction)
function multiroute_flow{T<:Real, R<:Real}(
  breakingpoints::Vector{Tuple{T, T, Int}},       # vector of breaking points
  routes::R                                       # keyword argument for routes
  )
  return intersection(breakingpoints, routes)
end
#2-b Output: flow value, flows(, labels)
function multiroute_flow{T1<:Real, T2<:Real, R<:Real}(
  breakingpoints::Vector{Tuple{T1, T1, Int}}, # vector of breaking points
  routes::R,                                # keyword argument for routes
  flow_graph::DiGraph,                      # the input graph
  source::Int,                              # the source vertex
  target::Int,                              # the target vertex
  capacity_matrix::AbstractArray{T2, 2} =    # edge flow capacities
    DefaultCapacity(flow_graph);
  flow_algorithm::AbstractFlowAlgorithm  =  # keyword argument for algorithm
    PushRelabelAlgorithm()
  )
  x, f = intersection(breakingpoints, routes)

  # For other cases, capacities need to be Floats
  if !(T2<:AbstractFloat)
    capacity_matrix = convert(AbstractArray{Float64, 2}, capacity_matrix)
  end

  return maximum_flow(flow_graph, source, target, capacity_matrix,
                      algorithm = flow_algorithm, restriction = x)
end

"""
The generic multiroute_flow function will output three kinds of results:

- When the number of routes is 0 or non-specified, the set of breaking points of
the multiroute flow is returned.
- When the input is limited to a set of breaking points and a route value k,
only the value of the k-route flow is returned
- Otherwise, a tuple with 1) the maximum flow and 2) the flow matrix. When the
max-flow subroutine is the Boykov-Kolmogorov algorithm, the associated mincut is
returned as a third output.

When the input is a network, it requires the following arguments:

- flow_graph::DiGraph                   # the input graph
- source::Int                           # the source vertex
- target::Int                           # the target vertex
- capacity_matrix::AbstractArray{T, 2}  # edge flow capacities with T<:Real
- flow_algorithm::AbstractFlowAlgorithm # keyword argument for flow algorithm
- mrf_algorithm::AbstractFlowAlgorithm  # keyword argument for multiroute flow algorithm
- routes::R<:Real                       # keyword argument for the number of routes

When the input is only the set of (breaking) points and the number of route,
it requires the following arguments:

- breakingpoints::Vector{Tuple{T, T, Int}},    # vector of breaking points
- routes::R<:Real,                             # number of routes

When the input is the set of (breaking) points, the number of routes,
and the network descriptors, it requires the following arguments:

- breakingpoints::Vector{Tuple{T1, T1, Int}} # vector of breaking points (T1<:Real)
- routes::R<:Real                            # number of routes
- flow_graph::DiGraph                        # the input graph
- source::Int                                # the source vertex
- target::Int                                # the target vertex
- capacity_matrix::AbstractArray{T2, 2}      # optional edge flow capacities (T2<:Real)
- flow_algorithm::AbstractFlowAlgorithm      # keyword argument for algorithm

The function defaults to the Push-relabel (classical flow) and Kishimoto
(multiroute) algorithms. Alternatively, the algorithms to be used can also
be specified through  keyword arguments. A default capacity of 1 is assumed
for each link if no capacity matrix is provided.

The mrf_algorithm keyword is inforced to Extended Multiroute Flow
in the following cases:

- The number of routes is non-integer
- The number of routes is 0 or non-specified

### Usage Example :
(please consult the  max_flow section for options about flow_algorithm
and capacity_matrix)

```julia

# Create a flow-graph and a capacity matrix
flow_graph = DiGraph(8)
flow_edges = [
    (1, 2, 10), (1, 3, 5),  (1, 4, 15), (2, 3, 4),  (2, 5, 9),
    (2, 6, 15), (3, 4, 4),  (3, 6, 8),  (4, 7, 16), (5, 6, 15),
    (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6),  (7, 8, 10)
]
capacity_matrix = zeros(Int, 8, 8)
for e in flow_edges
    u, v, f = e
    add_edge!(flow_graph, u, v)
    capacity_matrix[u, v] = f
end

# Run default multiroute_flow with an integer number of routes = 2
f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 2)

# Run default multiroute_flow with a noninteger number of routes = 1.5
f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 1.5)

# Run default multiroute_flow for all the breaking points values
points = multiroute_flow(flow_graph, 1, 8, capacity_matrix)
# Then run multiroute flow algorithm for any positive number of routes
f, F = multiroute_flow(points, 1.5)
f = multiroute_flow(points, 1.5, valueonly = true)

# Run multiroute flow algorithm using Boykov-Kolmogorov algorithm as max_flow routine
f, F, labels = multiroute_flow(flow_graph, 1, 8, capacity_matrix,
               algorithm = BoykovKolmogorovAlgorithm(), routes = 2)

```
"""

function multiroute_flow{T<:Real, R<:Real}(
  flow_graph::DiGraph,                    # the input graph
  source::Int,                            # the source vertex
  target::Int,                            # the target vertex
  capacity_matrix::AbstractArray{T, 2} =  # edge flow capacities
    DefaultCapacity(flow_graph);
  flow_algorithm::AbstractFlowAlgorithm  =    # keyword argument for algorithm
    PushRelabelAlgorithm(),
  mrf_algorithm::AbstractMultirouteFlowAlgorithm  =    # keyword argument for algorithm
    KishimotoAlgorithm(),
  routes::R = 0              # keyword argument for number of routes (0 = all values)
  )
  # a flow with a set of 1-disjoint pathes is a classical max-flow
  (routes == 1) &&
  return maximum_flow(flow_graph, source, target, capacity_matrix, flow_algorithm)

  # routes > λ (connectivity) → f = 0
  λ = maximum_flow(flow_graph, source, target, DefaultCapacity(flow_graph),
                                             algorithm = flow_algorithm)[1]
  (routes > λ) && return empty_flow(capacity_matrix, flow_algorithm)

  # For other cases, capacities need to be Floats
  if !(T<:AbstractFloat)
    capacity_matrix = convert(AbstractArray{Float64, 2}, capacity_matrix)
  end

  # Ask for all possible values (breaking points)
  (routes == 0) &&
  return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm)

  # The number of routes is a float → EMRF
  (R <: AbstractFloat) &&
  return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)

  # Other calls
  return multiroute_flow(flow_graph, source, target, capacity_matrix,
                               flow_algorithm, mrf_algorithm, routes)
end


text:
include("euler/euler.jl")
using .Numbers: is_s_number
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result = 0
    for i ∈ 2:n if is_s_number(i, i * i) result += i * i end end
    return result
end

compute(100)

compute(1_000_000)

@benchmark compute(1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using BSON: @load
using Flux
using Plots
# get the first layer parameters for influence analysis
@load "best.bson" model
beta = model.layers[1].W # get first layer betas
z = maximum(abs.(beta),dims=1);
heatmap(z, xlabel="statistic", title="Importance of inputs, bright=high, dark=low")
savefig("ImportanceOfStatistics.png")


text:
using Test
using SafeTestsets

@testset verbose=true "Qecsim.jl" begin
    @safetestset "error.jl" begin include("error.jl") end
    @safetestset "PauliTools.jl" begin include("PauliTools.jl") end
    @safetestset "Model.jl" begin include("Model.jl") end
    @safetestset "models/BasicModels.jl" begin include("models/BasicModels.jl") end
    @safetestset "models/GenericModels/simpleerrormodel.jl" begin
        include("models/generic/simpleerrormodel.jl")
    end
    @safetestset "models/GenericModels/naivedecoder.jl" begin
        include("models/generic/naivedecoder.jl")
    end
    @safetestset "App.jl" begin include("App.jl") end
end


text:
"""
    Domain

Abstract type for different domains.
"""
abstract type Domain end

"""
    ClosedIntervalDomain

Abstract type for interval domains.
"""
abstract type AbstractIntervalDomain <: Domain end

"""
    ClosedIntervalDomain(left, right)

Interval domain.
"""
struct ClosedIntervalDomain{T} <: AbstractIntervalDomain
    left::T
    right::T
end

"""
    PeriodicIntervalDomain(left, right)

Periodic interval domain.
"""
struct PeriodicIntervalDomain{T} <: AbstractIntervalDomain
    left::T
    right::T
end
PeriodicIntervalDomain(a, b) = PeriodicIntervalDomain(promote(a, b)...)


text:
@testset "Constructors" begin
    @testset "Vector" begin
        @testset "Dense Vector" begin
            test_frule(GBVector, rand(-10.0:0.05:10.0, 10))
            test_rrule(GBVector, rand(-10.0:0.05:10.0, 10))
        end
        @testset "Sparse Vector" begin
            test_frule(GBVector, [1,3,5] ⊢ NoTangent(), rand(-10.0:0.05:10.0, 3))
            test_rrule(GBVector, [1,4,5] ⊢ NoTangent(), rand(-10.0:0.05:10.0, 3))
        end
    end
    @testset "Matrix" begin
        @testset "Dense Matrix and Vector" begin
            M = rand(-10.0:0.05:10.0, 10, 20)
            v = rand(-10.0:0.05:10.0, 10)
            test_frule(GBMatrix, M)
            test_rrule(GBMatrix, M)
            test_frule(GBMatrix, v)
            test_rrule(GBMatrix, v)
        end
    end
end


text:
typealias VecVec{T} Vector{Vector{T}}
typealias VecVecVec{T} Vector{Vector{Vector{T}}}

nested_array_size{T<:AbstractVector}(x::Vector{T}) = [nested_array_size(v) for v in x]
nested_array_size(x::AbstractVector) = length(x)

nested_array_length{T<:AbstractVector}(x::Vector{T}) = sum([nested_array_length(x) for x in x])
nested_array_length(x::AbstractVector) = length(x)

nested_array_type(T::DataType, d::Int) = d == 0 ? T : Vector{nested_array_type(T,d-1)}

function nested_array_depth(T::DataType)
    d = 0
    while issubtype(T,AbstractArray)
        T = T.parameters[1]
        d += 1
    end
    d
end

nested_array_depth(x::AbstractArray) = nested_array_depth(typeof(x))

function nested_array_eltype(T::DataType)
    while issubtype(T,AbstractArray)
        T = eltype(T)
    end
    T::DataType
end

nested_array_eltype(x::AbstractArray) = nested_array_eltype(typeof(x))

function rand_nested_array(d,len)
    if d == 1
        return rand(rand(len))
    else
        T = nested_array_type(Float64,d-1)
        T[rand_nested_array(d-1,len) for _ in 1:rand(len)]
    end
end

reduce_size{T<:Integer}(shp::Vector{T}) = sum(shp)
reduce_size{T}(shp::VecVec{T}) = [reduce_size(s) for s in shp]

function reduce_ptr{T<:Integer}(ptr::VecVec{T})
    p = Array{Int}(length(ptr)+1)
    for i in 1:length(ptr)
        p[i] = ptr[i][1]
    end
    p[end] = ptr[end][end]
    p
end

reduce_ptr{T<:Integer}(ptr::VecVecVec{T}) = [reduce_ptr(p) for p in ptr]
reduce_ptr{T<:Integer}(ptr::Vector{T}) = 1

@generated function expand_size(shp)
    d = nested_array_depth(shp)
    quote
        shps = Any[shp]
        for _ in 1:$d
            shp = nested_array_size(shp)
            insert!(shps,1,shp)
        end
        tuple(shps...)
    end
end

expand_size(shp::Integer) = (shp,)

function index_nested_array!{T<:AbstractVector}(x::Vector{T}, val, k)
    ptr = Array{Int}(length(x)+1)
    ptr[1] = k[1]
    for (i,v) in enumerate(x)
        ptr[i+1] = (k[1] += length(v))
        val[ptr[i]:(ptr[i+1]-1)] = v
    end
    ptr
end

index_nested_array!{T<:AbstractVector}(x::VecVec{T},val,k) = [index_nested_array!(v,val,k) for v in x]
index_nested_array!(x,val) = index_nested_array!(x,val,[1])

@generated function index_nested_array(x)
    Tv = nested_array_eltype(x)
    quote
        val = Array{$Tv}(nested_array_length(x))
        ptr = index_nested_array!(x,val)
        val,ptr
    end
end

function combine_size(shp)
    n = length(shp)
    d = length(shp[1])
    T = typeof(shp[1]).types
    cshp = Any[length(shp)]
    for j in 1:d
        push!(cshp,T[j][s[j] for s in shp])
    end
    tuple(cshp...)
end

add_offset{T<:Integer}(x::Vector{T},i,k) = begin k[1] = x[end]; x+i-1 end
add_offset{T<:Integer}(x::VecVec{T},i,k) = [add_offset(v,i,k) for v in x]

function combine_ptr(ptr)
    k = [1]
    T = typeof(ptr[1])
    T[add_offset(p,k[1],k) for p in ptr]
end

type RaggedArray{Tv,N,Tp,Ts} <: DenseVector{Tv}
    val::Vector{Tv}
    ptr::Tp
    shp::Ts
end

@generated function RaggedArray(x::AbstractVector)
    Tv = nested_array_eltype(x)
    N = nested_array_depth(x)
    quote
        val,ptr = index_nested_array(x)
        shp = nested_array_size(x)
        ptrs = Any[ptr]
        shps = Any[expand_size(shp)]
        for _ in 1:$(N-1)
            ptr = reduce_ptr(ptr)
            shp = reduce_size(shp)
            insert!(ptrs,1,ptr)
            insert!(shps,1,expand_size(shp))
        end
        ptr = tuple(ptrs...)
        shp = tuple(shps...)
        RaggedArray{$Tv,$N,typeof(ptr),typeof(shp)}(val,ptr,shp)
    end
end

function RaggedArray{Tv,N}(x::RaggedArray{Tv,N}...)
    val = vcat([v.val for v in x]...)

    ptr = Any[combine_ptr([v.ptr[i] for v in x]) for i in 2:N]
    insert!(ptr,1,reduce_ptr(ptr[1]))
    insert!(ptr,1,1)
    ptr = tuple(ptr...)

    shp = Any[combine_size([v.shp[i] for v in x]) for i in 1:N]
    insert!(shp,1,(reduce_size(shp[1][end]),))
    shp = tuple(shp...)

    RaggedArray{Tv,N,typeof(ptr),typeof(shp)}(val,ptr,shp)
end

Base.eltype{Tv}(a::RaggedArray{Tv}) = Tv
Base.length(ra::RaggedArray) = length(ra.val)
Base.ndims(::RaggedArray) = 1
Base.eachindex(ra::RaggedArray) = eachindex(ra.val)
Base.linearindexing{T<:RaggedArray}(::Type{T}) = Base.LinearFast()
Base.size(ra::RaggedArray) = size(ra.val)
Base.size(ra::RaggedArray, i::Integer) = size(ra.val, i)

Base.similar{Tv,N,Tp,Ts}(ra::RaggedArray{Tv,N,Tp,Ts}) = RaggedArray{Tv,N,Tp,Ts}(similar(val),copy(ra.ptr),copy(ra.shp))

_ref(sym, i) = Expr(:ref, sym, i)  # e.g. A[i_1,...,i_N]
_nestedref(sym, i) = _ref(sym, i)
_nestedref(sym, i...) = _nestedref(_nestedref(sym, i[1]), i[2:end]...)  # A[i_1][i_2]...[i_n]
_nindices(sym, n) = [symbol(sym, :_, i) for i = 1:n]  # [i_1,...,i_n]

@generated function len{Tv,N}(ra::RaggedArray{Tv,N}, d::Integer, i::Integer...)
    n = length(i)                                               # number of indices
    ind = [_ref(:i,j) for j in 1:n]                             # [i_1,...,i_n]
    len = _nestedref(:(ra.shp), :d, n+1, ind...)                # ra.shp[d][2][i[1]], ra.shp[d][3][i[1]][i[2]], etc...
    quote
        $len::Int
    end
end

Base.getindex(ra::RaggedArray, i...) = error("not implemented")
Base.getindex(ra::RaggedArray, i::Integer) = ra.val[i]

@generated function Base.getindex{Tv,N}(ra::RaggedArray{Tv,N}, i::Integer...)
    d = length(i)                                           # number of indices
    d > N && throw(BoundsError)
    ind = [_ref(:i,j) for j in 1:d]                         # [i_1,...,i_d]
    offset = _nestedref(:(ra.ptr), d, ind[end:-1:2]...)     # ra.ptr[1][i[2]], ra.ptr[2][i[3]][i[2]], ra.ptr[3][i[4]][i[3]][i[2]], etc...
    quote
      ra.val[$offset + i[1] - 1]                            # no bounds check ?
    end
end

@generated function Base.getindex{Tv,N}(ra::RaggedArray{Tv,N}, ::Colon, i::Integer...)
    d = length(i)+1
    d > N && throw(BoundsError)
    ind = reverse([_ref(:i,j) for j in 1:length(i)])
    a = _nestedref(:(ra.ptr), d, ind...)
    ind[end] = :($(ind[end])+1)
    b = _nestedref(:(ra.ptr), d, ind...)
    quote
      ra.val[$a:$b-1]
    end
end

# # multiply by weights w and then sum over last dimension
# function sum{Tv,N,Tp,Ts}(ra::RaggedArray{Tv,N,Tp,Ts}, w::AbstractVector)
#     raw = RaggedArray{Tv,N,Tp,Ts}(ra.val.*w,ra.ptr,ra.shp)
#     [sum(raw[:,i]) for i in 1:len(raw,N)]
# end

max_stars_repo_name: johncsnyder/Element


text:
module CUDAnative

using CUDAapi
using CUDAdrv

using LLVM
using LLVM.Interop

using Adapt
using TimerOutputs
using DataStructures

using Libdl

## global state

const toolkit_version = Ref{VersionNumber}()
version() = toolkit_version[]

# version compatibility
const target_support = Ref{Vector{VersionNumber}}()
const ptx_support = Ref{Vector{VersionNumber}}()

# paths
const libdevice = Ref{Union{String,Dict{VersionNumber,String}}}()
const libcudadevrt = Ref{String}()
const nvdisasm = Ref{String}()
const ptxas = Ref{String}()

## source code includes

include("utils.jl")

# needs to be loaded _before_ the compiler infrastructure, because of generated functions
include("device/tools.jl")
include("device/pointer.jl")
include("device/array.jl")
include("device/cuda.jl")
include("device/llvm.jl")
include("device/runtime.jl")

include("init.jl")

include("compiler.jl")
include("execution.jl")
include("exceptions.jl")
include("reflection.jl")

include("deprecated.jl")

## initialization

const __initialized__ = Ref(false)
functional() = __initialized__[]

function __init__()
    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    silent = parse(Bool, get(ENV, "JULIA_CUDA_SILENT", "false")) || precompiling
    verbose = parse(Bool, get(ENV, "JULIA_CUDA_VERBOSE", "false"))

    # if any dependent GPU package failed, expect it to have logged an error and bail out
    if !CUDAdrv.functional()
        verbose && @warn "CUDAnative.jl did not initialize because CUDAdrv.jl failed to"
        return
    end

    try
        ## target support

        # LLVM.jl

        llvm_version = LLVM.version()
        llvm_targets, llvm_isas = llvm_support(llvm_version)

        # Julia

        julia_llvm_version = Base.libllvm_version
        if julia_llvm_version != llvm_version
            error("LLVM $llvm_version incompatible with Julia's LLVM $julia_llvm_version")
        end

        if v"8.0" <= llvm_version < v"9.0"
            # NOTE: corresponding functionality in irgen.jl
            silent || @warn "LLVM pre-9.0 detected, disabling debug info emission for CUDA kernels"
        end

        # CUDA

        toolkit_dirs = find_toolkit()
        toolkit_version[] = find_toolkit_version(toolkit_dirs)
        if toolkit_version[] <= v"9"
            silent || @warn "CUDAnative.jl only supports CUDA 9.0 or higher (your toolkit provides CUDA $(toolkit_version[]))"
        end

        cuda_targets, cuda_isas = cuda_support(CUDAdrv.version(), toolkit_version[])

        target_support[] = sort(collect(llvm_targets ∩ cuda_targets))
        isempty(target_support[]) && error("Your toolchain does not support any device target")

        ptx_support[] = sort(collect(llvm_isas ∩ cuda_isas))
        isempty(ptx_support[]) && error("Your toolchain does not support any PTX ISA")

        @debug("CUDAnative supports devices $(verlist(target_support[])); PTX $(verlist(ptx_support[]))")

        let val = find_libdevice(target_support[], toolkit_dirs)
            val === nothing && error("Your CUDA installation does not provide libdevice")
            libdevice[] = val
        end

        let val = find_libcudadevrt(toolkit_dirs)
            val === nothing && error("Your CUDA installation does not provide libcudadevrt")
            libcudadevrt[] = val
        end

        let val = find_cuda_binary("nvdisasm", toolkit_dirs)
            val === nothing && error("Your CUDA installation does not provide the nvdisasm binary")
            nvdisasm[] = val
        end

        let val = find_cuda_binary("ptxas", toolkit_dirs)
            val === nothing && error("Your CUDA installation does not provide the ptxas binary")
            ptxas[] = val
        end

        ## actual initialization

        __init_compiler__()

        resize!(thread_contexts, Threads.nthreads())
        fill!(thread_contexts, nothing)
        CUDAdrv.initializer(maybe_initialize)

        __initialized__[] = true
    catch ex
        # don't actually fail to keep the package loadable
        if !silent
            if verbose
                @error "CUDAnative.jl failed to initialize" exception=(ex, catch_backtrace())
            else
                @info "CUDAnative.jl failed to initialize, GPU functionality unavailable (set JULIA_CUDA_SILENT or JULIA_CUDA_VERBOSE to silence or expand this message)"
            end
        end
    end
end

verlist(vers) = join(map(ver->"$(ver.major).$(ver.minor)", sort(collect(vers))), ", ", " and ")

function llvm_support(version)
    @debug("Using LLVM v$version")

    # https://github.com/JuliaGPU/CUDAnative.jl/issues/428
    if version >= v"8.0" && VERSION < v"1.3.0-DEV.547"
        error("LLVM 8.0 requires a newer version of Julia")
    end

    InitializeAllTargets()
    haskey(targets(), "nvptx") ||
        error("""
            Your LLVM does not support the NVPTX back-end.

            This is very strange; both the official binaries
            and an unmodified build should contain this back-end.""")

    target_support = sort(collect(CUDAapi.devices_for_llvm(version)))

    ptx_support = CUDAapi.isas_for_llvm(version)
    push!(ptx_support, v"6.0") # JuliaLang/julia#23817
    ptx_support = sort(collect(ptx_support))

    @debug("LLVM supports devices $(verlist(target_support)); PTX $(verlist(ptx_support))")
    return target_support, ptx_support
end

function cuda_support(driver_version, toolkit_version)
    @debug("Using CUDA driver v$driver_version and toolkit v$toolkit_version")

    # the toolkit version as reported contains major.minor.patch,
    # but the version number returned by libcuda is only major.minor.
    toolkit_version = VersionNumber(toolkit_version.major, toolkit_version.minor)
    if toolkit_version > driver_version
        @warn("""CUDA $(toolkit_version.major).$(toolkit_version.minor) is not supported by
                 your driver (which supports up to $(driver_version.major).$(driver_version.minor))""")
    end

    driver_target_support = CUDAapi.devices_for_cuda(driver_version)
    toolkit_target_support = CUDAapi.devices_for_cuda(toolkit_version)
    target_support = sort(collect(driver_target_support ∩ toolkit_target_support))

    driver_ptx_support = CUDAapi.isas_for_cuda(driver_version)
    toolkit_ptx_support = CUDAapi.isas_for_cuda(toolkit_version)
    ptx_support = sort(collect(driver_ptx_support ∩ toolkit_ptx_support))

    @debug("CUDA driver supports devices $(verlist(driver_target_support)); PTX $(verlist(driver_ptx_support))")
    @debug("CUDA toolkit supports devices $(verlist(toolkit_target_support)); PTX $(verlist(toolkit_ptx_support))")

    return target_support, ptx_support
end

end


text:
using Base.Filesystem
import MIMEs

"""
    open_in_default_browser(url)

Open a URL in the ambient default browser.

Note: this was copied from Pluto.jl.
"""
function open_in_default_browser(url::AbstractString)::Bool
    try
        if Sys.isapple()
            Base.run(`open $url`)
            true
        elseif Sys.iswindows() || detectwsl()
            Base.run(`cmd.exe /s /c start "" /b $url`)
            true
        elseif Sys.islinux()
            Base.run(`xdg-open $url`)
            true
        else
            false
        end
    catch ex
        false
    end
end

"""
    update_and_close_viewers!(wss::Vector{HTTP.WebSockets.WebSocket})

Take a list of viewers, i.e. WebSocket connections from a client,
send a message with data "update" to each of them (to trigger a page reload),
then close the connection. Finally, empty the list since all connections are
closing anyway and clients will re-connect from the re-loaded page.
"""
function update_and_close_viewers!(wss::Vector{HTTP.WebSockets.WebSocket})
    foreach(wss) do wsi
        try
            write(wsi, "update")
            close(wsi)
        catch
        end
    end
    empty!(wss)
    return nothing
end

"""
    file_changed_callback(f_path::AbstractString)

Function reacting to the change of the file at `f_path`. Is set as callback for
the file watcher.
"""
function file_changed_callback(f_path::AbstractString)
    if VERBOSE[]
        println("ℹ [LiveUpdater]: Reacting to change in file '$f_path'...")
    end
    if endswith(f_path, ".html")
        # if html file, update viewers of this file only
        # check the viewer still exists otherwise may error
        # see issue https://github.com/asprionj/LiveServer.jl/issues/95
        if haskey(WS_VIEWERS, f_path)
            update_and_close_viewers!(WS_VIEWERS[f_path])
        end
    else
        # otherwise (e.g. modification to a CSS file), update all viewers
        foreach(update_and_close_viewers!, values(WS_VIEWERS))
    end
    return nothing
end

"""
    get_fs_path(req_path::AbstractString)

Return the filesystem path corresponding to a requested path, or an empty
String if the file was not found.
"""
function get_fs_path(req_path::AbstractString)::String
    uri = HTTP.URI(req_path)
    # first element after the split is **always** "/" --> 2:end
    r_parts = HTTP.URIs.unescapeuri.(split(lstrip(uri.path, '/'), '/'))
    fs_path = joinpath(r_parts...)

    if !isempty(CONTENT_DIR[])
        fs_path = joinpath(CONTENT_DIR[], fs_path)
    end

    isfile(fs_path) && return fs_path

    tmp = joinpath(fs_path, "index.html")
    isfile(tmp)     && return tmp

    # content of the dir will be shown
    isdir(fs_path)  && return fs_path

    # 404 will be shown
    return ""
end

"""
    append_slash(url::AbstractString) -> url′::AbstractString

Append `/` to the path part of `url`; i.e., transform `a/b` to `a/b/` and `/a/b?c=d` to
`/a/b/?c=d`.
"""
function append_slash(url_str::AbstractString)
    uri = HTTP.URI(url_str)
    return string(endswith(uri.path, "/") ? uri : merge(uri; path = uri.path * "/"))
end

"""
    lstrip_cdir(s)

Discard the 'CONTENT_DIR' part (passed via `dir=...`) of a path.
"""
lstrip_cdir(s::AbstractString) =
    lstrip(s[nextind(s, length(CONTENT_DIR[])):end], ['/', '\\'])

"""
    get_dir_list(dir::AbstractString) -> index_page::AbstractString

Generate list of content at path `dir`.
"""
function get_dir_list(dir::AbstractString)
    list   = readdir(dir; join=true, sort=true)
    io     = IOBuffer()
    predir = ifelse(isempty(CONTENT_DIR[]), "", "[$(append_slash(CONTENT_DIR[]))]")
    sdir   = predir * lstrip_cdir(dir)

    write(io, """
        <!DOCTYPE HTML>

            <meta charset="utf-8">
            <meta name="viewport" content="width=device-width, initial-scale=1">
            <link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/spcss">
            <title>Directory listing</title>

              Directory listing

              $(sdir)

        """
    )

    list_files = [f for f in list if isfile(f)]
    list_dirs  = [d for d in list if d ∉ list_files]

    for fname in list_files
        link = lstrip_cdir(fname)
        name = splitdir(fname)[end]
        post = ifelse(islink(fname), " @", "")
        write(io, """

$(name)$(post)

            """
        )
    end
    for fdir in list_dirs
        link = lstrip_cdir(fdir)
        name = splitdir(fdir)[end]
        pre  = "📂 "
        post = ifelse(islink(fdir), " @", "")
        write(io, """

$(pre)$(name)$(post)

            """
        )
    end
    write(io, """

            🏠 root

            💻 LiveServer.jl

        """
    )
    return String(take!(io))
end

"""
    serve_file(fw, req::HTTP.Request; inject_browser_reload_script = true)

Handler function for serving files. This takes a file watcher, to which files
to be watched can be added, and a request (e.g. a path entered in a tab of the
browser), and converts it to the appropriate file system path.

The cases are as follows:
1. the path corresponds exactly to a file. If it's a html-like file,
    LiveServer will try injecting the reloading `<script>` (see file
    `client.html`) at the end, just before the `</body>` tag.
2. the path corresponds to a directory in which there is an `index.html`,
    same action as (1) assuming the `index.html` is implicit.
3. the path corresponds to a directory in which there is not an `index.html`,
    list the directory contents.
4. not (1,2,3), a 404 is served.

All files served are added to the file watcher, which is responsible to check
whether they're already watched or not. Finally the file is served via a 200
(successful) response. If the file does not exist, a response with status 404
and message is returned.
"""
function serve_file(
            fw, req::HTTP.Request;
            inject_browser_reload_script::Bool = true,
            allow_cors::Bool = false
        )::HTTP.Response

    ret_code = 200
    fs_path  = get_fs_path(req.target)

    # if get_fs_path returns an empty string, there's two cases:
    # 1. the path is a directory without an `index.html` --> list dir
    # 2. otherwise serve a 404 (see if there's a dedicated 404 path,
    #     otherwise just use a basic one).
    if isempty(fs_path)

        if req.target == "/"
            index_page = get_dir_list(".")
            return HTTP.Response(200, index_page)
        end

        ret_code = 404
        # Check if /404/ or /404.html exists and serve that as a body
        for f in ("/404/", "/404.html")
            maybe_path = get_fs_path(f)
            if !isempty(maybe_path)
                fs_path = maybe_path
                break
            end
        end

        # If still not found a body, return a generic error message
        if isempty(fs_path)
            return HTTP.Response(404, """
                404: file not found. Perhaps you made a typo in the URL,
                or the requested file has been deleted or renamed.
                """
            )
        end
    end

    if isdir(fs_path)
        index_page = get_dir_list(fs_path)
        return HTTP.Response(200, index_page)
    end

    # In what follows, fs_path points to a file
    # --> html-like: try to inject reload-script
    # --> other: just get the browser to show it
    #
    ext     = lstrip(last(splitext(fs_path)), '.') |> string
    content = read(fs_path, String)

    # build the response with appropriate mime type (this is inspired from Mux
    # https://github.com/JuliaWeb/Mux.jl/blob/master/src/examples/files.jl)
    content_type = let
        mime_from_ext = MIMEs.mime_from_extension(ext, nothing)
        if mime_from_ext !== nothing
            MIMEs.contenttype_from_mime(mime_from_ext)
        else
            HTTP.sniff(content)
        end
    end

    # avoid overly-specific text types so the browser can try rendering
    # all text-like documents instead of offering to download all files
    m = match(r"^text\/(.*?);(.*)$", content_type)
    if m !== nothing
        text_type = m.captures[1]
        if text_type ∉ ("html", "javascript", "css")
            content_type = "text/plain;$(m.captures[2])"
        end
    end
    plain = "text/plain; charset=utf8"
    for p  in ("application/toml", "application/x-sh")
        content_type = replace(content_type, p => plain)
    end

    # if html-like, try adding the browser-sync script to it
    inject_reload = inject_browser_reload_script && (
            startswith(content_type, "text/html") ||
            startswith(content_type, "application/xhtml+xml")
        )

    if inject_reload
        end_body_match = match(r"", content)
        if end_body_match === nothing
            # no  tag found, trying to add the reload script at the
            # end. This may fail.
            content *= BROWSER_RELOAD_SCRIPT
        else
            end_body = prevind(content, end_body_match.offset)
            # reconstruct the page with the reloading script
            io = IOBuffer()
            write(io, SubString(content, 1:end_body))
            write(io, BROWSER_RELOAD_SCRIPT)
            write(io, SubString(content, nextind(content, end_body):lastindex(content)))
            content = take!(io)
        end
    end

    headers = ["Content-Type" => content_type]

    if allow_cors
        push!(headers, "Access-Control-Allow-Origin" => "*")
    end
    resp         = HTTP.Response(ret_code, content)
    resp.headers = HTTP.mkheaders(headers)

    # add the file to the file watcher
    watch_file!(fw, fs_path)

    # return the response
    return resp
end

"""
    ws_upgrade(http::HTTP.Stream)

Upgrade the HTTP request in the stream to a websocket.
"""
function ws_upgrade(http::HTTP.Stream)
    # adapted from HTTP.WebSockets.upgrade; note that here the upgrade will always
    # have  the right format as it always triggered by after a Response
    HTTP.setstatus(http, 101)
    HTTP.setheader(http, "Upgrade" => "websocket")
    HTTP.setheader(http, "Connection" => "Upgrade")
    key = HTTP.header(http, "Sec-WebSocket-Key")
    HTTP.setheader(http, "Sec-WebSocket-Accept" => HTTP.WebSockets.accept_hash(key))
    HTTP.startwrite(http)

    io = http.stream
    return HTTP.WebSockets.WebSocket(io; server=true)
end

"""
    ws_tracker(ws::HTTP.WebSockets.WebSocket, target::AbstractString)

Adds the websocket connection to the viewers in the global dictionary `WS_VIEWERS` to the entry
corresponding to the targeted file.
"""
function ws_tracker(ws::HTTP.WebSockets.WebSocket, target::AbstractString)
    # add to list of html files being "watched"
    # NOTE: this file always exists because the query is generated just after
    # serving it
    fs_path = get_fs_path(target)

    # if the file is already being viewed, add ws to it (e.g. several tabs)
    # otherwise add to dict
    if fs_path ∈ keys(WS_VIEWERS)
        push!(WS_VIEWERS[fs_path], ws)
    else
        WS_VIEWERS[fs_path] = [ws]
    end

    try
        # NOTE: browsers will drop idle websocket connections so this effectively
        # forces the websocket to stay open until it's closed by LiveServer (and
        # not by the browser) upon writing a `update` message on the websocket.
        # See update_and_close_viewers
        while isopen(ws.io)
            sleep(0.1)
        end
    catch err
        # NOTE: there may be several sources of errors caused by the precise moment
        # at which the user presses CTRL+C and after what events. In an ideal world
        # we would check that none of these errors have another source but for now
        # we make the assumption it's always the case (note that it can cause other
        # errors than InterruptException, for instance it can cause errors due to
        # stream not being available etc but these all have the same source).
        # - We therefore do not propagate the error but merely store the information
        # that there was a forcible interruption of the websocket so that the
        # interruption can be guaranteed to be propagated.
        WS_INTERRUPT[] = true
    end
    return nothing
end

"""
    serve(filewatcher; host="127.0.0.1", port=8000, dir="", verbose=false, coreloopfun=(c,fw)->nothing, inject_browser_reload_script::Bool = true, launch_browser::Bool = false, allow_cors::Bool = false)

Main function to start a server at `http://host:port` and render what is in the current
directory. (See also [`example`](@ref) for an example folder).

# Arguments

- `filewatcher` is a file watcher implementing the API described for [`SimpleWatcher`](@ref) (which also is the default) and messaging the viewers (via WebSockets) upon detecting file changes.
- `port` is an integer between 8000 (default) and 9000.
- `dir` specifies where to launch the server if not the current working directory.
- `verbose` is a boolean switch to make the server print information about file changes and connections.
- `coreloopfun` specifies a function which can be run every 0.1 second while the liveserver is going; it takes two arguments: the cycle counter and the filewatcher. By default the coreloop does nothing.
- `launch_browser=false` specifies whether to launch the ambient browser at the localhost URL or not.
- `allow_cors::Bool=false` will allow cross origin (CORS) requests to access the server via the "Access-Control-Allow-Origin" header.
- `preprocess_request=identity`: specifies a function which can transform a request before a response is returned; its only argument is the current request.

# Example

```julia
LiveServer.example()
serve(host="127.0.0.1", port=8080, dir="example", verbose=true, launch_browser=true)
```

You should then see the `index.html` page from the `example` folder being rendered. If you change the file, the browser will automatically reload the
page and show the changes.
"""
function serve(fw::FileWatcher=SimpleWatcher(file_changed_callback);
               host::String="127.0.0.1", port::Int=8000, dir::AbstractString="", verbose::Bool=false,
               coreloopfun::Function=(c, fw)->nothing,
               preprocess_request=identity,
               inject_browser_reload_script::Bool = true,
               launch_browser::Bool = false,
               allow_cors::Bool = false)

    8000 ≤ port ≤ 9000 || throw(ArgumentError("The port must be between 8000 and 9000."))
    setverbose(verbose)

    if !isempty(dir)
        isdir(dir) || throw(ArgumentError("The specified dir '$dir' is not recognised."))
        CONTENT_DIR[] = dir
    end

    start(fw)

    # make request handler
    req_handler = HTTP.RequestHandlerFunction() do req
        req = preprocess_request(req)
        serve_file(fw, req; inject_browser_reload_script = inject_browser_reload_script, allow_cors = allow_cors)
    end

    server = Sockets.listen(parse(IPAddr, host), port)
    url = "http://$(host == string(Sockets.localhost) ? "localhost" : host):$port"
    println("✓ LiveServer listening on $url/ ...\n  (use CTRL+C to shut down)")
    @async HTTP.listen(host, port;
                       server=server, readtimeout=0, reuse_limit=0) do http::HTTP.Stream
        # reuse_limit=0 ensures that there won't be an error if killing and restarting the server.
        if HTTP.WebSockets.is_upgrade(http.message)
            # upgrade to websocket
            ws = ws_upgrade(http)
            # add to list of viewers and keep open until written to
            ws_tracker(ws, http.message.target)
        else
            # handle HTTP request
            HTTP.handle(req_handler, http)
        end
    end

    launch_browser && open_in_default_browser(url)
    # wait until user interrupts the LiveServer (using CTRL+C).
    try
        counter = 1
        while true
            if WS_INTERRUPT.x || fw.status == :interrupted
                # rethrow the interruption (which may have happened during
                # the websocket handling or during the file watching)
                throw(InterruptException())
            end
            # do the auxiliary function if there is one (by default this does nothing)
            coreloopfun(counter, fw)
            # update the cycle counter and sleep (yields to other threads)
            counter += 1
            sleep(0.1)
        end
    catch err
        if !isa(err, InterruptException)
            throw(err)
        end
    finally
        # cleanup: close everything that might still be alive
        VERBOSE[] && println("\n⋮ shutting down LiveServer")
        # stop the filewatcher
        stop(fw)
        # close any remaining websockets
        for wss ∈ values(WS_VIEWERS), wsi ∈ wss
            close(wsi.io)
        end
        # empty the dictionary of viewers
        empty!(WS_VIEWERS)
        # shut down the server
        close(server)
        VERBOSE[] && println("\n✓ LiveServer shut down.")
        # reset environment variables
        CONTENT_DIR[] = ""
        WS_INTERRUPT[] = false
    end
    return nothing
end


text:
###################################################
## geo/stats/statistics.jl
## define methods to calculate errors on GeoTrips using DataSplits
###################################################

"""
	`tripsLogError`: compute log error on passed set of indices
	Can pass in precomputed times for faster runtime
"""
function tripsLogError(gt::GeoTimings, IDlist::Vector{Int}, times::Vector{Float64}=Float64[])
	error = 0.
	counter = 0
	precomputed = (length(times) == length(gt.trips))
	for ID in IDlist
		if precomputed
			timing = times[ID]
		else
			timing = estimateTime(gt, ID)
		end
		if timing > 0
			error += log(timing/gt.trips[ID].time)^2
			counter += 1
		end
	end
	error = error/counter
	return error
end
trTripsLogError(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[]) = tripsLogError(gt, trainSet(ds), times)
testTripsLogError(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[]) = tripsLogError(gt, testSet(ds), times)

"""
	`tripsLogBias`: compute log bias on trips for whatever set of indices is passed
	Log bias is defined as (log predicted time - log real time)
	Can pass in precomputed times for faster runtime
"""
function tripsLogBias(gt::GeoTimings, IDlist::Vector{Int}, times::Vector{Float64}=Float64[])
	bias = 0.
	counter = 0
	precomputed = (length(times) == length(gt.trips))
	for ID in IDlist
		if precomputed
			timing = times[ID]
		else
			timing = estimateTime(gt, ID)
		end
		if timing > 0
			bias += log(timing/gt.trips[ID].time)
			counter += 1
		end
	end
	bias = bias / counter
	return bias
end
testTripsLogBias(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[]) = tripsLogBias(gt, testSet(ds), times)
trTripsLogBias(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[]) = tripsLogBias(gt, trainSet(ds), times)

"""
	`tripsLogErrorByTime`: returns RMS on indicated subset of GeoData, subsetting errors in different time buckets
	Can pass in precomputed times for faster runtime
"""
function tripsLogErrorByTime(gt::GeoTimings, IDlist::Vector{Int}, times::Vector{Float64}=Float64[], timeBound::Array{Float64,1} = [270., 450., 720., 900., 100_000.])
	numBins = length(timeBound)
	error = 0. * collect(1:numBins)
	numInBin = 0 * collect(1:numBins)
	precomputed = (length(times) == length(gt.trips))
	for ID in IDlist
		idx = 1
		if precomputed
			timing = times[ID]
		else
			timing = estimateTime(gt, ID)
		end
		while timing > timeBound[idx] && idx < length(timeBound)
			idx += 1
		end
		if timing > 0
			error[idx] += log(timing/gt.trips[ID].time)^2
			numInBin[idx] += 1
		end
	end
	error = error ./ numInBin
	return error
end
trTripsLogErrorByTime(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[], timeBound::Array{Float64, 1} = [270., 450., 720., 900., 100_000.]) = tripsLogErrorByTime(gt, trainSet(ds), times, timeBound)
testTripsLogErrorByTime(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[], timeBound::Array{Float64, 1} = [270., 450., 720., 900., 100_000.]) = tripsLogErrorByTime(gt, testSet(ds), times, timeBound)

"""
	`tripsLogBiasByTime`: returns log bias on indicated subset of GeoData, subsetting errors in different time buckets
	Can pass in precomputed times for faster runtime
"""
function tripsLogBiasByTime(gt::GeoTimings, IDlist::Vector{Int}, times::Vector{Float64}=Float64[], timeBound::Array{Float64,1} = [270.,450., 720., 900., 100_000.])
	numBins = length(timeBound)
	error = 0. * collect(1:numBins)
	numInBin = 0 * collect(1:numBins)
	precomputed = (length(times) == length(gt.trips))
	for ID in IDlist
		idx = 1
		if precomputed
			timing = times[ID]
		else
			timing = estimateTime(gt, ID)
		end
		while timing > timeBound[idx] && idx < length(timeBound)
			idx += 1
		end
		if timing > 0
			error[idx] += log(timing/gt.trips[ID].time)
			numInBin[idx] += 1
		end
	end
	error = error ./ numInBin
	return error
end
trTripsLogBiasByTime(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[], timeBound::Array{Float64, 1} = [270., 450., 720., 900., 100_000.]) = tripsLogBiasByTime(gt, trainSet(ds), times, timeBound)
testTripsLogBiasByTime(gt::GeoTimings, ds::DataSplit, times::Vector{Float64}=Float64[], timeBound::Array{Float64, 1} = [270., 450., 720., 900., 100_000.]) = tripsLogBiasByTime(gt, testSet(ds), times, timeBound)


text:
#= test/permute.jl
=#

@testset "permute" begin
  n = 10
  p = randperm(n)

  @time for t in (RowPermutation, ColumnPermutation)
    println("  $t")

    A = t(p)
    P = Matrix(A)
    @test P == sparse(A)
    @test P == A*Matrix(I,n,n)
    @test P == Matrix(I,n,n)*A

    x = rand(n)
    @test A  *x == P  *x
    @test A' *x == P' *x
    @test transpose(A)*x == transpose(P)*x

    B = rand(n, n)
    @test A  *B   == P  *B
    @test A  *B'  == P  *B'
    @test A  *transpose(B) == P  *transpose(B)
    @test A' *B   == P' *B
    @test A' *B'  == P' *B'
    @test transpose(A)*B   == transpose(P)*B
    @test transpose(A)*transpose(B) == transpose(P)*transpose(B)
    @test B  *A   == B  *P
    @test B  *A'  == B  *P'
    @test B  *transpose(A) == B  *transpose(P)
    @test B' *A   == B' *P
    @test B' *A'  == B' *P'
    @test transpose(B)*A   == transpose(B)*P
    @test transpose(B)*transpose(A) == transpose(B)*transpose(P)
  end
end


text:
module Mat

using StaticArrays

function translation(::Type{T}, xyz::Real...) where {T <: Real}
    let x, y, z = xyz
        return @SMatrix T[
            1   0   0   x
            0   1   0   y
            0   0   1   z
            0   0   0   1
        ]
    end
end

translation(xyz::Real...) = translation(Float64, xyz...)

function scaling(::Type{T}, xyz::Real...) where {T <: Real}
    let x, y, z = xyz
        return @SMatrix T[
            x   0   0   0
            0   y   0   0
            0   0   z   0
            0   0   0   1
        ]
    end
end

scaling(xyz::Real...) = scaling(Float64, xyz...)

rotationX(ϕ::Real, T::DataType = Float64) = @SMatrix T[
    1   0        0       0
    0   cos(ϕ)  -sin(ϕ)  0
    0   sin(ϕ)   cos(ϕ)  0
    0   0        0       1
]

rotationY(ϕ::Real, T::DataType = Float64) = @SMatrix T[
    cos(ϕ)  0   sin(ϕ)  0
    0       1   0       0
   -sin(ϕ)  0   cos(ϕ)  0
    0       0   0       1
]

rotationZ(ϕ::Real, T::DataType = Float64) = @SMatrix T[
    cos(ϕ)  -sin(ϕ)  0  0
    sin(ϕ)   cos(ϕ)  0  0
    0        0       1  0
    0        0       0  1
]

function rotation(::Type{T}; x::Real = 0, y::Real = 0, z::Real = 0) where {T <: Real}
    (x == 0 ? 1 : rotationX(x, T)) *
    (y == 0 ? 1 : rotationY(y, T)) *
    (z == 0 ? 1 : rotationZ(z, T))
end

rotation(; xyz...) = rotation(Float64; xyz...)

end # Mat


text:
using ThompsonSampling, Base.Test

x = linspace(0,2pi,100)
y = linspace(0,2pi,100)

g(x,y) = sin(x)*sin(y)

gη(x) = g(x[1],x[2]) + 0.05*randn()

x0 = [π 3π/2;π π+1.0]
y0 = [gη(x0[:,i]) for i in 1:size(x0,2)]

t = ThompsonSampler(gη,x0,y0,100,[0.0,0.0],[2pi,2pi])

collect(Iterators.take(t,2))
@time collect(Iterators.take(t,50));