text:
import Base: getindex
import AutoGrad: recording

abstract type SlothLayer; end

struct Sequential
    layers::Array
end

function (model::Sequential)(x)
    for layer in model.layers
        x = layer(x)
    end
    return x
end

getindex(seq::Sequential, ind::Union{Int, UnitRange}) = seq.layers[ind]

struct Linear <: SlothLayer
    in_features::Int
    out_features::Int
    weight::SlothWeight
    bias::SlothBias
end

function (layer::Linear)(x)
    return layer.weight * x .+ layer.bias
end

function Linear(; in_features::Int, out_features::Int, init=xavier, bias=true)
    w = param(out_features, in_features; init=init)
    b = ifelse(bias, param0(out_features, 1), F(0.0))
    return Linear(in_features, out_features, w, b)
end

struct Conv <: SlothLayer
    in_channels::Int
    out_channels::Int
    kernel_size::IntHyperparam
    stride::IntHyperparam
    padding::IntHyperparam
    dilation::IntHyperparam
    mode::Int
    weight::SlothWeight
    bias::SlothBias
end

function (layer::Conv)(x)
    y = conv4(layer.weight, x;
        padding=layer.padding,
        stride=layer.stride,
        dilation=layer.dilation,
        mode=layer.mode)
    return y .+ layer.bias
end

function Conv(;
    in_channels::Int,
    out_channels::Int,
    kernel_size::IntHyperparam,
    stride::IntHyperparam=1,
    padding::IntHyperparam=0,
    dilation::IntHyperparam=1,
    mode::Int=0,
    bias::Bool=true,
    init=xavier)

    k = kernel_size
    k = ifelse(typeof(k) <: Int, (k, k), k)

    w = param(k..., in_channels, out_channels; init=init)
    b = ifelse(bias, param0(1, 1, out_channels, 1), F(0.0))

    return Conv(
        in_channels,
        out_channels,
        kernel_size,
        stride,
        padding,
        dilation,
        mode,
        w,
        b)
end

struct ConvTranspose <: SlothLayer
    in_channels::Int
    out_channels::Int
    kernel_size::IntHyperparam
    stride::IntHyperparam
    padding::IntHyperparam
    dilation::IntHyperparam
    mode::Int
    weight::SlothWeight
    bias::SlothBias
end

function (layer::ConvTranspose)(x)
    y = deconv4(layer.weight, x;
        padding=layer.padding,
        stride=layer.stride,
        dilation=layer.dilation,
        mode=layer.mode)
    y = y .+ layer.bias
end

function ConvTranspose(;
    in_channels::Int,
    out_channels::Int,
    kernel_size::IntHyperparam,
    stride::IntHyperparam=1,
    padding::IntHyperparam=0,
    dilation::IntHyperparam=1,
    mode::Int=0,
    bias::Bool=true,
    init=xavier)

    k = kernel_size
    k = ifelse(typeof(k) <: Int, (k, k), k)

    w = param(k..., out_channels, in_channels; init=init)
    b = ifelse(bias, param0(1, 1, out_channels, 1), F(0.0))

    return ConvTranspose(
        in_channels,
        out_channels,
        kernel_size,
        stride,
        padding,
        dilation,
        mode,
        w,
        b)
end

@enum PoolMode MAX_POOL=0 PADDED_AVG_POOL=1 AVG_POOL=2

struct Pool <: SlothLayer
    window::IntHyperparam
    stride::IntHyperparam
    padding::IntHyperparam
    dilation::IntHyperparam
    mode::PoolMode
end

function Pool(;
    window::IntHyperparam=2,
    stride::Union{Nothing, IntHyperparam}=nothing,
    padding::IntHyperparam=0,
    dilation::IntHyperparam=1,
    mode::PoolMode
)

    dilation != 1 && @warn "Knet does not support dilated pooling."
    stride = ifelse(isa(stride, Nothing), window, stride)
    return Pool(window, stride, padding, dilation, mode)
end

function (layer::Pool)(x)
    pool(x;
        window=layer.window,
        stride=layer.stride,
        padding=layer.padding,
        mode=Int(layer.mode))
end

MaxPool(; kwargs...) = Pool(; mode=MAX_POOL, kwargs...)
PaddedAvgPool(; kwargs...) = Pool(; mode=PADDED_AVG_POOL, kwargs...)
AvgPool(; kwargs...) = Pool(; mode=AVG_POOL, kwargs...)

struct BatchNorm{T<:AbstractFloat} <: SlothLayer
    num_features::Int
    momentum::T
    eps::T
    weight::SlothWeight
    moments::BNMoments
end

function (layer::BatchNorm)(x)
    batchnorm(x, layer.moments, layer.weight)
end

function BatchNorm(;
    num_features::Int,
    momentum::T=F(0.1),
    eps::T=F(1.0f-5)
) where T <: AbstractFloat
    weight, moments = Param(Knet.atype((bnparams(num_features)))), bnmoments()
    return BatchNorm(num_features, momentum, eps, weight, moments)
end

struct Embedding <: SlothLayer
    num_embeddings::Int
    embedding_dim::Int
    weight::SlothWeight
end

function (layer::Embedding)(x)
    return layer.weight[:, x]
end

function Embedding(; num_embeddings::Int, embedding_dim::Int, init=xavier)
    weight = param(embedding_dim, num_embeddings)
    return Embedding(num_embeddings, embedding_dim, weight)
end

struct Dropout <: SlothLayer
    p::AbstractFloat
    Dropout(; p::AbstractFloat=F(0.5)) = new(p)
end

function (layer::Dropout)(x; drop=recording())
    return dropout(x, layer.p; drop=drop)
end

struct Activation <: SlothLayer
    f
end

(l::Activation)(x) = l.f.(x)

Relu() = Activation(relu)
Tanh() = Activation(tanh)
Sigm() = Activation(sigm)

struct LeakyRelu
    α::AbstractFloat
    LeakyRelu(α::AbstractFloat=F(0.02)) = new(α)
end

function (layer::LeakyRelu)(x)
    pos = relu.(x)
    neg = min.(F(0.0), x)
    return pos + layer.α .* neg
end

struct Flatten <: SlothLayer; end

(l::Flatten)(x) = reshape(x, :,  size(x)[end])


text:
#####
##### Tests for basic manifold functionality
#####

types = [Vector{Float64},
         SizedVector{3, Float64},
         MVector{3, Float64},
         Vector{Float32},
         SizedVector{3, Float32},
         MVector{3, Float32}]

@testset "Euclidean" begin
    M = Euclidean(3)
    base_points = [randn(3) for _=1:3]
    t = rand()

    for T in types
        @testset "Type $T" begin
            x₀ = convert(T, base_points[1])
            v₀ = convert(T, base_points[2])
            p₀ = convert(T, base_points[3])

            x, p = cogeodesic_flow(M, x₀, p₀, t)
            @test x ≈ x₀ + t * p₀
            @test p ≈ p₀

            x, v = geodesic_flow(M, x₀, v₀, t)
            @test x ≈ x₀ + t * v₀
            @test v ≈ v₀
        end
    end
end

@testset "Sphere" begin
    M = Sphere(3)

    base_points = [normalize(randn(3))]
    Π = I - base_points[1] * base_points[1]'
    push!(base_points, normalize(Π * randn(3)))
    push!(base_points, normalize(Π * randn(3)))

    for T in types
        @testset "Type $T" begin
            x₀ = convert(T, base_points[1])
            v₀ = convert(T, base_points[2])
            p₀ = convert(T, base_points[3])
            t = eltype(T)(π)

            x, p = cogeodesic_flow(M, x₀, p₀, t)
            @test x ≈ -x₀
            @test p ≈ -p₀

            x, v = geodesic_flow(M, x₀, v₀, t)
            @test x ≈ -x₀
            @test v ≈ -v₀
        end
    end
end

ext: jl


text:
import Base.diag, Base.diagm
export diag, diagm

function diag(x::AffineExpr, num::Int64)
  if num <= -x.size[1] || num >= x.size[2]
    error("Out of bounds diagonal number for diag")
  end
  if num < 0
    start_row = -num + 1
    start_col = 1
  else
    start_row = 1
    start_col = num + 1
  end
  len = min(x.size[1] - start_row + 1, x.size[2] - start_col + 1)
  start_ind = (start_col - 1) * x.size[1] + start_row

  indexer = Constant(
    sparse(
      1 : len,
      start_ind : x.size[1] + 1 : start_ind + (len - 1) * (x.size[1] + 1),
      1.0,
      len,
      x.size[1] * x.size[2]
    )
  )

  vars_to_coeffs_map = Dict{UInt64, Constant}()
  for (v, c) in x.vars_to_coeffs_map
    vars_to_coeffs_map[v] = indexer * c
  end
  constant = indexer * x.constant
  this = AffineExpr(:diag, (x,), vars_to_coeffs_map, constant, (len, 1))
  this.evaluate = ()->diag(full(x.evaluate()), num)
  return this
end

diag(x::AffineExpr) = diag(x, 0)

function diagm(x::AffineExpr)
  if x.size[2] != 1
    error("Can only call diagm on column vectors")
  end

  indexer = Constant(
    sparse(
      1 : x.size[1] + 1 : 1 + (x.size[1] - 1) * (x.size[1] + 1),
      1 : x.size[1],
      1.0,
      x.size[1] * x.size[1],
      x.size[1]
    )
  )

  vars_to_coeffs_map = Dict{UInt64, Constant}()
  for (v, c) in x.vars_to_coeffs_map
    vars_to_coeffs_map[v] = indexer * c
  end
  constant = indexer * x.constant
  this = AffineExpr(:diagm, (x,), vars_to_coeffs_map, constant, (x.size[1], x.size[1]))
  this.evaluate = ()->diagm(vec(x.evaluate()))
  return this
end


text:
function sayhi(name)
    println("Hi $name, it's great to see you!")
end

function f(x)
    x^2
end

sayhi("C-3PO")

f(42)

sayhi2(name) = println("Hi $name, it's great to see you!")

f2(x) = x^2

sayhi2("R2D2")

f2(42)

sayhi3 = name -> println("Hi $name, it's great to see you!")

f3 = x -> x^2

sayhi3("Chewbacca")

f3(42)

sayhi(55595472)

A = rand(3, 3)
A

f(A)

f("hi")

v = rand(3)

f(v)

v = [3, 5, 2]

sort(v)

v

sort!(v)

v

A = [i + 3*j for j in 0:2, i in 1:3]

f(A)

B = f.(A)

A[2, 2]

A[2, 2]^2

B[2, 2]

v = [1, 2, 3]

f.(v)

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
# kminus,i,j = k*(j/Ns,i)
function kminusij(imax, Ns)

rows    = imax + 3
columns = Ns[imax] + 3
istart  = 3
idimer  = 4
iend    = rows - 1
j0      = 2
jmax    = columns - 1

kmij = zeros(Float64, rows, columns)

for i = idimer:iend
    monomers = i - 2
    for j = j0:Ns[monomers] + 2
        ligands = j - 2
        kmij[i,j] = (ligands)/Ns[monomers]
    end
end

return kmij

end


text:
export DenseMultiVector

"""
DenseMultiVector represents a dense multi-vector.  Note that all the vectors in a single DenseMultiVector are the same size
"""
mutable struct DenseMultiVector{Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer} <: MultiVector{Data, GID, PID, LID}
    data::Array{Data, 2} # data[1, 2] is the first element of the second vector
    numVectors::LID

    map::BlockMap{GID, PID, LID}
end

## Constructors ##

"""
    DenseMultiVector{Data}(::BlockMap{GID, PID, LID}, numVecs::Integer, zeroOut=true)

Creates a new DenseMultiVector based on the given map
"""
function DenseMultiVector{Data}(map::BlockMap{GID, PID, LID}, numVecs::Integer, zeroOut=true) where {Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    localLength = numMyElements(map)
    if zeroOut
        data = zeros(Data, (localLength, numVecs))
    else
        data = Array{Data, 2}(undef, localLength, numVecs)
    end
    DenseMultiVector{Data, GID, PID, LID}(data, numVecs, map)
end

"""
    DenseMultiVector(map::BlockMap{GID, PID, LID}, data::AbstractArray{Data, 2})

Creates a new DenseMultiVector wrapping the given array.  Changes to the DenseMultiVector or Array will affect the other
"""
function DenseMultiVector(map::BlockMap{GID, PID, LID}, data::AbstractArray{Data, 2}) where {Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    localLength = numMyElements(map)
    if size(data, 1) != localLength
        throw(InvalidArgumentError("Length of vectors does not match local length indicated by map"))
    end
    DenseMultiVector{Data, GID, PID, LID}(data, size(data, 2), map)
end

## External methods ##

numVectors(mVect::DenseMultiVector) = mVect.numVectors
getMap(mVect::DenseMultiVector) = mVect.map
getLocalArray(mVect::DenseMultiVector) = mVect.data

function Base.copy(vect::DenseMultiVector{Data, GID, PID, LID})::DenseMultiVector{Data, GID, PID, LID} where {Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    DenseMultiVector{Data, GID, PID, LID}(copy(vect.data), vect.numVectors, vect.map)
end

function Base.copyto!(dest::DenseMultiVector{Data, GID, PID, LID}, src::DenseMultiVector{Data, GID, PID, LID})::DenseMultiVector{Data, GID, PID, LID} where {Data, GID, PID, LID}
    copyto!(dest.data, src.data)
    dest.numVectors = src.numVectors
    dest.map = src.map

    dest
end


text:
# run our functions once so that they start precompiling
println("Warm-up start")

include("./main.jl")
@info "Warm-up: optimizing"
for tmax in 0.0:0.5:4.0
    # we go from infeasible to feasible, should trigger most of JuMP's relevant internals
    result = opt_controls_temp(;opt_parameters=Dict("temp_overshoot" => tmax, "temp_goal" => tmax))
    result = opt_controls_temp(;opt_parameters=Dict("temp_overshoot" => 999.0, "temp_goal" => tmax))
    MsgPack.pack(result)
end
forward_controls_temp()

@info "Warm-up: HTTP server"
t = @async run("127.0.0.1", 40404)

sleep(5)
download("http://127.0.0.1:40404/index.html")

# this is async because it blocks for some reason
@async Base.throwto(t, InterruptException())
sleep(2) # i am pulling these numbers out of thin air

println("Warm-up done")


text:
##############################################################################
##
## PooledDataArray type definition
##
## An AbstractDataArray with efficient storage when values are repeated. A
## PDA wraps an array of unsigned integers, which are used to index into a
## compressed pool of values. NA's are 0's in the refs array.  The unsigned
## integer type used for the refs array defaults to Uint32.  The `compact`
## function converts to a smallest integer size that will index the entire pool.
##
## TODO: Allow ordering of factor levels
## TODO: Add metadata for dummy conversion
##
##############################################################################

# This is used as a wrapper during PooledDataArray construction only, to distinguish
# arrays of pool indices from normal arrays
type RefArray{R<:Integer,N}
    a::Array{R,N}
end

type PooledDataArray{T, R<:Integer, N} <: AbstractDataArray{T, N}
    refs::Array{R, N}
    pool::Vector{T}

    function PooledDataArray(rs::RefArray{R, N},
                             p::Vector{T})
        # refs mustn't overflow pool
        if length(rs.a) > 0 && maximum(rs.a) > prod(size(p))
            error("Reference array points beyond the end of the pool")
        end
        new(rs.a,p)
    end
end
typealias PooledDataVector{T,R} PooledDataArray{T,R,1}
typealias PooledDataMatrix{T,R} PooledDataArray{T,R,2}

##############################################################################
##
## PooledDataArray constructors
##
# How do you construct a PooledDataArray from an Array?
# From the same sigs as a DataArray!
#
# Algorithm:
# * Start with:
#   * A null pool
#   * A pre-allocated refs
#   * A hash from T to Int
# * Iterate over d
#   * If value of d in pool already, set the refs accordingly
#   * If value is new, add it to the pool, then set refs
##############################################################################

# Echo inner constructor as an outer constructor
function PooledDataArray{T,R<:Integer,N}(refs::RefArray{R, N},
                                         pool::Vector{T})
    PooledDataArray{T,R,N}(refs, pool)
end

# A no-op constructor
PooledDataArray(d::PooledDataArray) = d

# Constructor from array, w/ pool, missingness, and ref type
function PooledDataArray{T,R<:Integer,N}(d::AbstractArray{T, N},
                                         pool::Vector{T},
                                         m::AbstractArray{Bool, N},
                                         r::Type{R} = DEFAULT_POOLED_REF_TYPE)
    if length(pool) > typemax(R)
        error("Cannot construct a PooledDataVector with type $R with a pool of size $(length(pool))")
    end

    newrefs = Array(R, size(d))
    poolref = Dict{T, R}()

    # loop through once to fill the poolref dict
    for i = 1:length(pool)
        poolref[pool[i]] = i
    end

    # fill in newrefs
    for i = 1:length(d)
        if m[i]
            newrefs[i] = 0
        else
            newrefs[i] = get(poolref, d[i], 0)
        end
    end
    return PooledDataArray(RefArray(newrefs), pool)
end

# Constructor from array, w/ missingness and ref type
function PooledDataArray{T,R<:Integer,N}(d::AbstractArray{T, N},
                                         m::AbstractArray{Bool, N},
                                         r::Type{R} = DEFAULT_POOLED_REF_TYPE)
    pool = sort(unique(d[!m]))
    PooledDataArray(d, pool, m, r)
end

# Construct an all-NA PooledDataVector of a specific type
PooledDataArray(t::Type, dims::Int...) = PooledDataArray(Array(t, dims...), trues(dims...))
PooledDataArray{R<:Integer}(t::Type, r::Type{R}, dims::Int...) = PooledDataArray(Array(t, dims...), trues(dims...), r)

# Construct an empty PooledDataVector of a specific type
PooledDataArray(t::Type) = PooledDataArray(similar(Array(t,1),0), trues(0))
PooledDataArray{R<:Integer}(t::Type, r::Type{R}) = PooledDataArray(similar(Array(t,1),0), trues(0), r)

# Convert a BitArray to an Array{Bool} (m = missingness)
# For some reason an additional method is needed but even that doesn't work
# For a BitArray a refs type of Uint8 will always be sufficient as the size of the pool is 0, 1 or 2
PooledDataArray{N}(d::BitArray{N}) = PooledDataArray(bitunpack(d), falses(size(d)), Uint8)
PooledDataArray{N}(d::BitArray{N}, m::AbstractArray{Bool, N}) = PooledDataArray(bitunpack(d), m, Uint8)

# Convert a DataArray to a PooledDataArray
PooledDataArray{T,R<:Integer}(da::DataArray{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(da.data, da.na, r)
PooledDataArray{T,R<:Integer}(da::DataArray{T},
                              pool::Vector{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(da.data, pool, da.na, r)

# Convert a Array{T} to a PooledDataArray
PooledDataArray{T,R<:Integer}(d::Array{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(d, falses(size(d)), r)
PooledDataArray{T,R<:Integer}(d::Array{T},
                              pool::Vector{T},
                              r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray(d, pool, falses(size(d)), r)

# Explicitly convert Ranges into a PooledDataVector
PooledDataArray{R<:Integer}(rs::Ranges,
                            r::Type{R} = DEFAULT_POOLED_REF_TYPE) = PooledDataArray([rs], falses(length(rs)), r)

# Initialized constructors with 0's, 1's
for (f, basef) in ((:pdatazeros, :zeros), (:pdataones, :ones))
    @eval begin
        ($f)(dims::Int...) = PooledDataArray(($basef)(dims...), falses(dims...))
        ($f)(t::Type, dims::Int...) = PooledDataArray(($basef)(t, dims...), falses(dims...))
        ($f){R<:Integer}(t::Type, r::Type{R}, dims::Int...) = PooledDataArray(($basef)(t, dims...), falses(dims...), r)
    end
end

##############################################################################
##
## Basic size properties of all Data* objects
##
##############################################################################

Base.size(pda::PooledDataArray) = size(pda.refs)
Base.length(pda::PooledDataArray) = length(pda.refs)
Base.endof(pda::PooledDataArray) = endof(pda.refs)

##############################################################################
##
## Copying Data* objects
##
##############################################################################

Base.copy(pda::PooledDataArray) = PooledDataArray(RefArray(copy(pda.refs)),
                                                  copy(pda.pool))
# TODO: Implement copy_to()

function Base.append!{T,R1,R2}(x::PooledDataArray{T,R1,1}, y::PooledDataArray{T,R2,1})
    pool = myunique(vcat(x.pool, y.pool))
    xn, yn = length(x), length(y)
    resize!(x.refs, xn + yn)

    xtidx::Array{R1} = findat(pool, x.pool)
    for i in 1:xn
        if x.refs[i] != 0
            x.refs[i] = xtidx[x.refs[i]]
        end
    end

    ytidx::Array{R1} = findat(pool, y.pool)
    for i in 1:yn
        x.refs[xn + i] = y.refs[i] == 0 ? 0 : ytidx[y.refs[i]]
    end

    x.pool = pool
    return x
end

function Base.append!{T,R}(x::PooledDataArray{T,R,1}, y::AbstractArray{T,1})
    poolref = Dict{T, R}()
    for i in 1:length(x.pool)
        poolref[x.pool[i]] = i
    end
    for v in y
        if !isna(v) && !haskey(poolref, v)
            poolref[v] = length(poolref) + 1
        end
    end

    pool = Array(T, length(poolref))
    for (v, i) in poolref
        pool[i] = v
    end

    xn, yn = length(x), length(y)
    resize!(x.refs, xn + yn)

    xtidx::Array{R} = findat(pool, x.pool)
    for i in 1:xn
        if x.refs[i] != 0
            x.refs[i] = xtidx[x.refs[i]]
        end
    end

    for i in 1:yn
        if !isna(y[i])
            x.refs[i+xn] = poolref[y[i]]
        end
    end
    x.pool = pool
    return x
end

##############################################################################
##
## Predicates, including the new isna()
##
##############################################################################

function Base.isnan(pda::PooledDataArray)
    PooledDataArray(RefArray(copy(pda.refs)), isnan(pda.pool))
end

function Base.isfinite(pda::PooledDataArray)
    PooledDataArray(RefArray(copy(pda.refs)), isfinite(pda.pool))
end

isna(pda::PooledDataArray) = pda.refs .== 0

#' @description
#'
#' Safe and type-stable way to determine if element `i` of an
#' PooledDataArray is `NA`.
#'
#' @param a::PooledDataArray The PooledDataArray whose missingness will
#'        be assessed.
#' @param i::Integer The index of the element to be checked for `NA`.
#'
#' @returns na::Bool Is the element `NA` or not?
#'
#' @examples
#'
#' a = @pdata([1, 2, 3])
#' isna(a, 1)
isna(pda::PooledDataArray, i::Real) = pda.refs[i] == 0 # -> Bool

##############################################################################
##
## PooledDataArray utilities
##
## TODO: Add methods with these names for DataArray's
##       Decide whether levels() or unique() is primitive. Make the other
##       an alias.
##  Tom: I don't think levels and unique are the same. R doesn't include NA's
##       with levels, but it does with unique. Having these different is
##       useful.
##
##############################################################################

function compact{T,R<:Integer,N}(d::PooledDataArray{T,R,N})
    sz = length(d.pool)

    REFTYPE = sz <= typemax(Uint8)  ? Uint8 :
              sz <= typemax(Uint16) ? Uint16 :
              sz <= typemax(Uint32) ? Uint32 :
                                      Uint64

    if REFTYPE == R
        return d
    end

    newrefs = convert(Array{REFTYPE, N}, d.refs)
    PooledDataArray(RefArray(newrefs), d.pool)
end

function Base.unique{T}(x::PooledDataArray{T})
    if anyna(x)
        n = length(x.pool)
        d = Array(T, n + 1)
        for i in 1:n
            d[i] = x.pool[i]
        end
        m = falses(n + 1)
        m[n + 1] = true
        return DataArray(d, m)
    else
        return DataArray(copy(x.pool), falses(length(x.pool)))
    end
end

levels{T}(pda::PooledDataArray{T}) = copy(pda.pool)

function PooledDataArray{S,R,N}(x::PooledDataArray{S,R,N},
                                newpool::Vector{S})
    # QUESTION: should we have a ! version of this? If so, needs renaming?
    tidx::Array{R} = findat(newpool, x.pool)
    refs = zeros(R, length(x))
    for i in 1:length(refs)
        if x.refs[i] != 0
            refs[i] = tidx[x.refs[i]]
        end
    end
    PooledDataArray(RefArray(refs), newpool)
end

myunique(x::AbstractVector) = x[sort(unique(findat(x, x)))]  # gets the ordering right
myunique(x::AbstractDataVector) = myunique(dropna(x))   # gets the ordering right; removes NAs

function setlevels{T,R}(x::PooledDataArray{T,R}, newpool::AbstractVector)
    pool = myunique(newpool)
    refs = zeros(R, length(x))
    tidx::Array{R} = findat(pool, newpool)
    tidx[isna(newpool)] = 0
    for i in 1:length(refs)
        if x.refs[i] != 0
            refs[i] = tidx[x.refs[i]]
        end
    end
    return PooledDataArray(RefArray(refs), pool)
end

function setlevels!{T,R}(x::PooledDataArray{T,R}, newpool::AbstractVector{T})
    if newpool == myunique(newpool) # no NAs or duplicates
        x.pool = newpool
        return x
    else
        x.pool = myunique(newpool)
        tidx::Array{R} = findat(x.pool, newpool)
        tidx[isna(newpool)] = 0
        for i in 1:length(x.refs)
            if x.refs[i] != 0
                x.refs[i] = tidx[x.refs[i]]
            end
        end
        return x
    end
end

function setlevels(x::PooledDataArray, d::Dict)
    newpool = copy(DataArray(x.pool))
    # An NA in `v` is put in the pool; that will cause it to become NA
    for (k,v) in d
        idx = findin(newpool, [k])
        if length(idx) == 1
            newpool[idx[1]] = v
        end
    end
    setlevels(x, newpool)
end

function setlevels!{T,R}(x::PooledDataArray{T,R}, d::Dict{T,T})
    for (k,v) in d
        idx = findin(x.pool, [k])
        if length(idx) == 1
            x.pool[idx[1]] = v
        end
    end
    x
end

function setlevels!{T,R}(x::PooledDataArray{T,R}, d::Dict{T,Any}) # this version handles NAs in d's values
    newpool = copy(DataArray(x.pool))
    # An NA in `v` is put in the pool; that will cause it to become NA
    for (k,v) in d
        idx = findin(newpool, [k])
        if length(idx) == 1
            newpool[idx[1]] = v
        end
    end
    setlevels!(x, newpool)
end

reorder(x::PooledDataArray) = PooledDataArray(x, sort(levels(x)))  # just re-sort the pool

reorder(x::PooledDataArray, y::AbstractVector...) = reorder(mean, x, y...)

### FIXME: this can't work because we don't know about DataFrames
# reorder(fun::Function, x::PooledDataArray, y::AbstractVector...) =
#     reorder(fun, x, DataFrame({y...}))

Base.reverse(x::PooledDataArray) = PooledDataArray(RefArray(reverse(x.refs)), x.pool)

##############################################################################
##
## similar()
##
##############################################################################

function Base.similar{T,R}(pda::PooledDataArray{T,R}, S::Type, dims::Dims)
    PooledDataArray(RefArray(zeros(R, dims)), S[])
end

##############################################################################
##
## find()
##
##############################################################################

Base.find(pdv::PooledDataVector{Bool}) = find(array(pdv, false))

##############################################################################
##
## getindex()
##
##############################################################################

# pda[SingleItemIndex]
function Base.getindex(pda::PooledDataArray, i::Real)
    if pda.refs[i] == 0
        return NA
    else
        return pda.pool[pda.refs[i]]
    end
end

# pda[MultiItemIndex]
function Base.getindex(pda::PooledDataArray, inds::AbstractDataVector{Bool})
    inds = find(inds)
    return PooledDataArray(RefArray(pda.refs[inds]), copy(pda.pool))
end
function Base.getindex(pda::PooledDataArray, inds::AbstractDataVector)
    inds = dropna(inds)
    return PooledDataArray(RefArray(pda.refs[inds]), copy(pda.pool))
end
function Base.getindex(pda::PooledDataArray, inds::Union(Vector, BitVector, Ranges))
    return PooledDataArray(RefArray(pda.refs[inds]), copy(pda.pool))
end

# pdm[SingleItemIndex, SingleItemIndex)
function Base.getindex(pda::PooledDataArray, i::Real, j::Real)
    if pda.refs[i, j] == 0
        return NA
    else
        return pda.pool[pda.refs[i, j]]
    end
end

# pda[SingleItemIndex, MultiItemIndex]
function Base.getindex(pda::PooledDataArray, i::Real, col_inds::AbstractDataVector{Bool})
    getindex(pda, i, find(col_inds))
end
function Base.getindex(pda::PooledDataArray, i::Real, col_inds::AbstractDataVector)
    getindex(pda, i, dropna(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
             i::Real,
             col_inds::Union(Vector, BitVector, Ranges))
    error("not yet implemented")
    PooledDataArray(RefArray(pda.refs[i, col_inds]), pda.pool[i, col_inds])
end

# pda[MultiItemIndex, SingleItemIndex]
function Base.getindex(pda::PooledDataArray, row_inds::AbstractDataVector{Bool}, j::Real)
    getindex(pda, find(row_inds), j)
end
function Base.getindex(pda::PooledDataArray, row_inds::AbstractVector, j::Real)
    getindex(pda, dropna(row_inds), j)
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
             row_inds::Union(Vector, BitVector, Ranges),
             j::Real)
    error("not yet implemented")
    PooledDataArray(RefArray(pda.refs[row_inds, j]), pda.pool[row_inds, j])
end

# pda[MultiItemIndex, MultiItemIndex]
function Base.getindex(pda::PooledDataArray,
             row_inds::AbstractDataVector{Bool},
             col_inds::AbstractDataVector{Bool})
    getindex(pda, find(row_inds), find(col_inds))
end
function Base.getindex(pda::PooledDataArray,
             row_inds::AbstractDataVector{Bool},
             col_inds::AbstractDataVector)
    getindex(pda, find(row_inds), dropna(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
             row_inds::AbstractDataVector{Bool},
             col_inds::Union(Vector, BitVector, Ranges))
    getindex(pda, find(row_inds), col_inds)
end
function Base.getindex(pda::PooledDataArray,
             row_inds::AbstractDataVector,
             col_inds::AbstractDataVector{Bool})
    getindex(pda, dropna(row_inds), find(col_inds))
end
function Base.getindex(pda::PooledDataArray,
             row_inds::AbstractDataVector,
             col_inds::AbstractDataVector)
    getindex(pda, dropna(row_inds), dropna(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
             row_inds::AbstractDataVector,
             col_inds::Union(Vector, BitVector, Ranges))
    getindex(pda, dropna(row_inds), col_inds)
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
             row_inds::Union(Vector, BitVector, Ranges),
             col_inds::AbstractDataVector{Bool})
    getindex(pda, row_inds, find(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
             row_inds::Union(Vector, BitVector, Ranges),
             col_inds::AbstractDataVector)
    getindex(pda, row_inds, dropna(col_inds))
end
# TODO: Make inds::AbstractVector
function Base.getindex(pda::PooledDataArray,
                  row_inds::Union(Vector, BitVector, Ranges),
                  col_inds::Union(Vector, BitVector, Ranges))
    error("not yet implemented")
    PooledDataArray(RefArray(pda.refs[row_inds, col_inds]), pda.pool[row_inds, col_inds])
end

##############################################################################
##
## setindex!() definitions
##
##############################################################################

function getpoolidx{T,R<:Union(Uint8, Uint16, Int8, Int16)}(pda::PooledDataArray{T,R}, val::Any)
    val::T = convert(T,val)
    pool_idx = findfirst(pda.pool, val)
    if pool_idx <= 0
        push!(pda.pool, val)
        pool_idx = length(pda.pool)
        if pool_idx > typemax(R)
            error("You're using a PooledDataArray with ref type $R, which can only hold $(int(typemax(R))) values,\n",
                  "and you just tried to add the $(typemax(R)+1)th reference.  Please change the ref type\n",
                  "to a larger int type, or use the default ref type ($DEFAULT_POOLED_REF_TYPE).")
        end
    end
    return pool_idx
end

function getpoolidx{T,R}(pda::PooledDataArray{T,R}, val::Any)
    val::T = convert(T,val)
    pool_idx = findfirst(pda.pool, val)
    if pool_idx <= 0
        push!(pda.pool, val)
        pool_idx = length(pda.pool)
    end
    return pool_idx
end

# x[SingleIndex] = NA
# TODO: Delete values from pool that no longer exist?
# Not a good idea.  Add another function called drop_unused_levels to do this.
# R has the convention that if f is a factor then factor(f) drops unused levels
function Base.setindex!(x::PooledDataArray, val::NAtype, ind::Real)
    x.refs[ind] = 0
    return x
end

# x[SingleIndex] = Single Item
# TODO: Delete values from pool that no longer exist?
function Base.setindex!{T,R}(x::PooledDataArray{T,R}, val::Any, ind::Real)
    val = convert(T, val)
    x.refs[ind] = getpoolidx(x, val)
    return x
end

# x[MultiIndex] = NA
# TODO: Find a way to delete the next four methods
function Base.setindex!(x::PooledDataArray{NAtype},
                val::NAtype,
                inds::AbstractVector{Bool})
    error("Don't use PooledDataVector{NAtype}'s")
end
function Base.setindex!(x::PooledDataArray{NAtype},
                val::NAtype,
                inds::AbstractVector)
    error("Don't use PooledDataVector{NAtype}'s")
end
function Base.setindex!(x::PooledDataArray, val::NAtype, inds::AbstractVector{Bool})
    inds = find(inds)
    x.refs[inds] = 0
    return x
end
function Base.setindex!(x::PooledDataArray, val::NAtype, inds::AbstractVector)
    x.refs[inds] = 0
    return x
end

# pda[MultiIndex] = Multiple Values
function Base.setindex!(pda::PooledDataArray,
                   vals::AbstractVector,
                   inds::AbstractVector{Bool})
    setindex!(pda, vals, find(inds))
end
function Base.setindex!(pda::PooledDataArray,
                   vals::AbstractVector,
                   inds::AbstractVector)
    for (val, ind) in zip(vals, inds)
        pda[ind] = val
    end
    return pda
end

# pda[SingleItemIndex, SingleItemIndex] = NA
function Base.setindex!{T,R}(pda::PooledDataMatrix{T,R}, val::NAtype, i::Real, j::Real)
    pda.refs[i, j] = zero(R)
    return pda
end
# pda[SingleItemIndex, SingleItemIndex] = Single Item
function Base.setindex!{T,R}(pda::PooledDataMatrix{T,R}, val::Any, i::Real, j::Real)
    val = convert(T, val)
    pda.refs[i, j] = getpoolidx(pda, val)
    return pda
end

##############################################################################
##
## show() and similar methods
##
##############################################################################

function Base.string(x::PooledDataVector)
    tmp = join(x, ", ")
    return "[$tmp]"
end

# # Need setindex!()'s to make this work
# This is broken now because the inner show returns to the outer show.
# function show(io::IO, pda::PooledDataArray)
#     invoke(show, (Any, AbstractArray), io, pda)
#     print(io, "\nlevels: ")
#     print(io, levels(pda))
# end

##############################################################################
##
## Replacement operations
##
##############################################################################

function replace!(x::PooledDataArray{NAtype}, fromval::NAtype, toval::NAtype)
    NA # no-op to deal with warning
end
function replace!(x::PooledDataArray, fromval::NAtype, toval::NAtype)
    NA # no-op to deal with warning
end
function replace!{S, T}(x::PooledDataArray{S}, fromval::T, toval::NAtype)
    fromidx = findfirst(x.pool, fromval)
    if fromidx == 0
        error("can't replace a value not in the pool in a PooledDataVector!")
    end

    x.refs[x.refs .== fromidx] = 0

    return NA
end
function replace!{S, T}(x::PooledDataArray{S}, fromval::NAtype, toval::T)
    toidx = findfirst(x.pool, toval)
    # if toval is in the pool, just do the assignment
    if toidx != 0
        x.refs[x.refs .== 0] = toidx
    else
        # otherwise, toval is new, add it to the pool
        push!(x.pool, toval)
        x.refs[x.refs .== 0] = length(x.pool)
    end

    return toval
end
function replace!{R, S, T}(x::PooledDataArray{R}, fromval::S, toval::T)
    # throw error if fromval isn't in the pool
    fromidx = findfirst(x.pool, fromval)
    if fromidx == 0
        error("can't replace a value not in the pool in a PooledDataArray!")
    end

    # if toval is in the pool too, use that and remove fromval from the pool
    toidx = findfirst(x.pool, toval)
    if toidx != 0
        x.refs[x.refs .== fromidx] = toidx
        #x.pool[fromidx] = None    TODO: what to do here??
    else
        # otherwise, toval is new, swap it in
        x.pool[fromidx] = toval
    end

    return toval
end

##############################################################################
##
## Sorting can use the pool to speed things up
##
##############################################################################

# TODO handle sortperm for non-sorted keys
Base.sortperm(pda::PooledDataArray) = groupsort_indexer(pda)[1]
function Base.sortperm(pda::PooledDataArray)
    if issorted(pda.pool)
        return groupsort_indexer(pda)[1]
    else
        return sortperm(reorder!(copy(pda)))
    end
end

Base.sortperm(pda::PooledDataArray, ::Base.Sort.ReverseOrdering) = reverse(sortperm(pda))
Base.sort(pda::PooledDataArray) = pda[sortperm(pda)]
Base.sort(pda::PooledDataArray, ::Base.Sort.ReverseOrdering) = pda[reverse(sortperm(pda))]
type FastPerm{O<:Base.Sort.Ordering,V<:AbstractVector} <: Base.Sort.Ordering
    ord::O
    vec::V
end
FastPerm{O<:Base.Sort.Ordering,V<:AbstractVector}(o::O,v::V) = FastPerm{O,V}(o,v)
Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ForwardOrdering,V}) = x[sortperm(o.vec)]
Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ReverseOrdering,V}) = x[reverse(sortperm(o.vec))]
Perm{O<:Base.Sort.Ordering}(o::O, v::PooledDataVector) = FastPerm(o, v)

##############################################################################
##
## PooledDataVecs: EXPLANATION SHOULD GO HERE
##
##############################################################################

function PooledDataVecs{S,Q<:Integer,R<:Integer,N}(v1::PooledDataArray{S,Q,N},
                                                   v2::PooledDataArray{S,R,N})
    pool = sort(unique([v1.pool, v2.pool]))
    sz = length(pool)

    REFTYPE = sz <= typemax(Uint8)  ? Uint8 :
              sz <= typemax(Uint16) ? Uint16 :
              sz <= typemax(Uint32) ? Uint32 :
                                      Uint64

    tidx1 = convert(Vector{REFTYPE}, findat(pool, v1.pool))
    tidx2 = convert(Vector{REFTYPE}, findat(pool, v2.pool))
    refs1 = zeros(REFTYPE, length(v1))
    refs2 = zeros(REFTYPE, length(v2))
    for i in 1:length(refs1)
        if v1.refs[i] != 0
            refs1[i] = tidx1[v1.refs[i]]
        end
    end
    for i in 1:length(refs2)
        if v2.refs[i] != 0
            refs2[i] = tidx2[v2.refs[i]]
        end
    end
    return (PooledDataArray(RefArray(refs1), pool),
            PooledDataArray(RefArray(refs2), pool))
end

function PooledDataVecs{S,R<:Integer,N}(v1::PooledDataArray{S,R,N},
                                        v2::AbstractArray{S,N})
    return PooledDataVecs(v1,
                          PooledDataArray(v2))
end
function PooledDataVecs{S,R<:Integer,N}(v1::PooledDataArray{S,R,N},
                                        v2::AbstractArray{S,N})
    return PooledDataVecs(v1,
                          PooledDataArray(v2))
end

####
function PooledDataVecs{S,R<:Integer,N}(v1::AbstractArray{S,N},
                                        v2::PooledDataArray{S,R,N})
    return PooledDataVecs(PooledDataArray(v1),
                          v2)
end

function PooledDataVecs(v1::AbstractArray,
                        v2::AbstractArray)

    ## Return two PooledDataVecs that share the same pool.

    ## TODO: allow specification of REFTYPE
    refs1 = Array(DEFAULT_POOLED_REF_TYPE, size(v1))
    refs2 = Array(DEFAULT_POOLED_REF_TYPE, size(v2))
    poolref = Dict{promote_type(eltype(v1), eltype(v2)), DEFAULT_POOLED_REF_TYPE}()
    maxref = 0

    # loop through once to fill the poolref dict
    for i = 1:length(v1)
        if !isna(v1[i])
            poolref[v1[i]] = 0
        end
    end
    for i = 1:length(v2)
        if !isna(v2[i])
            poolref[v2[i]] = 0
        end
    end

    # fill positions in poolref
    pool = sort(collect(keys(poolref)))
    i = 1
    for p in pool
        poolref[p] = i
        i += 1
    end

    # fill in newrefs
    zeroval = zero(DEFAULT_POOLED_REF_TYPE)
    for i = 1:length(v1)
        if isna(v1[i])
            refs1[i] = zeroval
        else
            refs1[i] = poolref[v1[i]]
        end
    end
    for i = 1:length(v2)
        if isna(v2[i])
            refs2[i] = zeroval
        else
            refs2[i] = poolref[v2[i]]
        end
    end

    return (PooledDataArray(RefArray(refs1), pool),
            PooledDataArray(RefArray(refs2), pool))
end

Base.convert{S,T,R1<:Integer,R2<:Integer,N}(::Type{PooledDataArray{S,R1,N}}, pda::PooledDataArray{T,R2,N}) =
    PooledDataArray(RefArray(convert(Array{R1,N}, pda.refs)), convert(Vector{S}, pda.pool))
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R,N}}, pda::PooledDataArray{T,R,N}) =
    PooledDataArray(RefArray(copy(pda.refs)), convert(Vector{S}, pda.pool))
Base.convert{T,R<:Integer,N}(::Type{PooledDataArray{T,R,N}}, pda::PooledDataArray{T,R,N}) = pda
Base.convert{S,T,R1<:Integer,R2<:Integer,N}(::Type{PooledDataArray{S,R1}}, pda::PooledDataArray{T,R2,N}) =
    convert(PooledDataArray{S,R1,N}, pda)
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S}}, pda::PooledDataArray{T,R,N}) =
    convert(PooledDataArray{S,R,N}, pda)
Base.convert{T,R<:Integer,N}(::Type{PooledDataArray}, pda::PooledDataArray{T,R,N}) = pda

Base.convert{T,R<:Integer,N}(::Type{PooledDataArray{T,R,N}}, a::DataArray{T,N}) =
    PooledDataArray(a.data, a.na, R)
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R,N}}, a::DataArray{T,N}) =
    convert(PooledDataArray{S,R,N}, PooledDataArray(a.data, a.na, R))
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R}}, a::DataArray{T,N}) =
    convert(PooledDataArray{S,R,N}, PooledDataArray(a.data, a.na, R))
Base.convert{S,T,N}(::Type{PooledDataArray{S}}, a::DataArray{T,N}) =
    convert(PooledDataArray{S}, PooledDataArray(a.data, a.na))
Base.convert(::Type{PooledDataArray}, a::DataArray) =
    PooledDataArray(a.data, a.na)

Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R,N}}, a::AbstractArray{T,N}) =
    PooledDataArray(convert(Array{S,N}, a), falses(size(a)), R)
Base.convert{S,T,R<:Integer,N}(::Type{PooledDataArray{S,R}}, a::AbstractArray{T,N}) =
    PooledDataArray(convert(Array{S,N}, a), falses(size(a)), R)
Base.convert{S,T,N}(::Type{PooledDataArray{S}}, a::AbstractArray{T,N}) =
    PooledDataArray(convert(Array{S,N}, a), falses(size(a)))
Base.convert(::Type{PooledDataArray}, a::AbstractArray) =
    PooledDataArray(a, falses(size(a)))

function Base.convert{S,T,R<:Integer,N}(::Type{DataArray{S,N}},
                                        pda::PooledDataArray{T,R,N})
    res = DataArray(Array(S, size(pda)), BitArray(size(pda)))
    for i in 1:length(pda)
        r = pda.refs[i]
        if r == 0 # TODO: Use zero(R)
            res.na[i] = true
        else
            res.na[i] = false
            res.data[i] = pda.pool[r]
        end
    end
    return res
end
Base.convert{S,T,R<:Integer,N}(::Type{DataArray{S}}, pda::PooledDataArray{T,R,N}) =
    convert(DataArray{T,N}, pda)
Base.convert{T,R<:Integer,N}(::Type{DataArray}, pda::PooledDataArray{T,R,N}) =
    convert(DataArray{T,N}, pda)

pdata(a::AbstractArray) = convert(PooledDataArray, a)

# Turn a PooledDataArray into an Array. Fail on NA
function array{T, R}(da::PooledDataArray{T, R})
    n = length(da)
    res = Array(T, size(da))
    for i in 1:n
        if da.refs[i] == zero(R)
            error(NAException())
        else
            res[i] = da.pool[da.refs[i]]
        end
    end
    return res
end

function array{T, R}(da::PooledDataArray{T, R}, replacement::T)
    n = length(da)
    res = Array(T, size(da))
    for i in 1:n
        if da.refs[i] == zero(R)
            res[i] = replacement
        else
            res[i] = da.pool[da.refs[i]]
        end
    end
    return res
end

function dropna{T}(pdv::PooledDataVector{T})
    n = length(pdv)
    res = Array(T, n)
    total = 0
    for i in 1:n
        if pdv.refs[i] > 0
            total += 1
            res[total] = pdv.pool[pdv.refs[i]]
        end
    end
    return res[1:total]
end


text:
@with_kw struct ModelParams
	n_agents = 10
	discretize_m = 100.
	jump_infect=true
end

@with_kw mutable struct Agent
    id::Int
	last_node::Int
	current_node::Int
	infected = false
	infection_time = -1
	total_route_len = 0.0
end

@with_kw struct VertexMeta
	ix::Int
	enu::ENU
	mapix::Union{Nothing,Int} = nothing
	osmid::Union{Nothing,Int} = nothing
	randomNode= (osmid != nothing)
	agents = Set{Int}()
end

@with_kw struct Simulation
    p::ModelParams
	g::SimpleGraph
	w::SparseMatrixCSC{Float64,Int64}  #edge weights
	vertex_meta::Vector{VertexMeta}
	agents::Vector{Agent}
	m::Union{OpenStreetMapX.MapData, Nothing}
	map_g_connected_points::Vector{Int}
	occupied_vertices::Set{Int}
	infected_agents::Set{Int}
	step_infected_agents_count::Vector{Int}
	step_total_avg_m_driven::Vector{Float64}
		#total number of meters driven by an average car
end

function Simulation(p::ModelParams, m::MapData; store_map::Bool=true)
	map_g_connected_points = sort(LightGraphs.strongly_connected_components(m.g),
	                    lt=(x,y)->length(x)<length(y), rev=true)[1]
	sort!(map_g_connected_points)
	map_g_connected_points_set = Set(map_g_connected_points)
	g = SimpleGraph()
	LightGraphs.add_vertices!(g,length(map_g_connected_points))
	vertex_meta = Vector{VertexMeta}(undef,length(map_g_connected_points))
	for i in 1:length(vertex_meta)
	    vertex_meta[i] = VertexMeta(;ix=i,
	        enu=m.nodes[m.n[map_g_connected_points[i]]],
	        mapix=map_g_connected_points[i],
	        osmid=m.n[map_g_connected_points[i]])
	end
	visited = Set{Pair{Int,Int}}() #used to not inlude bidirect edges

	wd = Dict{Tuple{Int,Int},Float64}() #edge weights as dict

	for e in edges(m.g)
	    (!(e.src in map_g_connected_points_set) || !(e.dst in map_g_connected_points_set)) && continue;
	    edgepair = Pair{Int,Int}(sort!([e.src, e.dst])...)
	    edgepair in visited && continue;
	    push!(visited, edgepair)
	    srcix = searchsortedfirst(map_g_connected_points,e.src)
	    dstix = searchsortedfirst(map_g_connected_points,e.dst)
	    d = distance(vertex_meta[srcix].enu,vertex_meta[dstix].enu)
	    add_vxs = max(0 , Int(ceil(d/p.discretize_m)-1) )
		piece_d = d/(add_vxs+1)
	    local lastix = srcix
	    for piece in 1:add_vxs
	        LightGraphs.add_vertex!(g)
	        local ix = nv(g)
	        LightGraphs.add_edge!(g, lastix, ix)
			wd[(lastix, ix)] = piece_d
	        lastix = ix
	        push!(vertex_meta,
	            VertexMeta(;ix=ix,
	                enu=enu_weighted(vertex_meta[srcix].enu,
	                               vertex_meta[dstix].enu,
	                               piece/(add_vxs+1) ))
	            )
	    end
	    LightGraphs.add_edge!(g, lastix, dstix)
		wd[(lastix, dstix)] = piece_d
	end
	w = SparseArrays.spzeros(length(vertex_meta),length(vertex_meta))

	for k in keys(wd)
		w[k[1], k[2]] = w[k[2], k[1]] = wd[k]
	end

	Simulation(	p=p, g=g, w=w, vertex_meta=vertex_meta,
	 			agents=Vector{Agent}(undef, p.n_agents),
				m=( store_map ? m : nothing ),
				map_g_connected_points=map_g_connected_points,
				occupied_vertices = Set{Int}(),
				infected_agents = Set{Int}(),
				step_infected_agents_count = Vector{Int}(),
				step_total_avg_m_driven = Vector{Float64}() )
end

function init!(s::Simulation; spawn_point::AbstractVector{Int} = 1:nv(s.g))
	empty!(s.occupied_vertices)
	empty!(s.infected_agents)
	empty!(s.step_infected_agents_count)
	empty!(s.step_total_avg_m_driven)
	for vm in s.vertex_meta
		empty!(vm.agents)
	end
	for i in 1:s.p.n_agents
		a = Agent(id=i,current_node=rand(spawn_point),last_node=-1)
		s.agents[i] = a
		push!(s.vertex_meta[a.current_node].agents, a.id)
		push!(s.occupied_vertices,a.current_node)
	end
	for aid in s.vertex_meta[s.agents[end].current_node].agents
		s.agents[aid].infected = true
		s.agents[aid].infection_time = 1
		push!(s.infected_agents, aid)
	end
	push!(s.step_infected_agents_count, length(s.infected_agents))
	push!(s.step_total_avg_m_driven, 0)
end

function enu_weighted(enu1::ENU, enu2::ENU, weight_perc::Float64)::ENU
	return ENU(enu1.east*(1-weight_perc) 	+ enu2.east*weight_perc,
		enu1.north*(1-weight_perc) 	+ enu2.north*weight_perc,
		enu1.up*(1-weight_perc) 	+ enu2.up*weight_perc)
end

function step!(s::Simulation)
	infected_nodes =  Set{Int}()
	from_infected_node_targets = Dict{Int, Set{Int}}()

	for node in s.occupied_vertices
		empty!(s.vertex_meta[node].agents)
	end
	empty!(s.occupied_vertices)
	for agent in s.agents
		vm = s.vertex_meta[agent.current_node]
		local nextnode::Int
		ns = LightGraphs.neighbors(s.g,agent.current_node)
		if vm.randomNode
			nextnode = rand(ns)
		else
			if agent.last_node!=ns[1]
				nextnode = ns[1]
			else
				nextnode = ns[2]
			end
		end

		agent.last_node = agent.current_node
		agent.current_node = nextnode
		agent.total_route_len += s.w[agent.last_node, agent.current_node]
		if agent.infected
			push!(infected_nodes, agent.current_node)
			if (s.p.jump_infect)
				push!(get!(from_infected_node_targets, agent.last_node, Set{Int}()), agent.current_node)
			end
		end
		push!(s.occupied_vertices, agent.current_node)
		push!(s.vertex_meta[agent.current_node].agents, agent.id)
	end
	if (s.p.jump_infect)
		for agent in s.agents
			if (!agent.infected) &&
					haskey(from_infected_node_targets, agent.current_node)
				# check if the agent just jumped over another infected agents
				if agent.last_node in from_infected_node_targets[agent.current_node]
					# the agent and all oter at the new node will get infected
					push!(infected_nodes, agent.current_node)
				end
			end
		end
	end

	step_no = length(s.step_infected_agents_count)+1
	for node in infected_nodes
		for aid in s.vertex_meta[node].agents
			if !s.agents[aid].infected
				s.agents[aid].infected = true
				s.agents[aid].infection_time = step_no
				push!(s.infected_agents, aid)
			end
		end
	end
	push!(s.step_infected_agents_count, length(s.infected_agents))
	push!(s.step_total_avg_m_driven,
		sum(a.total_route_len for a in s.agents)/length(s.agents))
end

function simulate_steps!(s::Simulation; max_idle=500)
	steps_without_change = 0
	last_infected_count = 0
	while (infected_count = length(s.infected_agents)) < length(s.agents)
		step!(s)
		if last_infected_count == infected_count
			steps_without_change += 1
			steps_without_change >= max_idle && break;
		else
			steps_without_change = 0
			last_infected_count = infected_count
		end
	end
end


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

mutable struct ApplicationGatewayCustomError <: SwaggerModel
    statusCode::Any # spec type: Union{ Nothing, String } # spec name: statusCode
    customErrorPageUrl::Any # spec type: Union{ Nothing, String } # spec name: customErrorPageUrl

    function ApplicationGatewayCustomError(;statusCode=nothing, customErrorPageUrl=nothing)
        o = new()
        validate_property(ApplicationGatewayCustomError, Symbol("statusCode"), statusCode)
        setfield!(o, Symbol("statusCode"), statusCode)
        validate_property(ApplicationGatewayCustomError, Symbol("customErrorPageUrl"), customErrorPageUrl)
        setfield!(o, Symbol("customErrorPageUrl"), customErrorPageUrl)
        o
    end
end # type ApplicationGatewayCustomError

const _property_map_ApplicationGatewayCustomError = Dict{Symbol,Symbol}(Symbol("statusCode")=>Symbol("statusCode"), Symbol("customErrorPageUrl")=>Symbol("customErrorPageUrl"))
const _property_types_ApplicationGatewayCustomError = Dict{Symbol,String}(Symbol("statusCode")=>"String", Symbol("customErrorPageUrl")=>"String")
Base.propertynames(::Type{ ApplicationGatewayCustomError }) = collect(keys(_property_map_ApplicationGatewayCustomError))
Swagger.property_type(::Type{ ApplicationGatewayCustomError }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayCustomError[name]))}
Swagger.field_name(::Type{ ApplicationGatewayCustomError }, property_name::Symbol) =  _property_map_ApplicationGatewayCustomError[property_name]

const _allowed_ApplicationGatewayCustomError_statusCode = ["HttpStatus403", "HttpStatus502"]

function check_required(o::ApplicationGatewayCustomError)
    true
end

function validate_property(::Type{ ApplicationGatewayCustomError }, name::Symbol, val)
    if name === Symbol("statusCode")
        Swagger.validate_param(name, "ApplicationGatewayCustomError", :enum, val, _allowed_ApplicationGatewayCustomError_statusCode)
    end
end


text:
include("generator_utils.jl")

# const repl_emoji_symbols_170
# const repl_emoji_symbols_below_170
include("repl_emoji_symbols.jl")

emoji_data     = get_emoji_data()
dict_170       = get_diff_dict(repl_emoji_symbols_170,       emoji_data)
dict_below_170 = get_diff_dict(repl_emoji_symbols_below_170, emoji_data)

# const repl_latex_symbols
include("repl_latex_symbols.jl")

gen_doc("../docs/src/additional_symbols_170.md",
        """# VERSION >= v\"1.7.0-DEV.849\"""",
        dict_170)
gen_doc("../docs/src/additional_symbols_below_170.md",
        """# VERSION <  v\"1.7.0-DEV.849\"""",
        dict_below_170)
gen_doc("../docs/src/emoji_symbols_in_repl_170.md",
        """# VERSION >= v\"1.7.0-DEV.849\"""",
        repl_emoji_symbols_170)
gen_doc("../docs/src/emoji_symbols_in_repl_below_170.md",
        """# VERSION <  v\"1.7.0-DEV.849\"""",
        repl_emoji_symbols_below_170)
gen_doc("../docs/src/latex_symbols_in_repl.md",
        """# LaTeX symbols in REPL""",
        repl_latex_symbols)


text:
"""
    AbstractMagnetismModel

Abstract type for various magnetism models.

Subtypes
* `Ising`
* `Potts`
* `XY`
"""
abstract type AbstractMagnetismModel end

function show(io::IO, model::AbstractMagnetismModel)
    println(typeof(model))
    show(model.params)
end

"""
    Ising(
        L         ::Int,
        β         ::Real;
        n_sweeps  ::Int    =10^5,
        cutoff    ::Float64=0.2,
        start_type::Symbol =:cold,
        seed      ::Int    =8,
        )

Type representing a 2D Ising model.

Arguments
* `L`: The side length of the lattice, i.e. for a square lattice there will be `L^2` spins
* `β`: `1 / kT`, i.e. the inverse temperature

Keyword Arguments
* `n_sweeps`  : The number of update sweeps to run over the lattice
* `cutoff`    : The percentage of measurements to throw away, i.e. throw away values during
    thermalization
* `start_type`: `:cold` starts the lattice with all spins aligned, `:hot` all spins random
* `seed`      : Seed value for the random number generator

Returns
* A new instance of `Ising`

Attributes
* `σ`          : Lattice containing the spins, values are ±1
* `params`     : Parameters type, contains relevant parameters for the simulation
* `observables`: Observables type, contains relevant observables data for the simulation
* `rng`        : Random number generator

# Examples
```julia-repl
julia> Ising(32, log(1 + √2) / 2)
Ising
L          = 32
β          = 0.3
n_sweeps   = 100000
start_type = cold
seed       = 8
q          = 0
```
"""
struct Ising <: AbstractMagnetismModel
    σ::Array{Int8,2}
    params::Parameters
    observables::Observables
    rng::AbstractRNG

    function Ising(
        L::Int,
        β::Real;
        n_sweeps::Int = 10^5,
        cutoff::Float64 = 0.2,
        start_type::Symbol = :cold,
        seed::Int = 8,
    )

        params = Parameters(L, β, n_sweeps, cutoff, start_type, seed, 0)
        rng = MersenneTwister(params.seed)
        σ = initialize_σ_Ising(params.L, params.start_type, rng)
        observables = Observables(Int)

        return new(σ, params, observables, rng)
    end
end

"""
    Potts(
        q         ::Int,
        L         ::Int,
        β         ::Real;
        n_sweeps  ::Int    =10^5,
        cutoff    ::Float64=0.2,
        start_type::Symbol =:cold,
        seed      ::Int    =8,
        )

Type representing a 2D Potts q-state model.

Arguments
* `q`: The number of different states each spin can take on
* `L`: The side length of the lattice, i.e. for a square lattice there will be `L^2` spins
* `β`: `1 / kT`, i.e. the inverse temperature

Keyword Arguments
* `n_sweeps`  : The number of update sweeps to run over the lattice
* `cutoff`    : The percentage of measurements to throw away, i.e. throw away values during
    thermalization
* `start_type`: `:cold` starts the lattice with all spins aligned, `:hot` all spins random
* `seed`      : Seed value for the random number generator

Returns
* A new instance of `Potts`

Attributes
* `σ`          : Lattice containing the spins, values are ∈ {1, ..., q}
* `params`     : Parameters type, contains relevant parameters for the simulation
* `observables`: Observables type, contains relevant observables data for the simulation
* `rng`        : Random number generator

# Examples
```julia-repl
julia> Potts(8, 32, log(1 + √8))
Potts
L          = 32
β          = 0.3
n_sweeps   = 100000
start_type = cold
seed       = 8
q          = 8
```
"""
struct Potts <: AbstractMagnetismModel
    σ::Array{Int8,2}
    params::Parameters
    observables::Observables
    rng::AbstractRNG

    function Potts(
        q::Int,
        L::Int,
        β::Real;
        n_sweeps::Int = 10^5,
        cutoff::Float64 = 0.2,
        start_type::Symbol = :cold,
        seed::Int = 8,
    )

        params = Parameters(L, β, n_sweeps, cutoff, start_type, seed, q)
        rng = MersenneTwister(seed)
        σ = initialize_σ_Potts(q, L, start_type, rng)
        observables = Observables(Int)

        return new(σ, params, observables, rng)
    end
end

"""
    XY(
        L         ::Int,
        β         ::Real;
        n_sweeps  ::Int    =10^5,
        cutoff    ::Float64=0.2,
        start_type::Symbol =:cold,
        seed      ::Int    =8,
        )

Type representing an XY model.

Arguments
* `L`: The side length of the lattice, i.e. for a square lattice there will be `L^2` spins
* `β`: `1 / kT`, i.e. the inverse temperature

Keyword Arguments
* `n_sweeps`  : The number of update sweeps to run over the lattice
* `cutoff`    : The percentage of measurements to throw away, i.e. throw away values during
    thermalization
* `start_type`: `:cold` starts the lattice with all spins aligned, `:hot` all spins random
* `seed`      : Seed value for the random number generator

Returns
* A new instance of `XY`

Attributes
* `σ`          : Lattice containing the spins, values are 2D unit vectors
* `params`     : Parameters type, contains relevant parameters for the simulation
* `observables`: Observables type, contains relevant observables data for the simulation
* `rng`        : Random number generator

# Examples
```julia-repl
julia> XY(32, 1)
XY
L          = 32
β          = 1.0
n_sweeps   = 100000
start_type = cold
seed       = 8
q          = 0
```
"""
struct XY <: AbstractMagnetismModel
    σ::Array{NTuple{2,Float64},2}
    params::Parameters
    observables::Observables
    rng::AbstractRNG

    function XY(
        L::Int,
        β::Real;
        n_sweeps::Int = 10^5,
        cutoff::Float64 = 0.2,
        start_type::Symbol = :cold,
        seed::Int = 8,
    )

        params = Parameters(L, β, n_sweeps, cutoff, start_type, seed, 0)
        rng = MersenneTwister(seed)
        σ = initialize_σ_XY(L, start_type, rng)
        observables = Observables(Float64)

        return new(σ, params, observables, rng)
    end
end


text:
#function ancestors(node::Array{Any,1},nd::Int)
#    anc=Array(Int,0)
#    parents=node[nd].parents
#    append!(anc,parents)
#    done=false
#    while !done
#        newparents=Array(Int,0)
#        for i in parents
#            append!(newparents,node[i].parents)
#        end
#        if isempty(newparents)
#            done=true
#        else
#            append!(anc,newparents)
#            parents=newparents
#        end
#    end
#    return sort(union(anc);rev=true)
#
#end
#export ancestors

function ancestors(node::Array{Any,1},nd::Int)
    anc=Array(Int,0)
    parents=node[nd].parents
    append!(anc,parents)
    done=false
    while !done
        newparents=Array(Int,0)
        for i in parents
            append!(newparents,node[i].parents)
        end
        newparents=setdiff(newparents,anc)
        if isempty(newparents)
            done=true
        else
            append!(anc,newparents)
            parents=newparents
        end
    end
    return sort(union(anc);rev=true)

end
export ancestors

function Parameters(net)
#    intersect(ancestors(net.node,net.FunctionNode),find(map(x->x.returnderivative,net.node))) # node indices that are parameters
#    intersect(ancestors(net.node,net.FunctionNode),find(map(x->x.returnderivative,net.node[net.validnodes]))) # node indices that are parameters

    good=[]
    for i in net.validnodes
        if net.node[i].returnderivative
            push!(good,i)
        end
    end
    intersect(ancestors(net.node,net.FunctionNode),good) # node indices that are parameters
end

    export Parameters

    function HasParents(node::ADnode)
        if length(node.parents)>0
            return true
        else
            return false
        end
    end
    export HasParents

function ForwardPassList!(net;ExcludeNodes=[])
    forwardlist=Array(Int,0)
    node=net.node
    N=length(node)
    exclude=Array(Int,length(ExcludeNodes))
    for i=1:length(ExcludeNodes)
        if isa(ExcludeNodes[i],ADnode)
            exclude[i]=ExcludeNodes[i].index
        else
            exclude[i]=ExcludeNodes[i]
        end
    end

    for i in intersect(setdiff(1:N,exclude),net.validnodes)
        thisnode=node[i]
        if length(findin(thisnode.parents,exclude))>0 # can be precomputed
            push!(exclude,i)# check
        else
#            if !(thisnode.input)
            if HasParents(thisnode)
                push!(forwardlist,i)
            end
        end
    end
    net.ForwardPassList=forwardlist
#    return net
end
export ForwardPassList!

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# include("hyperelasticity.jl")

ts = ExplicitSolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globaldata, ts)
Fext = compute_external_force(domain, globaldata, ts)

d0 = zeros(2domain.nnodes)
v0 = zeros(2domain.nnodes)
a0 = zeros(2domain.nnodes)

# general neural network
function nn_law(ε)
  ae(ε, [20,20,20,3])
end

# free energy
function nn_law(ε)
  φ = squeeze(ae(ε, [20,20,20,1]))
  tf.gradients(φ, ε)
end

# parameter
c1 = Variable(1.0) * 1e-1
c2 = Variable(1.0) * 1e-3
function nn_law(ε)
  compute_stress_rivlin_saunders(ε, c1, c2)
end

d, v, a= ExplicitSolver(globaldata, domain, d0, v0, a0, Δt, NT, nn_law, Fext, ubd, abd; strain_type="finite")

idx = @. (0:n)*m + 1
idx = [idx; (@. idx + domain.nnodes)]
loss = sum((d[:, idx] - d_[:, idx])^2)
sess = Session(); init(sess)
@info run(sess, loss)
# error()

BFGS!(sess, loss)

# for visualizing the parameters
# sol = []
# cb = (vs, iter, loss)->begin
#   push!(sol, vs)
#   printstyled("[#iter $iter] a = $vs, loss=$loss\n", color=:green)
# end
# BFGS!(sess, loss, callback = cb, vars = [c1, c2])

# plot([x[1] for x in sol], "+--", color="red", label="\$C1\$")
# plot(ones(length(sol))*1e-1, "-", color="k", alpha=0.3, label="Reference Value (C1 and C2)")
# plot([x[2] for x in sol], "x--", color ="green", label="\$C2\$")
# xlabel("Iterations")
# legend()


text:
"""
sample!([rng,] a, k; exclude = ())

Sample `k` element from array `a` without repetition and eventually excluding elements in `exclude`.
Pay attention, it changes the order of the elements in `a`.
"""
function sample!(rng::AbstractRNG, a::AbstractArray, k::Integer; exclude = ())
    length(a) < k + length(exclude) && error("Array too short.")
    res = Vector{eltype(a)}()
    sizehint!(res, k)
    n = length(a)
    i = 1
    while length(res) < k
        r = rand(rng, 1:n-i+1)
        if !(a[r] in exclude)
            push!(res, a[r])
            a[r], a[n-i+1] = a[n-i+1], a[r]
            i += 1
        end
    end
    return res
end

sample!(a::AbstractArray, k::Integer; exclude = ()) = sample!(getRNG(), a, k; exclude = exclude)

getRNG(seed::Integer = -1) = seed >= 0 ? MersenneTwister(seed) : Random.GLOBAL_RNG

# modified from http://stackoverflow.com/questions/25678112/insert-item-into-a-sorted-list-with-julia-with-and-without-duplicates
# returns true if insert succeeded, false if it was a duplicate
_insert_and_dedup!(v::Vector, x) = isempty(splice!(v, searchsorted(v,x), x))

function countfirst(itr, v)
    c = 0
    found = false
    @inbounds for x in itr
        c +=1
        x == v && (found=true; break)
    end
    return found ? c : 0
end

# for generic iterables with length
myrand(itr) = nth(itr, rand(1:length(itr)))

#used in flow and dismantling
struct GreaterThan2 <: Base.Order.Ordering end
Base.Order.lt(c::GreaterThan2, x, y) = x[2] > y[2]

struct LessThan2 <: Base.Order.Ordering end
Base.Order.lt(c::LessThan2, x, y) = x[2] < y[2]

signedtype(::Type{T}) where {T<:Integer} = typeof(signed(T(0)))
signedtype(::Type{T}) where {T<:AbstractFloat} = T


text:
function discretize(h5::String; inplace = 0, name = "特征", mpiarg = ``)
    julia = joinpath(Sys.BINDIR, Base.julia_exename())
    h5disc = joinpath(@__DIR__, "h5disc.jl")
    project = joinpath(@__DIR__, "..")
    run(`mpiexec $mpiarg $julia --project=$project $h5disc --name $name --inplace $inplace $h5`)
    h5′ = replace(h5, ".h5" => "_UInt8.h5")
    inplace == 1 ? h5 : h5′
end

discretize!(a...; ka...) = discretize(a...; inplace = 1, ka...)

function discretize(x::AbstractArray; ka...)
    h5 = randstring() * ".h5"
    h5write(h5, "x", x)
    discretize!(h5; name = "x", ka...)
    x8 = h5read(h5, "x")
    bin_edges = h5read(h5, "bin_edges")
    rm(h5, force = true)
    return x8, bin_edges
end

discretize(x, edge::Vector) = map(z -> searchsortedfirst(edge, z), x)

function discretize(x, bin_edges::Matrix)
    x′ = reshape(x, size(x, 1), :)
    x8′ = zeros(UInt8, size(x′))
    for f in 1:size(x8′, 1)
        x8′[f, :] .= discretize(x′[f, :], bin_edges[:, f])
    end
    reshape(x8′, size(x))
end

function undiscretize(x8, edge::Vector)
    @inbounds map(x8) do z
        if z == 0x00
            edge[1]
        elseif z == 0xff
            edge[end]
        else
            l = edge[z]
            r = edge[z + 1]
            (l + r) / 2
        end
    end
end

function undiscretize(x8, bin_edges::Matrix)
    x8′ = reshape(x8, size(x8, 1), :)
    x′ = zeros(Float32, size(x8′))
    for f in 1:size(x′, 1)
        x′[f, :] .= undiscretize(x8′[f, :], bin_edges[:, f])
    end
    reshape(x′, size(x8))
end


text:
abstract AbstractNewtonSolver{T} <: NonlinearSolver{T}

@define newton_solver_variables begin
    x::Vector{T}

    Fparams::TF
    Jparams::TJ

    linear::TL

    params::NonlinearSolverParameters{T}
    status::NonlinearSolverStatus{T}
end

function setInitialConditions!{T}(s::AbstractNewtonSolver{T}, x₀::Vector{T})
    s.x[:] = x₀
end


text:
function add_milp_relaxation(
    f::Function,
    partition::Vector{Float64},
    error_tolerance::Float64,
    m::JuMP.AbstractModel;
    name = "1",
)
    relaxation, _ = construct_milp_relaxation(
        f,
        partition,
        length_tolerance = 1e-5,
        derivative_tolerance = 1e-5,
        error_tolerance = error_tolerance,
    )
    num_variables = get_num_variables(relaxation)
    lb, ub = get_variable_bounds(relaxation)
    names = get_variable_names(relaxation)
    binary = relaxation.binary
    x = @variable(
        m,
        [i = 1:num_variables],
        lower_bound = lb[i],
        upper_bound = ub[i],
        binary = Bool(binary[i]),
        base_name = names[i] * "_" * name
    )
    A, b = get_eq_constraint_matrices(relaxation)
    @constraint(m, A * x .== b)
    A, b = get_leq_constraint_matrices(relaxation)
    @constraint(m, A * x .<= b)
    return x, relaxation.x_index, relaxation.y_index
end

function add_lp_relaxation(
    f::Function,
    partition::Vector{Float64},
    error_tolerance::Float64,
    m::JuMP.AbstractModel;
    name = "1",
)
    relaxation, _ = construct_lp_relaxation(
        f,
        partition,
        length_tolerance = 1e-5,
        derivative_tolerance = 1e-5,
        error_tolerance = error_tolerance,
    )
    num_variables = get_num_variables(relaxation)
    lb, ub = get_variable_bounds(relaxation)
    names = get_variable_names(relaxation)
    x = @variable(
        m,
        [i = 1:num_variables],
        lower_bound = lb[i],
        upper_bound = ub[i],
        base_name = names[i] * "_" * name
    )
    A, b = get_eq_constraint_matrices(relaxation)
    @constraint(m, A * x .== b)
    return x, relaxation.x_index, relaxation.y_index
end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
module supply_allocation

using JuMP
using Cbc
using CSV

function solve_supply_allocation()

    #------
    # DATA
    #------

    n = 4;
    m = 3;
    Q = [30 15 20 35];
    d = [40 10 35];

    c = [
        40	39	40;
        35	36	31;
        33	43	42;
        37	35	38
    ];

    #------
    # MODEL
    #------

    model = Model(with_optimizer(Cbc.Optimizer));

    @variable(model, x[1:n,1:m] >= 0)

    @objective(model, Min, sum(c[i,j] * x[i,j] for i in 1:n, j in 1:m));

    @constraint(model,[j = 1:m],sum(x[i,j] for i in 1:n) == d[j]);

    @constraint(model,[i = 1:n],sum(x[i,j] for j in 1:m) <= Q[i]);

    f = open("supply_allocation.lp", "w")
    print(f, model)
    close(f)

    #-------
    # SOLVE
    #-------

    optimize!(model)

    #------------------------
    # WRITE SOLUTION TO FILE
    #------------------------

    f = open("supply_allocation.csv","w");

    print(f,",");
    print(f,"Munich,Salzburg,Prague");
    println(f,"");

    ports = ["Rotterdam" "Wilhelmshaven" "Bremerhaven" "Hamburg" ];

    for i = 1:n
        print(f,ports[i],",");
        for j = 1:m
            print(f,value(x[i,j]),",");
        end
        println(f,"");
    end

    close(f);

end

solve_supply_allocation();

end


text:
export freqprior

struct Beta <: Stats{Bool}
  α::Float64
  β::Float64
end
update(d::Beta,x::Bool,w=1.0) = x ? Beta(d.α+w,d.β) : Beta(d.α,d.β+w)
downdate(d::Beta,x::Bool,w=1.0) = x ? Beta(d.α-w,d.β) : Beta(d.α,d.β-w)
mult(d::Beta,x) = d*x
Base.:(*)(d::Beta,x) = Beta(d.α*x,d.β*x)
Base.:(+)(a::Beta,b::Beta) = Beta(a.α+b.α,a.β+b.β)
pdf(d::Beta,x) = x ? d.α / (d.α + d.β) : d.β / (d.α + d.β)
logpdf(d::Beta,x) = log(pdf(d,x))

struct BinomialCond{T}
  at::Dict{T,Beta}
end
BinomialCond(pairs::Pair{T}...) where T = BinomialCond{T}(Dict(pairs...))
update!(d::BinomialCond,cond,x::Bool,w=1.0) =
  d.at[cond] = update(d.at[cond],x,w)
downdate!(d::BinomialCond,cond,x::Bool,w=1.0) =
  d.at[cond] = downdate(d.at[cond],x,w)
function mult!(d::BinomialCond,x)
  for cond in keys(d.at)
    d.at[cond] = mult(d.at[cond],x)
  end
end
Base.:(+)(a::BinomialCond,b::BinomialCond) =
  BinomialCond(merge(+,a.at,b.at))
function logpdf(d::BinomialCond,cond,x)
  p = d.at[cond].α / (d.at[cond].α + d.at[cond].β)
  x ? log(p) : log(1-p)
end

function freqprior(bias,N)
  Beta(N*(1-(bias/2+0.5)),N)
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


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

# %%
function pimc(n)
    c = 0
    for _ in 1:n
        c += rand()^2 + rand()^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc(10^9)
@time pimc(10^9)
@time pimc(10^9)

# %%
using Random: default_rng

function pimc2(n)
    rng = default_rng()
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc2(10^9)
@time pimc2(10^9)
@time pimc2(10^9)

# %%
using Random: MersenneTwister

function pimc3(n)
    rng = MersenneTwister()
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc3(10^9)
@time pimc3(10^9)
@time pimc3(10^9)

# %%
function pimc(n)
    c = 0
    for _ in 1:n
        c += rand()^2 + rand()^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc(10^9)
@time pimc(10^9)
@time pimc(10^9)

# %%
using Random: default_rng

function pimc2(n)
    rng = default_rng()
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc2(10^9)
@time pimc2(10^9)
@time pimc2(10^9)

# %%
using Random: MersenneTwister

function pimc3(n)
    rng = MersenneTwister()
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc3(10^9)
@time pimc3(10^9)
@time pimc3(10^9)

# %%
function pimc(n)
    c = 0
    for _ in 1:n
        c += rand()^2 + rand()^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc(10^9)
@time pimc(10^9)
@time pimc(10^9)

# %%
using Random: default_rng

function pimc2(n)
    rng = default_rng()
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc2(10^9)
@time pimc2(10^9)
@time pimc2(10^9)

# %%
using Random: MersenneTwister

function pimc3(n)
    rng = MersenneTwister()
    c = 0
    for _ in 1:n
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/n
end

@show VERSION
@time pimc3(10^9)
@time pimc3(10^9)
@time pimc3(10^9)

# %%

max_stars_repo_path: 0004/rand bench.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0004/rand bench.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0004/rand bench.jl

max_forks_repo_name: genkuroki/public


text:
# ---
# title: 427. Construct Quad Tree
# id: problem427
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/construct-quad-tree/description/>
# hidden: true
# ---
#
# Given a `n * n` matrix `grid` of `0's` and `1's` only. We want to represent
# the `grid` with a Quad-Tree.
#
# Return _the root of the Quad-Tree_ representing the `grid`.
#
# Notice that you can assign the value of a node to **True** or **False** when
# `isLeaf` is **False** , and both are **accepted** in the answer.
#
# A Quad-Tree is a tree data structure in which each internal node has exactly
# four children. Besides, each node has two attributes:
#
#   * `val`: True if the node represents a grid of 1's or False if the node represents a grid of 0's.
#   * `isLeaf`: True if the node is leaf node on the tree or False if the node has the four children.
#
#
#
#     class Node {
#         public boolean val;
#         public boolean isLeaf;
#         public Node topLeft;
#         public Node topRight;
#         public Node bottomLeft;
#         public Node bottomRight;
#     }
#
# We can construct a Quad-Tree from a two-dimensional area using the following
# steps:
#
#   1. If the current grid has the same value (i.e all `1's` or all `0's`) set `isLeaf` True and set `val` to the value of the grid and set the four children to Null and stop.
#   2. If the current grid has different values, set `isLeaf` to False and set `val` to any value and divide the current grid into four sub-grids as shown in the photo.
#   3. Recurse for each of the children with the proper sub-grid.
#
# ![](https://assets.leetcode.com/uploads/2020/02/11/new_top.png)
#
# If you want to know more about the Quad-Tree, you can refer to the
# [wiki](https://en.wikipedia.org/wiki/Quadtree).
#
# **Quad-Tree  format:**
#
# The output represents the serialized format of a Quad-Tree using level order
# traversal, where `null` signifies a path terminator where no node exists
# below.
#
# It is very similar to the serialization of the binary tree. The only
# difference is that the node is represented as a list `[isLeaf, val]`.
#
# If the value of `isLeaf` or `val` is True we represent it as **1** in the list
# `[isLeaf, val]` and if the value of `isLeaf` or `val` is False we represent it
# as **0**.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/11/grid1.png)
#
#
#
#     Input: grid = [[0,1],[1,0]]
#     Output: [[0,1],[1,0],[1,1],[1,1],[1,0]]
#     Explanation: The explanation of this example is shown below:
#     Notice that 0 represnts False and 1 represents True in the photo representing the Quad-Tree.
#     ![](https://assets.leetcode.com/uploads/2020/02/12/e1tree.png)
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/12/e2mat.png)
#
#
#
#     Input: grid = [[1,1,1,1,0,0,0,0],[1,1,1,1,0,0,0,0],[1,1,1,1,1,1,1,1],[1,1,1,1,1,1,1,1],[1,1,1,1,0,0,0,0],[1,1,1,1,0,0,0,0],[1,1,1,1,0,0,0,0],[1,1,1,1,0,0,0,0]]
#     Output: [[0,1],[1,1],[0,1],[1,1],[1,0],null,null,null,null,[1,0],[1,0],[1,1],[1,1]]
#     Explanation: All values in the grid are not the same. We divide the grid into four sub-grids.
#     The topLeft, bottomLeft and bottomRight each has the same value.
#     The topRight have different values so we divide it into 4 sub-grids where each has the same value.
#     Explanation is shown in the photo below:
#     ![](https://assets.leetcode.com/uploads/2020/02/12/e2tree.png)
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[1,1],[1,1]]
#     Output: [[1,1]]
#
#
# **Example 4:**
#
#
#
#     Input: grid = [[0]]
#     Output: [[1,0]]
#
#
# **Example 5:**
#
#
#
#     Input: grid = [[1,1,0,0],[1,1,0,0],[0,0,1,1],[0,0,1,1]]
#     Output: [[0,1],[1,1],[1,0],[1,0],[1,1]]
#
#
#
#
# **Constraints:**
#
#   * `n == grid.length == grid[i].length`
#   * `n == 2^x` where `0 <= x <= 6`
#
#
## @lc code=start
using LeetCode

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


text:
macro llvmir2svg(expr,filepath="")
    if expr.head == :call
        assigns = []
        vars = []
        for i in expr.args[2:end]
            var = gensym()
            push!(vars,var)
            push!(assigns,:($var = $(esc(i))))
        end
        return quote
            $(assigns...)
            llvm2graphfile($(esc(expr.args[1])),Base.typesof($(vars...)),$(esc(filepath)),"svg")
        end
    end
end

macro llvmir2html(expr,filepath="")
    if expr.head == :call
        assigns = []
        vars = []
        for i in expr.args[2:end]
            var = gensym()
            push!(vars,var)
            push!(assigns,:($var = $(esc(i))))
        end
        return quote
            $(assigns...)
            llvm2graphfile($(esc(expr.args[1])),Base.typesof($(vars...)),$(esc(filepath)),"html")
        end
    end
end

max_stars_repo_name: ChenNingCong/IRViewer

max_issues_repo_name: ChenNingCong/IRViewer

max_forks_repo_name: ChenNingCong/IRViewer


text:
using Pkg
using SHA
using TimeZones
using Base: AbstractVecOrTuple

timezones_built = false

##############
# CONVERSION #
##############

# NOTE: `jl` referes to a subset of Julia that can be serialized to or
# deserialized from JSON with ease

jl_to_json(j) = JSON.json(j)

json_to_jl(j) = JSON.parse(j)

key_to_jl(key) = reinterpret(UInt8, hash(string(key))) |> String
axis_to_jl(axis) = reinterpret(UInt8, hash(string(key))) |> String

total_memory_usage(val) =
    begin
        size = Base.summarysize(val)
        # TODO: Maybe make this larger
        if size ≤ 128
            0
        else
            size
        end
    end

# NOTE: This function is shared between the client library and the PT library
function indexapply(op, objs...; index::Integer=1)
    lists = [obj for obj in objs if (obj isa AbstractVector || obj isa Tuple)]
    length(lists) > 0 || throw(ArgumentError("Expected at least one tuple as input"))
    index = index isa Colon ? length(first(lists)) : index
    operands = [((obj isa AbstractVector || obj isa Tuple) ? obj[index] : obj) for obj in objs]
    indexres = op(operands...)
    res = first(lists)
    if first(lists) isa Tuple
        res = [res...]
        res[index] = indexres
        Tuple(res)
    else
        res = copy(res)
        res[index] = indexres
        res
    end
end

# converts give time as String to local timezone and returns DateTime
function parse_time(time)
    global timezones_built
    if !timezones_built
        # TODO: Add a try-catch here and then if building fails, @warn that we
        # can't convert to local timezone and simply return
        TimeZones.build()
        timezones_built = true
    end
    DateTime(astimezone(ZonedDateTime(time * "0000", "yyyy-mm-dd-HH:MM:SSzzzz"), localzone()))
end

function s3_bucket_arn_to_name(s3_bucket_arn)
    # Get s3 bucket name from arn
    s3_bucket_name = last(split(s3_bucket_arn, ":"))
    if endswith(s3_bucket_name, "/")
        s3_bucket_name = s3_bucket_name[1:end-1]
    elseif endswith(s3_bucket_name, "/*")
        s3_bucket_name = s3_bucket_name[1:end-2]
    elseif endswith(s3_bucket_name, "*")
        s3_bucket_name = s3_bucket_name[1:end-1]
    end
    return s3_bucket_name
end

function s3_bucket_name_to_arn(s3_bucket_name)
    # Get s3 bucket arn from name
    s3_bucket_arn = s3_bucket_name
    if endswith(s3_bucket_arn, "/")
        s3_bucket_arn = s3_bucket_arn[1:end-1]
    elseif endswith(s3_bucket_arn, "/*")
        s3_bucket_arn = s3_bucket_arn[1:end-2]
    elseif endswith(s3_bucket_arn, "*")
        s3_bucket_arn = s3_bucket_arn[1:end-1]
    end
    return s3_bucket_arn
end

##################
# AUTHENTICATION #
##################

# Process-local configuration for the account being used. It wouldn't be hard
# to but there shouldn't be any reason to make this thread-local (since only
# one account should be being used per workstation or per server where
# Banyan.jl may be being used). However, wrapping this in a mutex to ensure
# synchronized mutation in this module would be a good TODO.
global banyan_config = nothing
global aws_config_in_usage = nothing

function load_config(banyanconfig_path=nothing)
    global banyan_config

    if banyanconfig_path == nothing
        banyanconfig_path = joinpath(homedir(), ".banyan", "banyanconfig.toml")
    end
    if isfile(banyanconfig_path)
        banyan_config = TOML.parsefile(banyanconfig_path)
    end
end

function write_config(banyanconfig_path=nothing)
    global banyan_config

    # Write to banyanconfig.toml
    if banyanconfig_path == nothing
        banyanconfig_path = joinpath(homedir(), ".banyan", "banyanconfig.toml")
    end
    mkpath(joinpath(homedir(), ".banyan"))
    f = open(banyanconfig_path, "w")
    TOML.print(f, banyan_config)
    close(f)
end

if_in_or(key, obj, el = nothing) =
    if key in keys(obj) && !isnothing(obj[key])
        obj[key]
    else
        el
    end

function configure(; kwargs...)
    # This function allows for users to configure their authentication.
    # Authentication details are then saved in
    # `$HOME/.banyan/banyanconfig.toml` so they don't have to be entered in again
    # each time a program using the Banyan client library is run

    # Credentials are checked in the following locations in this order:
    #   1) function arguments, specified in kwargs
    #   2) environment variables
    #   3) `$HOME/.banyan/banyanconfig.toml`

    # Load arguments
    # If an argument is optional (e.g., `ec2_key_pair_name`), we default
    # to 0 to indicate that the argument was not specified. This is so that
    # we can differentiate between a user explicitly providing `nothing` as
    # the value for an arg, versus a default.
    kwargs = Dict(kwargs)
    user_id = if_in_or(:user_id, kwargs)
    api_key = if_in_or(:api_key, kwargs)
    ec2_key_pair_name = if_in_or(:ec2_key_pair_name, kwargs, 0)
    banyanconfig_path = if_in_or(:banyanconfig_path, kwargs)

    # Load config
    load_config(banyanconfig_path)
    global banyan_config

    # Check environment variables
    if isnothing(user_id) && haskey(ENV, "BANYAN_USER_ID")
        user_id = ENV["BANYAN_USER_ID"]
    end
    if isnothing(api_key) && haskey(ENV, "BANYAN_API_KEY")
        api_key = ENV["BANYAN_API_KEY"]
    end
    if ec2_key_pair_name == 0 && haskey(ENV, "BANYAN_EC2_KEY_PAIR_NAME")
        api_key = ENV["BANYAN_EC2_KEY_PAIR_NAME"]
    end

    # Check banyanconfig file
    if isnothing(user_id) && haskey(banyan_config, "banyan") && haskey(banyan_config["banyan"], "user_id")
        user_id = banyan_config["banyan"]["user_id"]
    end
    if isnothing(api_key) && haskey(banyan_config, "banyan") && haskey(banyan_config["banyan"], "api_key")
        api_key = banyan_config["banyan"]["api_key"]
    end

    # Initialize
    is_modified = false
    is_valid = true

    # Ensure a configuration has been created or can be created. Otherwise,
    # return nothing
    if isnothing(banyan_config) || isempty(banyan_config)
        if !isnothing(user_id) && !isnothing(api_key)
            banyan_config = Dict(
                "banyan" =>
                    Dict("user_id" => user_id, "api_key" => api_key),
                "aws" => Dict(),
            )
            is_modified = true
        else
            error("Your user ID and API key must be specified using either keyword arguments, environment variables, or banyanconfig.toml")
        end
    end

    # Check for changes in required
    if !isnothing(user_id) &&
        (user_id != banyan_config["banyan"]["user_id"])
         banyan_config["banyan"]["user_id"] = user_id
         is_modified = true
     end
    if !isnothing(api_key) && (api_key != banyan_config["banyan"]["api_key"])
        banyan_config["banyan"]["api_key"] = api_key
        is_modified = true
    end

    # Check for changes in other args

    # aws.ec2_key_pair_name
    if isnothing(ec2_key_pair_name)
        delete!(banyan_config["aws"], "ec2_key_pair_name")
    elseif (ec2_key_pair_name != 0) && (
        !(haskey(banyan_config["aws"], "ec2_key_pair_name")) ||
        ec2_key_pair_name != banyan_config["aws"]["ec2_key_pair_name"]
    )
        banyan_config["aws"]["ec2_key_pair_name"] = ec2_key_pair_name
        is_modified = true
    end

    # # aws.region
    # if !isnothing(region) && (
    #     !(haskey(banyan_config["aws"], "region")) ||
    #     region != banyan_config["aws"]["region"]
    # )
    #     banyan_config["aws"]["region"] = region
    #     is_modified = true
    # end

    # Update config file if it was modified
    if is_modified
        write_config(banyanconfig_path)
    end

    return banyan_config
end

function get_aws_config()
    global aws_config_in_usage

    # Get AWS configuration
    if isnothing(aws_config_in_usage)
        # Get region according to ENV, then credentials, then config files
        profile = get(ENV, "AWS_DEFAULT_PROFILE", get(ENV, "AWS_DEFAULT_PROFILE", "default"))
        region = get(ENV, "AWS_DEFAULT_REGION", "")
        if region == ""
            try
                configfile = read(Inifile(), joinpath(homedir(), ".aws", "config"))
                region = _get_ini_value(configfile, profile, "region", default_value="")
            catch
            end
        end
        if region == ""
            try
                credentialsfile = read(Inifile(), joinpath(homedir(), ".aws", "credentials"))
                region = _get_ini_value(credentialsfile, profile, "region", default_value="")
            catch
            end
        end

        if region == ""
            throw(ErrorException("Could not discover AWS region to use from looking at AWS_PROFILE, AWS_DEFAULT_PROFILE, AWS_DEFAULT_REGION, HOME/.aws/credentials, and HOME/.aws/config"))
        end

        aws_config_in_usage = Dict(
            :creds => AWSCredentials(),
            :region => region
        )
    end

    # # Use default location if needed
    # if !haskey(aws_config_in_usage, :region)
    #     @warn "Using default AWS region of us-west-2 in \$HOME/.banyan/banyanconfig.toml"
    #     aws_config_in_usage[:region] = "us-west-2"
    # end

    # Convert to dictionary and return

    aws_config_in_usage
end

get_aws_config_region() = get_aws_config()[:region]

#########################
# ENVIRONMENT VARIABLES #
#########################

is_debug_on() = get(ENV, "JULIA_DEBUG", "") == "Banyan"

macro in_env(key)
    return :(string("BANYAN_", getpid(), "_", $key) in keys(ENV))
end

macro env(key)
    return :(ENV[string("BANYAN_", getpid(), "_", $key)])
end

macro delete_in_env(key)
    return :(delete!(ENV, string("BANYAN_", getpid(), "_", $key)))
end

################
# API REQUESTS #
################

method_to_string(method) = begin
    if method == :create_cluster
        "create-cluster"
    elseif method == :destroy_cluster
        "destroy-cluster"
    elseif method == :describe_clusters
        "describe-clusters"
    elseif method == :create_job
        "create-job"
    elseif method == :destroy_job
        "destroy-job"
    elseif method == :describe_jobs
        "describe-jobs"
    elseif method == :evaluate
        "evaluate"
    elseif method == :update_cluster
        "update-cluster"
    elseif method == :set_cluster_ready
        "set-cluster-ready"
    end
end

"""
Sends given request with given content
"""
function request_body(url::AbstractString; kwargs...)
    global downloader
    resp = nothing
    body = sprint() do output
        resp = request(url; output=output, throw=false, downloader=downloader, kwargs...)
    end
    return resp, body
end

function request_json(url::AbstractString; kwargs...)
    resp, body = request_body(url; kwargs...)
    return resp, JSON.parse(body)
end

# Sends an HTTP request to the Banyan API and returns the
# parsed response. Sends the provided content as the body of
# the message and additionally adds the User ID and API Key,
# which are required on all requests for authentication. Additionally,
# a debug flag is sent. An exception is thrown is the HTTP
# requests returns  a 403, 500, or 504 HTTP error cdode. If the
# request times out, a warning message is printed out and `nothing`
# is returned.
# It is up to the caller to handle the case where the HTTP request
# times out and `nothing` is returned.
function send_request_get_response(method, content::Dict)
    # Prepare request
    configuration = load_config()
    user_id = configuration["banyan"]["user_id"]
    api_key = configuration["banyan"]["api_key"]
    content["debug"] = is_debug_on()
    url = string(BANYAN_API_ENDPOINT, method_to_string(method))
    headers = [
        "content-type" => "application/json",
        "Username-APIKey" => "$user_id-$api_key",
    ]
    resp, data = request_json(
	    url, input=IOBuffer(JSON.json(content)), method="POST", headers=headers
    )
    if resp.status == 403
        error("Please use a valid user ID and API key. Sign into the dashboard to retrieve these credentials.")
    elseif resp.status == 504
        # HTTP request timed out, for example
        if isa(data, Dict) && haskey(data, "message")
            data = data["message"]
        end
        @error data
        return nothing
    elseif resp.status == 500 || resp.status == 504
        error(data)
    elseif resp.status == 502
        error("Sorry there has been an error. Please contact support.")
    end
    return data

end

function send_request_get_response_using_http(method, content::Dict)
    # Prepare request
    # content = convert(Dict{Any, Any}, content)
    configuration = load_config()
    user_id = configuration["banyan"]["user_id"]
    api_key = configuration["banyan"]["api_key"]
    content["debug"] = is_debug_on()
    url = string(BANYAN_API_ENDPOINT, method_to_string(method))
    headers = (
        ("content-type", "application/json"),
        ("Username-APIKey", "$user_id-$api_key"),
    )

    # Post and return response
    try
        response = HTTP.post(url, headers, JSON.json(content))
        body = String(response.body)
        return JSON.parse(body)
    catch e
        if e isa HTTP.ExceptionRequest.StatusError
            if e.response.status == 403
                throw(
                    ErrorException(
                        "Please set a valid api_key. Sign in to the dashboard to retrieve your api key.",
                    ),
                )
            end
            if e.response.status != 504
                throw(ErrorException(String(take!(IOBuffer(e.response.body)))))
            # elseif method == :create_cluster
            #     # println(
            #     #     "Cluster creation in progress. Please check dashboard to view status.",
            #     # )
            # elseif method == :create_job
            #     # println(
            #     #     "Job creation in progress. Please check dashboard to view status.",
            #     # )
            # elseif method == :evaluate
            #     # println(
            #     #     "Evaluation is in progress. Please check dashboard to view status.",
            #     # )
            end
            rethrow()
        else
            rethrow()
        end
    end
end

#########
# FILES #
#########

function load_json(path::String)
    if startswith(path, "file://")
        if !isfile(path[8:end])
            error("File $path does not exist")
        end
        JSON.parsefile(path[8:end])
    elseif startswith(path, "s3://")
        error("S3 path not currently supported")
        # JSON.parsefile(S3Path(path, config=get_aws_config()))
    elseif startswith(path, "http://") || startswith(path, "https://")
	    JSON.parse(String(HTTP.get(path).body))
    else
        error("Path $path must start with \"file://\", \"s3://\", or \"http(s)://\"")
    end
end

# Loads file into String and returns
function load_file(path::String)
    if startswith(path, "file://")
        if !isfile(path[8:end])
            error("File $path does not exist")
        end
        String(read(open(path[8:end])))
    elseif startswith(path, "s3://")
        error("S3 path not currently supported")
        String(read(S3Path(path)))
    elseif startswith(path, "http://") || startswith(path, "https://")
        String(HTTP.get(path).body)
    else
        error("Path $path must start with \"file://\", \"s3://\", or \"http(s)://\"")
    end
end

####################################
# JULIA VERSION/PACKAGE MANAGEMENT #
####################################

function get_julia_version()
    return string(VERSION)
end

function get_loaded_packages()
    modules = map(
        m -> string(m),
        filter(
            m -> Main.eval(m) isa Module && !(m in [:Main, :Base, :Core, :InteractiveUtils, :IJulia]),
            names(Main, imported=true)
        )
    )
end

# Returns the directory in which the Project.toml file is located
function get_julia_environment_dir()
    return replace(Pkg.project().path, "Project.toml" => "")
end

# Returns SHA 256 of a string
function get_hash(s)
    return bytes2hex(sha256(s))
end

function format_bytes(bytes, decimals = 2)
    bytes == 0 && return "0 Bytes"
    k = 1024
    dm = decimals < 0 ? 0 : decimals
    sizes = ["Bytes", "KB", "MB", "GB", "TB", "PB", "EB", "ZB", "YB"]
    i = Base.convert(Int, floor(log(bytes) / log(k)))
    return string(round((bytes / ^(k, i)), digits = dm)) * " " * sizes[i+1]
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods."
mutable struct AbstractParametersSPARK{IT, DT, TT, D, S, R, P, equType <: NamedTuple, tabType <: AbstractTableau} <: Parameters{DT,TT}
    equs::equType
    tab::tabType
    Δt::TT

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function AbstractParametersSPARK{IT,DT,D,S,R,P}(equs::equType, tab::tabType, Δt::TT) where {IT,DT,TT,D,S,R,P,equType,tabType}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new{IT, DT, TT, D, S, R, P, equType, tabType}(equs, tab, Δt, zero(TT), q, p, λ)
    end

    function AbstractParametersSPARK{IT,DT,D}(equs::NamedTuple, tab::AbstractTableauSPARK, Δt::Number) where {IT,DT,D}
        AbstractParametersSPARK{IT, DT, D, tab.s, tab.r, tab.ρ}(equs, tab, Δt)
    end

    function AbstractParametersSPARK{IT,DT,D}(equs::NamedTuple, tab::Union{AbstractTableau,Tableau}, Δt::Number) where {IT,DT,D}
        AbstractParametersSPARK{IT, DT, D, tab.s, tab.r, 0}(equs, tab, Δt)
    end
end

function update_params!(params::AbstractParametersSPARK, sol::AtomicSolutionPDAE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
    params.λ .= sol.λ
end

@inline equation(int::AbstractIntegratorSPARK, i::Symbol) = parameters(int).equs[i]
@inline equations(int::AbstractIntegratorSPARK) = parameters(int).equs
@inline timestep(int::AbstractIntegratorSPARK) = parameters(int).Δt
@inline tableau(int::AbstractIntegratorSPARK) = parameters(int).tab

function Integrators.IntegratorCache{ST}(params::AbstractParametersSPARK{IT,DT,TT,D,S,R}; kwargs...) where {IT,ST,DT,TT,D,S,R}
    IntegratorCacheSPARK{ST,D,S,R}(; kwargs...)
end

@inline Integrators.CacheType(ST, params::AbstractParametersSPARK{IT,DT,TT,D,S,R}) where {IT,DT,TT,D,S,R} = IntegratorCacheSPARK{ST,D,S,R}


text:
#=
   Much of the activity bound to rational numbers and their arithmetic is tethered to integer multiplication.
   Rationals are a projectively exact representation within their computationally bounded domain.
   The ratio held as numerator and denominator is an exact representation of the rational value it projects.
   When applied to solve dynamical problems or to resolve the computational tractable in another manner,
   the coarse gapping that is given across steps of discrete rational availability is where approximation
   enters. This duality is the origin of their effectiveness.  We may play the representational equivalence
   off of the metrical coarseness and so apply this tool where it best shine.
=#

"""
    propinquity(x::`T`, y::`T`)::Union{`T`, widen(`T`)}

**x * y, without possiblity of overflow**

> `T` in {[U]Int8, [U]Int16, [U]Int32}

Implict provision of the return type selects either the args' shared type or the nearest larger type.
Using `Int16` args allows their product to be computed either as an `Int16` or as an `Int32`.
A wider type is used when necessary to preclude overflow, which itself protects the computation from
that sort of error.  When the originating type properly holds the product, the same type is used.
"""
propinquity(x::T, y::T) where {T<:FastInteger} =
    usewidemul(x,y) ? widemul(x,y) : x*y

"""
    Rapid determination of the working Type

```
We are given two range-limited integers of the same
performant system type, either signed or unsigned.
We are asked : "Do you assert that this type would
hold their product precisly, with full accuracy?"

`true` means
   It is certain that their product is within the standard discrete domain of the given type.

`false` means
   No assertion is made regarding the whether the given type would hold their product properly.
   It is asserted that at least of (a) or (b) hold.
     (a) It is known that use of with the given type nay promulage inaccuracies.
     (b) It is not known whether the precision of the given type covers this product.
```
""" usewidemul

@inline function usewidemul(x::T, y::T) where {T<:FastSigned}
    iszero( (leading_zeros(abs(x)) + leading_zeros(abs(y))) >> sizeof(T) )
end

@inline function usewidemul(x::T, y::T) where {T<:FastUnsigned}
    iszero( (leading_zeros(x) + leading_zeros(y)) >> sizeof(T) )
end