text:
"""
    _nearest_point(kdtree,sm_interpolant,x)

Helper function to calculate the distance between a sample point `x` and
all samples point in the kdtree.
"""
function _nearest_point(kdtree,sm_interpolant,x)
    x = permutedims(x') #Arrange in column vector
    v_tmp = [x; median(sm_interpolant(x))]
    v = SVector{length(v_tmp)}(v_tmp)
    _, dists = knn(kdtree, v, 1, true)

    return dists[1]
end

"""
    distance_infill(plan,samples,sm_interpolant)

Returns an anonymous function that calculates the distance between
sample point `x` and the closest point in `plan`.
The function value is included in the distance calculated by using the samples
and the estimated function value.
The distance is multiplied by -1 in order to maximise the distance between sample
point `x` and all other points.
"""
function distance_infill(plan,samples,sm_interpolant)

    #KD-tree for the plan as well as the samples from the plan
    combined_plan_samples = [plan;samples]
    kdtree = KDTree(combined_plan_samples)

    #The distance is calculated by picking a point in the design space, evaluating the
    #surrogate model at that point and checking the overall distance while
    #taking into consideration the function output. I.e constrained to the function surface

    return x -> -1*_nearest_point(kdtree,sm_interpolant,x)
end

"""
    minimum_infill(sm_interpolant)

Returns an anonymous function that calculates the minimum value of the surrogate
model ensemble predictions at position `x`.
"""
function minimum_infill(sm_interpolant)
    function (x)
        out = minimum(sm_interpolant(x))
        return out
    end
end

"""
    median_infill(sm_interpolant)

Returns an anonymous function that calculates the median value of the surrogate
model ensemble predictions at position `x`.
"""
function median_infill(sm_interpolant)
    function (x)
        out = median(sm_interpolant(x))
        return out
    end
end

"""
    mean_infill(sm_interpolant)

Returns an anonymous function that calculates the mean value of the surrogate
model ensemble predictions at position `x`.
"""
function mean_infill(sm_interpolant)
    function (x)
        out = mean(sm_interpolant(x))
        return out
    end
end

"""
    std_infill(sm_interpolant)

Returns an anonymous function that calculates the standard deviation of the surrogate
model ensemble predictions at position `x`. The std is multiplied by -1
in order to find the sample point `x` with the largest standard deviation.
"""
function std_infill(sm_interpolant)
    function (x)
        out = -1*std(sm_interpolant(x))
        return out
    end
end

"""
    weighted_med_std_infill(sm_interpolant)

Returns an anonymous function that calculates the weighted sum of the standard
deviation of the surrogate and the median location.
"""
function weighted_med_std_infill(sm_interpolant,w)

    df = std_infill(sm_interpolant)
    med = median_infill(sm_interpolant)

    function (x)
        out = w*med(x) + (1-w)*df(x)
        return out
    end
end

"""
    infill_objective(sm_interpolant,plan,samples,infill_funcs::Array{Symbol,1})

Returns the infill objective function based on the names supplied
in `infill_funcs`. Choose from `[:min,:median,:mean,:dist,:std]`.
"""
function infill_objective(sm_interpolant,plan,samples,infill_funcs::Array{Symbol,1})
    #Infill function options
    min_infill_fun = minimum_infill(sm_interpolant)
    median_infill_fun = median_infill(sm_interpolant)
    mean_infill_fun = mean_infill(sm_interpolant)
    dist_infill_fun = distance_infill(plan,samples,sm_interpolant)
    std_infill_fun = std_infill(sm_interpolant)
    w_std_med_infill_fun005 = weighted_med_std_infill(sm_interpolant,0.05)
    w_std_med_infill_fun01 = weighted_med_std_infill(sm_interpolant,0.1)
    w_std_med_infill_fun015 = weighted_med_std_infill(sm_interpolant,0.15)
    w_std_med_infill_fun02 = weighted_med_std_infill(sm_interpolant,0.2)
    w_std_med_infill_fun025 = weighted_med_std_infill(sm_interpolant,0.25)
    w_std_med_infill_fun03 = weighted_med_std_infill(sm_interpolant,0.3)
    w_std_med_infill_fun035 = weighted_med_std_infill(sm_interpolant,0.35)
    w_std_med_infill_fun04 = weighted_med_std_infill(sm_interpolant,0.4)
    w_std_med_infill_fun045 = weighted_med_std_infill(sm_interpolant,0.45)
    w_std_med_infill_fun05 = weighted_med_std_infill(sm_interpolant,0.5)
    w_std_med_infill_fun055 = weighted_med_std_infill(sm_interpolant,0.55)
    w_std_med_infill_fun06 = weighted_med_std_infill(sm_interpolant,0.6)
    w_std_med_infill_fun065 = weighted_med_std_infill(sm_interpolant,0.65)
    w_std_med_infill_fun07 = weighted_med_std_infill(sm_interpolant,0.7)
    w_std_med_infill_fun075 = weighted_med_std_infill(sm_interpolant,0.75)
    w_std_med_infill_fun08 = weighted_med_std_infill(sm_interpolant,0.8)
    w_std_med_infill_fun085 = weighted_med_std_infill(sm_interpolant,0.85)
    w_std_med_infill_fun09 = weighted_med_std_infill(sm_interpolant,0.9)
    w_std_med_infill_fun095 = weighted_med_std_infill(sm_interpolant,0.95)

    #Get the infill objective functions
    call(f, x) = f(x)
    library = Dict(
        :min => x -> min_infill_fun(x),
        :median => x -> median_infill_fun(x),
        :mean => x -> mean_infill_fun(x),
        :dist => x -> dist_infill_fun(x),
        :std => x -> std_infill_fun(x),
        :wstdmed005 => x -> w_std_med_infill_fun005(x),
        :wstdmed01 => x -> w_std_med_infill_fun01(x),
        :wstdmed015 => x -> w_std_med_infill_fun015(x),
        :wstdmed02 => x -> w_std_med_infill_fun02(x),
        :wstdmed025 => x -> w_std_med_infill_fun025(x),
        :wstdmed03 => x -> w_std_med_infill_fun03(x),
        :wstdmed035 => x -> w_std_med_infill_fun035(x),
        :wstdmed04 => x -> w_std_med_infill_fun04(x),
        :wstdmed045 => x -> w_std_med_infill_fun045(x),
        :wstdmed05 => x -> w_std_med_infill_fun05(x),
        :wstdmed055 => x -> w_std_med_infill_fun055(x),
        :wstdmed06 => x -> w_std_med_infill_fun06(x),
        :wstdmed065 => x -> w_std_med_infill_fun065(x),
        :wstdmed07 => x -> w_std_med_infill_fun07(x),
        :wstdmed075 => x -> w_std_med_infill_fun075(x),
        :wstdmed08 => x -> w_std_med_infill_fun08(x),
        :wstdmed085 => x -> w_std_med_infill_fun085(x),
        :wstdmed09 => x -> w_std_med_infill_fun09(x),
        :wstdmed095 => x -> w_std_med_infill_fun095(x)
    )
    functions_to_call = Tuple([library[s] for s in infill_funcs])
    infill_obj_fun = function (x)
        call.(functions_to_call, Ref(x))
    end
    return infill_obj_fun
end

"""
    closest_index(x_val, vals)

Based on the allowed categorical values, returns the index to the point closest
to `x_val`.
"""
function closest_index(x_val, vals)
    ibest = first(eachindex(vals))
    dxbest = abs(vals[ibest]-x_val)
    for I in eachindex(vals)
        dx = abs(vals[I]-x_val)
        if dx < dxbest
            dxbest = dx
            ibest = I
        end
    end
    ibest
end

"""
    colinmat(mat,vec)

Check if the vector `vec` exits as a column in `mat`. Greedy implementation
potentially not fit for high speed applications of large matrices.
"""
function colinmat(mat,vec)
    for column in eachcol(mat)
        all(column .== vec) && return true
    end
    false
end

"""
    infill_opt(search_range,infill_iterations,num_infill_points,infill_obj_fun,infill_funcs,plan,sm_interpolant,options)

Optimizes the infill points based on the supplied `infill_funcs` criteria.
Does not allow duplicated points.
"""
function infill_opt(search_range,infill_iterations,num_infill_points,infill_obj_fun,infill_funcs,plan,sm_interpolant,options)
    # Try generating pareto optimal infill points. Wrapped in try block due to
    # the algorithm sometimes failing to generate a result.
    infill_incomplete = true
    j = 0
    res_bboptim = nothing
    while infill_incomplete && j <= 100
        try
            if num_infill_points > 1
                #Multi-objective optimisation when number of infill points is larger than 1
                res_bboptim = bboptimize(infill_obj_fun; Method=:borg_moea,
                        FitnessScheme=ParetoFitnessScheme{length(infill_funcs[1:num_infill_points])}(is_minimizing=true),
                        SearchRange=search_range, ϵ=0.000001,
                        MaxFuncEvals=infill_iterations,
                        MaxStepsWithoutProgress=infill_iterations,TraceMode=:silent);
            else
                res_bboptim = bboptimize(x->infill_obj_fun(x)[1];
                        SearchRange=search_range,
                        MaxFuncEvals=infill_iterations,
                        MaxStepsWithoutProgress=infill_iterations,TraceMode=:silent);
            end
            infill_incomplete = false
        catch ex
            j += 1
            (j >= 100) && rethrow(ex)  #Show error if consistently failing
        end
    end

    # Pick the infill points
    infill_plan = Array{Float64,2}(undef,size(plan,1),0)
    infill_type = Array{Symbol,1}()
    infill_prediction = Array{Float64,1}()

    # Take the num_infill_points number of infill points from the cycled list.
    infill_obj_funs = Iterators.take(Base.Iterators.cycle(1:length(infill_funcs[1:num_infill_points])), num_infill_points)

    for i in infill_obj_funs
        # pf = pareto_frontier(res_bboptim)
        # best_obj1, idx_obj1 = findmin(map(elm -> fitness(elm)[i], pf))
        # bo1_solution = BlackBoxOptim.params(pf[idx_obj1]) # get the solution candidate itself...
        bo1_solution = best_candidate(res_bboptim)

        # Add the infill point if it does not exist in the plan or infill_plan
        v = copy(permutedims(bo1_solution'))
        if !colinmat(infill_plan,v) && !colinmat(plan,v)
            push!(infill_prediction,median(sm_interpolant(vec(v))))
            infill_plan = [infill_plan v]
            push!(infill_type,infill_funcs[i])
        end
    end

    # Cycle through the objective functions and update options with new list.
    ifuncs=circshift(infill_funcs,-num_infill_points)
    circshift_perm=[findfirst(isequal(x),ifuncs) for x in infill_funcs]
    options = Options(options; infill_funcs=ifuncs)

    return infill_plan,infill_type,infill_prediction,res_bboptim,options
end

"""
    infill_add(sm_interpolant,samples,plan,infill_prediction,search_range,infill_plan,infill_type,num_infill_points,infill_obj_fun,infill_funcs,res_bboptim)

Adds additional sample points if the `infill_opt` can't supply the number of
requested infill points. First a pareto optimal point based on the infill
functions is selected. The pareto point furthest away from all existing points
is selected. If this fails the last points are added as randomly selected points
to ensure that the requested number of infill points is added. Does not allow
duplicated points.
"""
function infill_add(sm_interpolant,samples,plan,infill_prediction,search_range,infill_plan,infill_type,num_infill_points,infill_obj_fun,infill_funcs,res_bboptim)

    try # Solution to handle the situation when pareto_frontier is not available
    for i = 1:(num_infill_points-size(infill_plan,2))
        infill_samples = Array{Float64,2}(undef,1,size(infill_plan,2))
        for i = 1:size(infill_plan,2)
            infill_samples[i] = median(sm_interpolant(infill_plan[:,i]))
        end

        # Find the point in the pareto front which is located furthest away from previous samples
        dist_infill_fun = distance_infill([plan infill_plan],[samples infill_samples],sm_interpolant)
        dist_obj = Inf
        par_f = pareto_frontier(res_bboptim)
        best_add_infill_solution = BlackBoxOptim.params(par_f[1])
        for pf in par_f
            bo1_solution = BlackBoxOptim.params(pf)
            if (cur_val = dist_infill_fun(bo1_solution)) < dist_obj
                dist_obj = cur_val
                best_add_infill_solution = copy(bo1_solution)
            end
        end

        # Add the infill point if it does not exist in the plan or infill_plan
        v = copy(permutedims(best_add_infill_solution'))
        if !colinmat(infill_plan,v) && !colinmat(plan,v)
            push!(infill_prediction,median(sm_interpolant(vec(v))))
            infill_plan = [infill_plan v]
            push!(infill_type,:pareto)
        end
    end
    catch
    end

    # Add random points if there aren't enough in the pareto front
    while (num_infill_points-size(infill_plan,2)) != 0

        v = Array{Float64,2}(undef,size(infill_plan,1),1)
        for i = 1:length(v)
            v[i] = scale(rand(),search_range[i][1],search_range[i][2],old_min=0,old_max=1)
        end
        # Add the infill point if it does not exist in the plan or infill_plan
        if !colinmat(infill_plan,v) && !colinmat(plan,v)
            push!(infill_prediction,median(sm_interpolant(vec(v))))
            infill_plan = [infill_plan v]
            push!(infill_type,:rand)
        end
    end

    return infill_plan,infill_type,infill_prediction
end

print_infill_head() = println("Finding new infill samples")

function print_infill_tail(infill_type,num_infill_points,infill_plan,infill_prediction)
    print("    ")
    for i = (length(infill_type)-num_infill_points+1):length(infill_type)
        print(@sprintf("%-15s",infill_type[i]))
    end
    print("\n")
    println("    ---------------------------------------------------------------")
    for j = 1:size(infill_plan,1)
        print("    ")
        for i = 1:size(infill_plan,2)
            print(@sprintf("%-15.7g",infill_plan[j,i]))
        end
        print("\n")
    end
    println("    ---------------------------------------------------------------")
    _, min_loc = findmin(infill_prediction)
    print("    ")
    for i = 1:size(infill_plan,2)
        if i == min_loc
            printstyled(@sprintf("%-15.7g",infill_prediction[i]); bold=true)
        else
            printstyled(@sprintf("%-15.7g",infill_prediction[i]))
        end
    end
    println("\t prediction\n")
end


text:
export absolute_degree, absolute_discriminant

################################################################################
#
#  Base field
#
################################################################################

@doc Markdown.doc"""
    base_field(L::NumField) -> NumField

Given a number field $L/K$ this function returns the base field $K$.
For absolute extensions this returns $\mathbf{Q}$.
"""
base_field(::NumField)

_base_ring(K::NumField) = base_field(K)

_base_ring(::FlintRationalField) = FlintQQ

################################################################################
#
#  Predicates
#
################################################################################

export isabsolute

@doc Markdown.doc"""
    isabsolute(L::NumField) -> Bool

Returns whether $L$ is an absolute extension, that is, whether the base field
of $L$ is $\mathbf{Q}$.
"""
isabsolute(::NumField)

isabsolute(::NumField) = false

isabsolute(::NumField{fmpq}) = true

@doc Markdown.doc"""
    elem_type(L::NumField) -> Type

Returns the type of the elements of $L$.
"""
elem_type(::NumField)

################################################################################
#
#  Degree
#
################################################################################

@doc Markdown.doc"""
    degree(L::NumField) -> Int

Given a number field $L/K$, this function returns the degree of $L$ over $K$.
"""
degree(A::NumField)

dim(K::NumField) = degree(K)

################################################################################
#
#  Absolute degree
#
################################################################################

@doc Markdown.doc"""
    absolute_degree(L::NumField) -> Int

Given a number field $L/K$, this function returns the degree of $L$ over
$\mathbf Q$.
"""
function absolute_degree(A::NumField)
  return absolute_degree(base_field(A)) * degree(A)
end

function absolute_degree(K::NumField{fmpq})
  return degree(K)
end

################################################################################
#
#  Is simple extension
#
################################################################################

@doc Markdown.doc"""
    issimple(L::NumField) -> Bool

Given a number field $L/K$ this function returns whether $L$ is simple, that is,
whether $L/K$ is defined by a univariate polynomial$.
"""
issimple(a::NumField)

################################################################################
#
#  Number field creation
#
################################################################################

@doc Markdown.doc"""
    NumberField(f::Poly{NumFieldElem}, s::String;
                cached::Bool = false, check::Bool = false) -> NumField, NumFieldElem

Given an irreducible polynomial $f \in K[x]$ over some number field $K$, this
function creates the simple number field $L = K[x]/(x)$ and returns $(L, b)$,
where $b$ is the class of $x$ in $L$. The string `s` is used only for printing
the primitive element $b$.

Testing that $f$ is irreducible can be disabled by setting the keyword argument
`check` to `false`.
"""
NumberField(f::PolyElem{<:NumFieldElem}, s::String;
            cached::Bool = false, check::Bool = false)

################################################################################
#
#  Is commutative
#
################################################################################

iscommutative(K::NumField) = true
iscommutative(::FlintRationalField) = true

################################################################################
#
#  Normal basis
#
################################################################################

@doc Markdown.doc"""
    normal_basis(L::NumField) -> NumFieldElem

Given a normal number field $L/K$, this function returns an element $a$ of $L$,
such that the orbit of $a$ under the Galois group of $L/K$ is an $K$-basis
of $L$.
"""
function normal_basis(L::NumField)
  n = degree(L)
  K = base_field(L)
  Aut = automorphisms(L)

  length(Aut) != degree(L) && error("The field is not normal over the rationals!")

  A = zero_matrix(K, n, n)
  r = one(L)
  while true
    r = rand(basis(L), -n:n)
    for i = 1:n
      y = Aut[i](r)
      for j = 1:n
        A[i,j] = coeff(y, j - 1)
      end
    end
    if rank(A) == n
      break
    end
  end
  return r
end

function isnormal_basis_generator(a::NumFieldElem)
  L = parent(a)
  n = degree(L)
  K = base_field(L)
  Aut = automorphisms(L)

  length(Aut) != degree(L) && error("The field is not normal over the rationals!")

  A = zero_matrix(K, n, n)
  for i = 1:n
    y = Aut[i](a)
    for j = 1:n
      A[i,j] = coeff(y, j - 1)
    end
  end
  return rank(A) == n
end

#trivia to make life easier

gens(L::SimpleNumField{T}) where {T} = [gen(L)]

function gen(L::SimpleNumField{T}, i::Int) where {T}
  i == 1 || error("index must be 1")
  return gen(L)
end

gen(L::NonSimpleNumField{T}, i::Int) where {T} = L(gen(parent(L.pol[1]), i))

function Base.getindex(L::SimpleNumField{T}, i::Int) where {T}
  if i == 0
    return one(L)
  elseif i == 1
    return gen(L)
  else
    error("index has to be 0 or 1")
  end
end

ngens(L::SimpleNumField{T}) where {T} = 1

function Base.getindex(L::NonSimpleNumField{T}, i::Int) where {T}
  i == 0 && return one(L)
  return gen(L, i)
end

function iscached(L::AnticNumberField)
  if haskey(Nemo.AnticNumberFieldID, (parent(L.pol), L.pol, L.S))
    return Nemo.AnticNumberFieldID[parent(L.pol), L.pol, L.S] === L
  end
  return false
end

function iscached(L::NfRel)
  if haskey(NfRelID, (parent(L.pol), L.pol, L.S))
    return NfRelID[parent(L.pol), L.pol, L.S] === L
  end
  return false
end

iscached(L::NonSimpleNumField) = false

export set_var!, set_vars!

#the Symbol is part of the key for caching, hence it should be be changed
@doc Markdown.doc"""
    set_var!(L::SimpleNumField, s::String)
    set_var!(L::SimpleNumField, s::Symbol)

Sets the name used when printing the primitive element of $L$.
This can only be set on fields constructed using `cached = false`.
"""
function set_var!(L::SimpleNumField{T}, s::String) where {T}
  iscached(L) && error("cannot set the name in a cached field")
  L.S = Symbol(s)
  nothing
end

function set_var!(L::SimpleNumField{T}, s::Symbol) where {T}
  iscached(L) && error("cannot set the name in a cached field")
  L.S = s
  nothing
end
@doc Markdown.doc"""
    set_vars!(L::NonSimpleNumField{T}, a::String)
    set_vars!(L::NonSimpleNumField{T}, a::Symbol)

Sets the string printed for each generator of the field. If the
string contains '#', then the hash-character is replaced by the index,
otherwise, the index is appended to the string.
Eg. `set_vars!(L, "g[#]")` will make the generators print like
array elements.
"""
function set_vars!(L::NonSimpleNumField{T}, a::Symbol) where {T}
  return set_vars!(L, String(a))
end

function set_vars!(L::NonSimpleNumField{T}, a::String) where {T}
  n = ngens(L)
  if occursin('#', a)
    S = [replace(a, '#'=>"$i") for i=1:n]
  else
    S = ["$a$i" for i=1:n]
  end
  return set_vars!(L, S)
end

@doc Markdown.doc"""
    set_vars!(L::NonSimpleNumField{T}, a::Array{String, 1})
    set_vars!(L::NonSimpleNumField{T}, a::Array{Symbol, 1})

Set the printing names for the generators to the string specified in
the array. The length has to be exactly `ngens(L)`.
"""
function set_vars!(L::NonSimpleNumField{T}, a::Array{String, 1}) where {T}
  length(a) == ngens(L) || error("need to have as many strings as generators")
  L.S = [Symbol(s) for s = a]
  nothing
end

function set_vars!(L::NonSimpleNumField{T}, a::Array{Symbol, 1}) where {T}
  length(a) == ngens(L) || error("need to have as many strings as generators")
  L.S = a
  nothing
end

iscyclotomic_type(K::NonSimpleNumField{T}) where {T} = false, fmpz(1)
iscyclotomic_type(K::NfRel) = false, fmpz(1)
function iscyclotomic_type(L::AnticNumberField)
  f = get_special(L, :cyclo)
  if f === nothing
    return false, fmpz(1)
  end
  return true, f
end

isquadratic_type(K::NonSimpleNumField{T}) where {T} = false, fmpz(1)
isquadratic_type(K::NfRel) = false, fmpz(1)
function isquadratic_type(L::AnticNumberField)
  f = get_special(L, :show)
  if f === Hecke.show_quad
    return true, numerator(-coeff(L.pol, 0))
  end
  return false, fmpz(1)
end

################################################################################
#
#  Absolute basis
#
################################################################################
@doc Markdown.doc"""
    absolute_basis(K::NumField) -> Vector{NumFieldElem}

Returns an array of elements that form a basis of $K$ (as a vector space)
over the rationals.
"""
absolute_basis(::NumField)

function absolute_basis(K::NumField)
  k = base_field(K)
  kabs = absolute_basis(k)
  B = basis(K)
  res = Vector{elem_type(K)}(undef, absolute_degree(K))
  ind = 1
  for b in basis(K)
    for bb in kabs
      res[ind] = bb * b
      ind += 1
    end
  end
  return res
end

function absolute_basis(K::NumField{fmpq})
  return basis(K)
end

################################################################################
#
#  Discriminant sign
#
################################################################################

@doc Markdown.doc"""
    discriminant_sign(K::NumField) -> Int

Returns the sign of the discriminant of the maximal order of $K$.
"""
function discriminant_sign(K::NumField)
# Compute the sign using Brill's theorem
  _, s = signature(K)
  if mod(s, 2) == 0
    return 1
  else
    return -1
  end
end

ext: jl


text:
# Autogenerated wrapper script for MSTM_jll for powerpc64le-linux-gnu-libgfortran5
export mstm

using MPICH_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("MSTM")
JLLWrappers.@declare_executable_product(mstm)
function __init__()
    JLLWrappers.@generate_init_header(MPICH_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        mstm,
        "bin/mstm",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct KruskalHeapEntry{E, T<:Real}
    edge::E
    dist::T
end

isless(e1::KruskalHeapEntry, e2::KruskalHeapEntry) = e1.dist < e2.dist

""" Performs [Quick-Find algorithm](https://en.wikipedia.org/wiki/Disjoint-set_data_structure)
on a given pair of nodes `p`and `q`, and makes a connection between them
in the vector `nodes`.
"""
function quick_find!(nodes::Vector, p, q)
    pid = nodes[p]
    qid = nodes[q]
    for i in 1:length(nodes)
        if nodes[i] == pid
            nodes[i] = qid
        end
    end
end

"""
    minimum_spanning_tree{T<:Real}(
        g, distmx::AbstractMatrix{T} = weights(g)
    )

Performs [Kruskal's algorithm](https://en.wikipedia.org/wiki/Kruskal%27s_algorithm)
on a connected, undirected graph `g`, having adjacency matrix `distmx`,
and computes minimum spanning tree. Returns a `Vector{KruskalHeapEntry}`,
that contains the containing edges and its weights.
"""
function minimum_spanning_tree(
    g::AGraph,
    distmx::AbstractMatrix{T} = weights(g)
) where T<:Real
    E = edgetype(g)
    V = vertextype(g)
    edge_list = Vector{KruskalHeapEntry{E,T}}()
    mst = Vector{E}()
    connected_nodes = Vector{V}(1:nv(g))

    sizehint!(edge_list, ne(g))
    sizehint!(mst, ne(g))

    for e in edges(g)
        heappush!(edge_list, KruskalHeapEntry(e, distmx[src(e),dst(e)]))
    end

    while !isempty(edge_list) && length(mst) < nv(g) - 1
        heap_entry = heappop!(edge_list)
        v = src(heap_entry.edge)
        w = dst(heap_entry.edge)

        if connected_nodes[v] != connected_nodes[w]
            quick_find!(connected_nodes, v, w)
            push!(mst, heap_entry.edge)
        end
    end

    return mst
end


text:
struct EPDensity{T<:AbstractMatrix,S<:AbstractVector}
    d::T
    detJ::T
    energy::S
    pitch::S
end

function local_distribution(M::AxisymmetricEquilibrium, grid::OrbitGrid, f::Vector, r, z;
                            energy=1.0:80.0, pitch=-1.0:0.02:1.0, nearest=false,kwargs...)
    f3d = map(grid,f)
    nenergy = length(energy)
    npitch = length(pitch)
    if nearest
        inds = filter(x->grid.orbit_index[x] != 0, CartesianIndices(grid.orbit_index))
        data = hcat(([grid.energy[I[1]], grid.pitch[I[2]], grid.r[I[3]]] for I in inds)...)
        tree = KDTree(data)
    end

    d = zeros(length(energy),length(pitch))
    detJ = zeros(length(energy),length(pitch))
    for i=1:nenergy, j=1:npitch
        for rr in r, zz in z
            v, dJ = eprz_to_eprt(M, energy[i], pitch[j], rr, zz; adaptive=false, kwargs...)
            dJ == 0 && continue
            if nearest
                idxs, dists = knn(tree,v[1:3],1,false)
                ii,jj,kk = Tuple(inds[idxs[1]])
            else
                !(grid.energy[1] <= v[1] <= grid.energy[end]) && continue
                !(grid.pitch[1] <= v[2] <= grid.pitch[end]) && continue
                !(grid.r[1] <= v[3] <= grid.r[end]) && continue
                ii = argmin(abs.(v[1] .- grid.energy))
                jj = argmin(abs.(v[2] .- grid.pitch))
                kk = argmin(abs.(v[3] .- grid.r))
            end
            d[i,j] = d[i,j] + dJ*f3d[ii,jj,kk]
            detJ[i,j] = dJ
        end
    end
    d = d/(length(r)*length(z))
    return EPDensity(d,detJ,energy,pitch)
end

function local_distribution(M::AxisymmetricEquilibrium, OS::OrbitSystem, f, orbs, sigma, r, z;
                            Js=Vector{Vector{Float64}}[],
                            energy=range(1.0,80.0,length=25),
                            pitch=range(-0.99,0.99,length=25),
                            distributed=false, atol=1e-3,
                            covariance=:local, norms=S3(1.0,1.0,1.0), kwargs...)

    nenergy = length(energy)
    npitch = length(pitch)
    m = orbs[1].coordinate.m
    q = orbs[1].coordinate.q
    #lorbs = Array{Orbit}(undef,nenergy,npitch)
    #@showprogress for i=1:nenergy
    #    for j = 1:npitch
    #        lorbs[i,j] = get_orbit(M, GCParticle(energy[i],pitch[j],r,z,m,q); kwargs...)
    #    end
    #end
    #lorbs = reshape(lorbs,nenergy*npitch)
    lorbs = reshape([get_orbit(M, GCParticle(energy[i],pitch[j],r,z,m,q); kwargs...) for i=1:nenergy,j=1:npitch],nenergy*npitch)
    if distributed
        lJs = pmap(o->get_jacobian(M,o), lorbs, on_error = ex->zeros(2))
                   #batch_size=round(Int, nenergy*npitch/(5*nprocs())))
    else
        lJs = [get_jacobian(M, o) for o in lorbs]
    end

    if covariance == :local
        Si = get_covariance_matrix(M, lorbs, orbs, sigma, Js_1=lJs, Js_2=Js,
                                   distributed=distributed, atol=atol)
    else
        Si = get_global_covariance_matrix(lorbs, orbs, sigma, norms=norms)
    end

    f_ep = reshape(max.(Si*(OS.S_inv*f),0.0),nenergy,npitch)
    detJ = reshape([length(j) != 0 ? j[1] : 0.0 for j in lJs],nenergy,npitch)
    f_ep .= f_ep./detJ
    f_ep[detJ .== 0.0] .= 0.0
    w = reshape([o.class in (:lost, :incomplete, :unknown) for o in lorbs],nenergy, npitch)
    f_ep[w] .= 0.0

    return EPDensity(f_ep,1.0./detJ,energy,pitch)
end

struct RZDensity{T<:AbstractMatrix,S<:AbstractVector}
    d::T
    r::S
    z::S
end

function rz_profile(M::AxisymmetricEquilibrium, OS::OrbitSystem, f::Vector, orbs, sigma;
                            Js=Vector{Vector{Float64}}[],
                            energy=range(1.0,80.0,length=25),
                            pitch=range(-0.99,0.99,length=25),
                            r = range(extrema(M.r)...,length=25),
                            z = range(extrema(M.z)...,length=25),
                            distributed=false, atol=1e-3, domain_check= (xx,yy) -> true,
                            covariance=:local, norms=S3(1.0,1.0,1.0),
                            checkpoint=true, warmstart=false,file="rz_progress.jld2", kwargs...)

    nenergy = length(energy)
    npitch = length(pitch)
    nr = length(r)
    nz = length(z)
    inds = CartesianIndices((nr,nz))
    m = orbs[1].coordinate.m
    q = orbs[1].coordinate.q

    f_rz = zeros(nr,nz)

    if warmstart && isfile(file) && (filesize(file) != 0)
        @load file f_rz last_ind
    else
        last_ind = inds[1]
    end

    if checkpoint
        touch(file)
    end

    for I in inds
        (I != inds[1] && I < last_ind) && continue
        i = I[1]
        j = I[2]
        rr = r[i]
        zz = z[j]
        !(domain_check(rr,zz)) && continue

        lorbs = reshape([get_orbit(M, GCParticle(energy[k],pitch[l],rr,zz,m,q); kwargs...) for k=1:nenergy,l=1:npitch],nenergy*npitch)
        if distributed
            lJs = pmap(o->get_jacobian(M,o), lorbs, on_error = ex->zeros(2))
                       #batch_size=round(Int, nenergy*npitch/(5*nprocs())))
        else
            lJs = [get_jacobian(M, o) for o in lorbs]
        end

        if covariance == :local
            Si = get_covariance_matrix(M, lorbs, orbs, sigma, Js_1=lJs, Js_2=Js,
                                       distributed=distributed, atol=atol)
        else
            Si = get_global_covariance_matrix(lorbs, orbs, sigma, norms=norms)
        end

        f_ep = reshape(max.(Si*(OS.S_inv*f),0.0),nenergy,npitch)
        detJ = reshape([length(j) != 0 ? j[1] : 0.0 for j in lJs],nenergy,npitch)
        f_ep .= f_ep./detJ
        f_ep[detJ .== 0.0] .= 0.0
        w = reshape([o.class in (:lost, :incomplete, :unknown) for o in lorbs],nenergy, npitch)
        f_ep[w] .= 0.0
        f_rz[i,j] = sum(f_ep)*step(energy)*step(pitch)
        if checkpoint
            @save file f_rz last_ind
        end
        last_ind = I
    end
    return RZDensity(f_rz,r,z)
end

struct EPRZDensity{T<:AbstractArray,S<:AbstractVector}
    d::T
    energy::S
    pitch::S
    r::S
    z::S
end

function eprz_distribution(M::AxisymmetricEquilibrium, OS::OrbitSystem, f::Vector, orbs, sigma;
                            Js=Vector{Vector{Float64}}[],
                            energy=range(1.0,80.0,length=25),
                            pitch=range(-0.99,0.99,length=25),
                            r = range(extrema(M.r)...,length=25),
                            z = range(extrema(M.z)...,length=25),
                            distributed=false, atol=1e-3, domain_check= (xx,yy) -> true,
                            covariance=:local, norms=S3(1.0,1.0,1.0),
                            checkpoint=true, warmstart=false,file="eprz_progress.jld2", kwargs...)

    nenergy = length(energy)
    npitch = length(pitch)
    nr = length(r)
    nz = length(z)
    inds = CartesianIndices((nr,nz))
    m = orbs[1].coordinate.m
    q = orbs[1].coordinate.q

    f_eprz = zeros(nenergy,npitch,nr,nz)

    if warmstart && isfile(file) && (filesize(file) != 0)
        @load file f_eprz last_ind
    else
        last_ind = inds[1]
    end

    if checkpoint
        touch(file)
    end

    for I in inds
        (I != inds[1] && I < last_ind) && continue
        i = I[1]
        j = I[2]
        rr = r[i]
        zz = z[j]
        !(domain_check(rr,zz)) && continue

        lorbs = reshape([get_orbit(M, GCParticle(energy[k],pitch[l],rr,zz,m,q); kwargs...) for k=1:nenergy,l=1:npitch],nenergy*npitch)
        if distributed
            lJs = pmap(o->get_jacobian(M,o), lorbs, on_error = ex->zeros(2))
                       #batch_size=round(Int, nenergy*npitch/(5*nprocs())))
        else
            lJs = [get_jacobian(M, o) for o in lorbs]
        end

        if covariance == :local
            Si = get_covariance_matrix(M, lorbs, orbs, sigma, Js_1=lJs, Js_2=Js,
                                       distributed=distributed, atol=atol)
        else
            Si = get_global_covariance_matrix(lorbs, orbs, sigma, norms=norms)
        end

        f_ep = reshape(max.(Si*(OS.S_inv*f),0.0),nenergy,npitch)
        detJ = reshape([length(j) != 0 ? j[1] : 0.0 for j in lJs],nenergy,npitch)
        f_ep .= f_ep./detJ
        f_ep[detJ .== 0.0] .= 0.0
        w = reshape([o.class in (:lost, :incomplete, :unknown) for o in lorbs],nenergy, npitch)
        f_ep[w] .= 0.0
        f_eprz[:,:,i,j] .= f_ep
        if checkpoint
            @save file f_eprz last_ind
        end
        last_ind = I
    end

    return EPRZDensity(f_eprz,energy,pitch,r,z)
end


text:
import StochasticIntegrals.ItoSet

"""
    generate_random_path(dimensions::Integer, ticks::Integer;
                         syncronous::Bool = false,
                         twister::MersenneTwister = MersenneTwister(1),
                         vol_dist::Distribution = Uniform(0.1/sqrt(252 * 8 * 3600), 0.5/sqrt(252 * 8 * 3600)),
                         refresh_rate_dist::Distribution    = Uniform(0.5, 5.0),
                         time_period_per_unit::Dates.Period = Second(1),
                         micro_noise_dist::Distribution     = Uniform(vol_dist.a * sqrt(time_period_ratio(Minute(5), time_period_per_unit)), vol_dist.b * sqrt(time_period_ratio(Minute(5), time_period_per_unit))),
                         assets::Union{Vector,Missing}      = missing,
                         brownian_corr_matrix::Union{Hermitian,Missing} = missing,
                         vols::Union{Vector,Missing}        = missing)

Generate a random path of price updates with a specified number of dimensions and ticks. There are options for whether the data is syncronous or asyncronous, the volatility of the price
processes, the refresh rate on the (exponential) arrival times of price updates, the minimum and the maximum microstructure noises.

Note the defaults are chosen to reflect a highcap stock with annualised volatility between 10% and 50%.
The standard deviation of microstructure noise is of the same order of magnitude as
5 minutes standard deviation of return. `vol * sqrt(60*5)` if vol is in seconds.
Refreshed ticks every 0.5-5 seconds (in expectation).

### Inputs
* `dimensions` - The number of assets.
* `ticks` - The number of ticks to produce.
* `syncronous` - Should ticks be syncronous (for each asset) or asyncronous.
* `twister` - The MersenneTwister used for RNG.
* `vol_dist` - The distribution to draw volatilities from (only used if vols is missing).
* `refresh_rate_dist` - The distribution to draw refresh rates (exponential distribution rates) from.
* `time_period_per_unit` - What time period should the time column correspond to.
* `micro_noise_dist`  - The distribution to draw assetwise microstructure noise standard deviations are drawn from.
* `assets` - The names of the assets that you want to use. The length of this must be equal to the `dimensions` input.
* `brownian_corr_matrix` - The correlation matrix to use. This is sampled from the Inverse Wishart distribution if none is input.
* `vols` - The volatilities to use. These are sampled  from the uniform distribution between `min_noise_var` and `max_noise_var`.
### Returns
* A `SortedDataFrame` of tick data.
* A `CovarianceMatrix` representing the true data generation process used in making the tick data.
* A `Dict` of microstructure noise variances for each asset.
* A `Dict` of update rates for each asset.
"""
function generate_random_path(dimensions::Integer, ticks::Integer;
                              syncronous::Bool = false,
                              twister::MersenneTwister = MersenneTwister(1),
                              vol_dist::Distribution = Uniform(0.1/sqrt(252 * 8 * 3600), 0.5/sqrt(252 * 8 * 3600)),
                              refresh_rate_dist::Distribution    = Uniform(0.5, 5.0),
                              time_period_per_unit::Dates.Period = Second(1),
                              micro_noise_dist::Distribution     = Uniform(vol_dist.a * sqrt(time_period_ratio(Minute(5), time_period_per_unit)), vol_dist.b * sqrt(time_period_ratio(Minute(5), time_period_per_unit))),
                              assets::Union{Vector,Missing}      = missing,
                              brownian_corr_matrix::Union{Hermitian,Missing} = missing,
                              vols::Union{Vector,Missing}        = missing)
    if (ismissing(assets) == false) && (dimensions != length(assets))
        error("If you input asset names then the number of asset names must be of length equal to the dimensions input.")
    end

    if ismissing(brownian_corr_matrix)
        wish = InverseWishart(dimensions, Matrix(Float64.(I(dimensions))))
        brownian_corr_matrix, _ = cov2cor(Hermitian(rand(twister, wish)))
    end

    vols = ismissing(vols) ? rand(twister, vol_dist, dimensions) : vols
    assets = ismissing(assets) ? Symbol.(:asset_, 1:dimensions) : assets

    ito_integrals = Dict(assets .=> map(i -> ItoIntegral(assets[i], PE_Function(vols[i], 0.0, 0.0, 0)  ), 1:dimensions))
    ito_set_ = ItoSet(brownian_corr_matrix, assets, ito_integrals)

    covar = StochasticIntegrals.SimpleCovariance(ito_set_, 0.0, 1.0; calculate_inverse = false, calculate_determinant = false)
    stock_processes = Dict(assets .=> map(a -> ItoProcess(0.0, 0.0, PE_Function(0.00, 0.0, 0.0, 0), ito_integrals[a]), assets))
    update_rates = Dict(assets .=> Exponential.(rand(twister, refresh_rate_dist, dimensions)))
    microstructure_noise = Dict(assets .=> rand(twister, micro_noise_dist, dimensions) .^ 2)
    ts = syncronous ? StochasticIntegrals.make_ito_process_syncronous_time_series(stock_processes, covar, mean(a-> update_rates[a].θ, assets), Int(ceil(ticks/dimensions)); ito_twister = twister) : StochasticIntegrals.make_ito_process_non_syncronous_time_series(stock_processes, covar, update_rates, ticks; timing_twister = twister, ito_twister = twister)
    ts = SortedDataFrame(ts, :Time, :Name, :Value, time_period_per_unit)

    standard_normal_draws = rand(twister, Normal(), nrow(ts.df))

    normal_draws = standard_normal_draws .* map(a -> sqrt(microstructure_noise[a]), ts.df[:,ts.grouping])
    ts.df[:,ts.value] += normal_draws
    return ts, CovarianceMatrix(brownian_corr_matrix, vols, assets, time_period_per_unit), microstructure_noise, update_rates
end

"""
    ItoSet(covariance_matrix::CovarianceMatrix{<:Real})

Convert a `CovarianceMatrix` into an `ItoSet` from the StochasticIntegrals package.
This package can then be used to do things like generate draws from the Multivariate
Gaussian corresponding to the covariance matrix and other things.
### Inputs
* `covariance_matrix` - The `CovarianceMatrix` that you want to convert into an `StochasticIntegrals.ItoSet`
### Returns
* A `StochasticIntegrals.ItoSet` struct.
"""
function ItoSet(covariance_matrix::CovarianceMatrix{<:Real})
    itos = Dict{Symbol,ItoIntegral}()
    for i in 1:length(covariance_matrix.labels)
        lab = covariance_matrix.labels[i]
        voll = covariance_matrix.volatility[i]
        itos[lab] =  ItoIntegral(lab, voll)
    end
    ito_set_ = StochasticIntegrals.ItoSet(covariance_matrix.correlation , covariance_matrix.labels, itos)
    return ito_set_
end


text:
#=
    sources:

    Journal of the Optical Society of America

    Günter Wyszecki & G. H. Fiedler: New Color-Matching Ellipses, 1971

    Color matching data were obtained by the following observers:

    DLM: David L. MacAdam (after Brown & MacAdam)
    PGN: P.G. Nutting Jr.
    WRB: W. R. J. Brown (after Brown & MacAdam)
    GW: 42 years old observer
    AR: 27 years old observer
    GF: 32 years old observer
    B(original name 12B): Average ellipses of W. R. J. Brown’s group of 12 observers
    FMC: predicted ellipses by the Friele-MacAdam-Chickering formula optimized for best fit with PGN ellipses.
=#
#=

    All data are given in Yxy color space coordinates!

=#

# experiment luminance data in [cd/m²]
#  background | test field
const ELLIPSE_LUMINANCE_DATA =
    #PGN_bg, PGN_lum
   [24.0     48.0;
    #WSS_bg, WSS_lum
    24.0     48.0;
    #FMC_bg, FMC_lum
    24.0     48.0;
    #WRB_bg, WRB_lum
    0.0     999.0;
    #DLM
    0.0     999.0;
    #12B
    8.9     999.0;
    #GF
    6.0      12.0;
    #AR
    6.0      12.0;
    #GW
    6.0      12.0]

# MacAdam color matching ellipses obtained by observer
# P.G. Nutting Jr., which refer to 2° test colors.
# L0 = Y0 ≈ 48 cd/m² surrounded by a 42° white surround
# with a luminance of ≈ 24 cd/m².

# Ellipses refer to standard deviations of color matches
# Estimated by MacAdam to be ~1/3 of the corresponding JND

    #  x0  y0    θ(deg)    a        b
const PGN_ellipses =
   [0.160  0.057   62.5  0.85e-03  0.35e-03;
    0.187  0.118   77.0  2.20e-03  0.55e-03;
    0.253  0.125   55.5  2.50e-03  0.50e-03;
    0.150  0.680  105.0  9.60e-03  2.30e-03;
    0.131  0.521  112.5  4.70e-03  2.00e-03;
    0.212  0.550  100.0  5.80e-03  2.30e-03;
    0.258  0.450   92.0  5.00e-03  2.00e-03;
    0.152  0.365  110.0  3.80e-03  1.90e-03;
    0.280  0.385   75.5  4.00e-03  1.50e-03;
    0.380  0.498   70.0  4.40e-03  1.20e-03;
    0.160  0.200  104.0  2.10e-03  0.95e-03;
    0.228  0.250   72.0  3.10e-03  0.90e-03;
    0.305  0.323   58.0  2.30e-03  0.90e-03;
    0.385  0.303   65.5  3.80e-03  1.60e-03;
    0.472  0.399   51.0  3.20e-03  1.40e-03;
    0.527  0.350   20.0  2.60e-03  1.30e-03;
    0.475  0.300   28.5  2.90e-03  1.10e-03;
    0.510  0.236   29.5  2.40e-03  1.20e-03;
    0.596  0.283   13.0  2.60e-03  1.30e-03;
    0.344  0.284   60.0  2.30e-03  0.90e-03;
    0.390  0.237   47.0  2.50e-03  1.00e-03;
    0.441  0.198   34.5  2.80e-03  0.95e-03;
    0.278  0.223   57.5  2.40e-03  0.55e-03;
    0.300  0.163   54.0  2.90e-03  0.60e-03;
    0.365  0.153   40.0  3.60e-03  0.95e-03]

    #WSS ellipses are threshold loci based on a 50% chance of seeing a color difference.

    #  x0  y0    θ(deg)    a        b
const WSS_ellipses =
   [0.160  0.057   85.5   2.4e-03   1.4e-03;
    0.187  0.118   76.8   6.3e-03   2.9e-03;
    0.253  0.125   53.4   7.4e-03   2.8e-03;
    0.150  0.680  116.6  24.2e-03  11.4e-03;
    0.131  0.521  105.3  22.6e-03  10.3e-03;
    0.212  0.550  105.7  19.4e-03  10.9e-03;
    0.258  0.450   91.3  16.3e-03   9.9e-03;
    0.152  0.365   97.0  17.8e-03   7.6e-03;
    0.280  0.385   79.9  14.8e-03   8.8e-03;
    0.380  0.498  110.0  12.8e-03   9.6e-03;
    0.160  0.200   92.0  11.2e-03   4.5e-03;
    0.228  0.250   76.7  12.6e-03   5.9e-03;
    0.305  0.323   68.0  13.5e-03   7.4e-03;
    0.385  0.303   65.8  12.2e-03   8.6e-03;
    0.472  0.399   55.4   9.5e-03   8.0e-03;
    0.527  0.350   38.4   8.8e-03   6.1e-03;
    0.475  0.300   41.5  10.9e-03   5.3e-03;
    0.510  0.236   34.3  10.9e-03   3.2e-03;
    0.596  0.283   31.0   8.3e-03   3.5e-03;
    0.344  0.284   56.7  12.6e-03   6.2e-03;
    0.390  0.237   45.3  11.8e-03   4.5e-03;
    0.441  0.198   36.4  11.4e-03   2.9e-03;
    0.278  0.223   62.1  11.5e-03   5.1e-03;
    0.300  0.163   49.6   9.7e-03   3.4e-03;
    0.365  0.153   38.8  10.1e-03   2.5e-03]

    #  x0  y0    θ(deg)    a        b
const FMC_ellipses =
   [0.160  0.057   76.4  0.85e-03  0.35e-03;
    0.187  0.118   72.0  2.00e-03  0.53e-03;
    0.253  0.125   53.5  2.30e-03  0.49e-03;
    0.150  0.680  100.5  7.30e-03  3.10e-03;
    0.131  0.521  100.3  5.60e-03  2.60e-03;
    0.212  0.550   91.1  5.20e-03  2.40e-03;
    0.258  0.450   80.4  4.10e-03  1.80e-03;
    0.152  0.365   92.4  4.30e-03  1.80e-03;
    0.280  0.385   73.3  3.60e-03  1.40e-03;
    0.380  0.498   68.9  4.00e-03  1.70e-03;
    0.160  0.200   85.7  3.00e-03  1.00e-03;
    0.228  0.250   72.8  3.20e-03  0.99e-03;
    0.305  0.323   65.5  3.30e-03  1.10e-03;
    0.385  0.303   61.2  3.30e-03  1.30e-03;
    0.472  0.399   53.1  3.10e-03  1.40e-03;
    0.527  0.350   43.5  2.80e-03  1.40e-03;
    0.475  0.300   43.4  2.80e-03  1.20e-03;
    0.510  0.236   32.8  2.70e-03  1.20e-03;
    0.596  0.283   27.8  2.40e-03  1.50e-03;
    0.344  0.284   56.9  3.10e-03  0.97e-03;
    0.390  0.237   46.1  2.90e-03  0.94e-03;
    0.441  0.198   35.5  2.90e-03  1.00e-03;
    0.278  0.223   61.2  3.00e-03  0.79e-03;
    0.300  0.163   49.8  2.70e-03  0.64e-03;
    0.365  0.153   38.1  2.80e-03  0.77e-03]

#=  l0 = 0.2 * log10(Y0), with Y0 ≡ L0 in cd/m²
    Y0 = 10^(5 * l0) cd/m²

    based on ds² =  1
    Colors of varying lightness presented in a dark surround.
=#

    # x0     y0    l0     θ       a       b
const WRB_ellipses =
   [0.580  0.341  0.253  171  3.0e-03  1.6e-03;
    0.422  0.330  0.260   40  2.8e-03  1.3e-03;
    0.364  0.505  0.243   67  2.4e-03  1.1e-03;
    0.370  0.218  0.195   35  2.7e-03  1.0e-03;
    0.492  0.446  0.212   32  1.7e-03  1.3e-03;
    0.308  0.424  0.201   72  2.8e-03  0.8e-03;
    0.446  0.450  0.260   49  1.6e-03  1.1e-03;
    0.319  0.581  0.230   90  2.3e-03  1.5e-03;
    0.682  0.302  0.264  156  1.3e-03  0.6e-03;
    0.650  0.292  0.265    1  1.6e-03  0.8e-03;
    0.314  0.183  0.097   39  2.8e-03  0.8e-03;
    0.221  0.148  0.055   46  1.2e-03  0.4e-03;
    0.668  0.317  0.185  164  1.9e-03  1.1e-03;
    0.484  0.267  0.113   29  3.0e-03  1.1e-03;
    0.436  0.366  0.211   17  2.7e-03  1.1e-03;
    0.525  0.390  0.199    0  2.1e-03  1.5e-03;
    0.552  0.290  0.240   17  2.7e-03  1.1e-03;
    0.249  0.204  0.100   61  2.5e-03  1.0e-03;
    0.258  0.325  0.165   79  3.1e-03  1.3e-03;
    0.268  0.322  0.168   66  3.0e-03  0.9e-03;
    0.310  0.288  0.160   61  2.8e-03  1.1e-03;
    0.255  0.255  0.131   69  2.5e-03  0.9e-03;
    0.383  0.489  0.280   90  1.8e-03  1.2e-03;
    0.323  0.558  0.265   93  3.2e-03  1.4e-03;
    0.354  0.358  0.267   43  2.7e-03  1.1e-03;
    0.416  0.236  0.231   40  1.8e-03  1.1e-03;
    0.312  0.317  0.265  176  1.8e-03  0.9e-03;
    0.473  0.317  0.313   27  1.6e-03  1.0e-03;
    0.358  0.230  0.220   42  2.3e-03  0.8e-03;
    0.321  0.491  0.245   92  2.7e-03  1.3e-03;
    0.294  0.367  0.258   78  2.1e-03  1.1e-03;
    0.294  0.367  0.006   69  5.5e-03  1.6e-03;
    0.385  0.400  0.275   62  2.0e-03  0.9e-03;
    0.199  0.524  0.139   97  5.1e-03  2.3e-03;
    0.212  0.642  0.135   98  8.3e-03  1.8e-03;
    0.203  0.419  0.146  107  3.9e-03  1.1e-03;
    0.171  0.189  0.188   92  2.3e-03  0.6e-03;
    0.191  0.329  0.152   86  3.6e-03  1.0e-03]

    # x0    y0     l0     θ      a         b
const DLM_ellipses =
   [0.580  0.341  0.253  159  3.0e-03  1.5e-03;
    0.422  0.330  0.260   25  2.9e-03  1.4e-03;
    0.364  0.505  0.243   58  2.8e-03  1.8e-03;
    0.370  0.218  0.195   24  2.8e-03  1.0e-03;
    0.492  0.446  0.212   18  2.0e-03  1.0e-03;
    0.308  0.424  0.201   56  2.3e-03  0.9e-03;
    0.446  0.450  0.260    3  1.7e-03  1.0e-03;
    0.319  0.581  0.230    0    1e-08    1e-08;
    0.682  0.302  0.264  150  1.3e-03  0.8e-03;
    0.650  0.292  0.265  174  1.7e-03  0.8e-03;
    0.314  0.183  0.097   32  2.5e-03  0.7e-03;
    0.221  0.148  0.055   46  1.2e-03  0.4e-03;
    0.668  0.317  0.185  162  2.0e-03  0.9e-03;
    0.484  0.267  0.113   13  2.7e-03  1.5e-03;
    0.436  0.366  0.211   22  1.3e-03  0.7e-03;
    0.525  0.390  0.199  168  2.3e-03  1.0e-03;
    0.552  0.290  0.240   16  2.6e-03  1.1e-03;
    0.249  0.204  0.100   56  2.3e-03  0.8e-03;
    0.258  0.325  0.165   56  2.1e-03  1.1e-03;
    0.268  0.322  0.168   59  1.8e-03  0.9e-03;
    0.310  0.288  0.160   51  2.0e-03  0.8e-03;
    0.255  0.255  0.131   60  2.4e-03  0.9e-03;
    0.383  0.489  0.280   78  2.0e-03  1.5e-03;
    0.323  0.558  0.265   83  2.4e-03  1.4e-03;
    0.354  0.358  0.267   43  2.1e-03  0.8e-03;
    0.416  0.236  0.231   35  3.2e-03  1.4e-03;
    0.312  0.317  0.265  143  2.1e-03  1.0e-03;
    0.473  0.317  0.313   18  2.3e-03  1.1e-03;
    0.358  0.230  0.220   35  3.0e-03  1.0e-03;
    0.321  0.491  0.245   71  2.2e-03  1.1e-03;
    0.294  0.367  0.258   61  2.2e-03  0.8e-03;
    0.294  0.367  0.006   69  3.5e-03  1.2e-03;
    0.385  0.400  0.275   38  1.6e-03  1.3e-03;
    0.199  0.524  0.139   90  3.0e-03  1.3e-03;
    0.212  0.642  0.135   83  5.2e-03  2.3e-03;
    0.203  0.419  0.146  106  4.2e-03  1.6e-03;
    0.171  0.189  0.188   89  2.3e-03  0.4e-03;
    0.191  0.329  0.152  100  2.3e-03  1.1e-03]

#The color-matching ellipses are cross sections through the centers of (weighted) mean color-matching ellipsoids determined by a group of 12 observers for 10° test colors of different luminances (l0 = 0.2 log10 Y0), with Y0 = L0 in cd/m²) presented in a large near-white surround with a luminance of 8.9 cd/m² (that is, l0 = 0.190). For details see Ref. 15 or 3. Test color (13) is identical to test color 13 except that its angular subtense is reduced to 2°. Test colors marked with a prime, e.g. 1', 2', etc. are within the chromaticity gamut of Fig. 2 of this paper.

    # x0     y0    l0     Θ    a    b
const B_ellipses =
   [0.308  0.280  0.240   67  2.1e-03  1.1e-03;
    0.339  0.493  0.268   88  3.7e-03  1.6e-03;
    0.426  0.219  0.238   43  2.3e-03  1.3e-03;
    0.536  0.411  0.264  169  2.9e-03  1.4e-03;
    0.457  0.368  0.248    4  1.9e-03  1.0e-03;
    0.288  0.630  0.247  104  5.4e-03  2.0e-03;
    0.254  0.456  0.236  112  1.8e-03  1.0e-03;
    0.674  0.308  0.252  169  1.8e-03  0.7e-03;
    0.209  0.137  0.202   77  2.3e-03  0.7e-03;
    0.340  0.178  0.211   65  2.3e-03  1.1e-03;
    0.422  0.483  0.236   86  3.4e-03  2.1e-03;
    0.262  0.343  0.241   47  1.8e-03  1.5e-03;
    0.368  0.352  0.267  165  1.9e-03  1.0e-03;
    0.643  0.304  0.230  175  2.4e-03  0.9e-03;
    0.385  0.546  0.241  128  4.3e-03  2.9e-03;
    0.270  0.549  0.247  127  4.5e-03  2.0e-03;
    0.308  0.307  0.263   35  1.8e-03  1.4e-03;
    0.562  0.252  0.227   14  2.0e-03  0.9e-03;
    0.217  0.606  0.255  128  6.0e-03  2.4e-03;
    0.196  0.293  0.245    6  2.4e-03  1.5e-03;
    0.161  0.410  0.248  158  3.0e-03  1.6e-03;
    0.257  0.202  0.234   40  2.8e-03  1.0e-03]

#=
    GF, WR, GW ellipses

    ds² = 7.81 chosen according to χ² for 3 deg. of freedom.
    On average, ellipsoids contain 95% of any random set
    of color matches.

=#

    # x0      y0    g11    g12   g22    g23    g33   g13    Θ      a       b
const GF_ellipses =
   [0.4786  0.4327   492   -275   717    114   459    39   34  5.0e-03  2.9e-03;
    0.3909  0.4375   690   -678   999   -568   853   270   39  7.2e-03  2.2e-03;
    0.4932  0.3462   146    123  2637  -1151  1678  -215  177  7.5e-03  1.7e-03;
    0.3079  0.3595  1999   -572  1014   -474  1372  -428   65  3.2e-03  1.9e-03;
    0.4458  0.4809   567   -355   372     74   542    37   53  8.8e-03  3.0e-03;
    0.4382  0.3417   379    109  2250   -531  1391  -368  177  4.6e-03  1.9e-03;
    0.3822  0.3372   698   -409  1634   -162  1142  -587   21  3.8e-03  2.1e-03;
    0.2295  0.2987  1017   -556  1447   -331   991   109   34  3.5e-03  2.1e-03;
    0.3310  0.2749  1292   -645   844   -526  1092   402   55  4.5e-03  2.1e-03;
    0.2766  0.2551  1928  -1350  2262   -651  1516   280   41  3.3e-03  1.5e-03;
    0.2054  0.2544  1134   -499  1372   -592  1369  -179   38  3.2e-03  2.1e-03;
    0.5868  0.3360   102     56   901   -196   416   -24  176  8.9e-03  2.9e-03;
    0.3545  0.4518  1162   -411   306   -144   656   253   68  7.6e-03  2.4e-03;
    0.2618  0.2912  3006  -1389  1085   -156   578    50   62  4.7e-03  1.4e-03;
    0.3101  0.3163  2118   -910  1280     34   949   -99   59  3.3e-03  1.7e-03;
    0.3251  0.3325  3514   -312  1674   -494  1722   593   81  2.2e-03  1.5e-03;
    0.2959  0.3009  2572   -510  1560   -812  1947  -124   67  2.4e-03  1.7e-03;
    0.3144  0.3355  4825  -1690  3441   -232  1746  1786   56  1.8e-03  1.1e-03;
    0.3425  0.5280   525   -132   131    -79  1426   207   73  9.3e-03  3.7e-03;
    0.3022  0.4725  1946   -489   296     -4   781    14   75  6.9e-03  1.9e-03;
    0.4873  0.3994   408   -108  1422   -246  1942  -294    6  4.4e-03  2.3e-03;
    0.5428  0.3776   223    -29  1107   -151  1981   -75    2  5.9e-03  2.7e-03;
    0.2670  0.4240   991   -238   330    214  1239  -212   72  5.6e-03  2.7e-03;
    0.3753  0.3867  1285   -868  1007   -162  1016  -134   50  5.4e-03  2.0e-03;
    0.5242  0.3140   157    171  1008   -244   513   -14  169  8.0e-03  2.7e-03;
    0.2730  0.3701  1490   -654   951    311   798  -392   56  3.9e-03  2.0e-03;
    0.4650  0.3122   265    142   950   -275   556  -221  169  5.7e-03  2.8e-03;
    0.4058  0.2909   425   -174  1536   -407   926  -286    9  4.4e-03  2.2e-03]

    # x0      y0    g11    g12   g22   g23    g33   g13   Θ      a       b
    const AR_ellipses=
   [0.4786  0.4327   345   -389   866   152   483  -149  28  7.5e-03  2.7e-03;
    0.3909  0.4375   646   -344   479   -53   269     6  52  6.1e-03  2.9e-03;
    0.4932  0.3462  1105  -1088  2251    41   635  -147  31  2.2e-03  1.6e-03;
    0.3079  0.3595  2454   -788   933  -100   668   516  67  3.6e-03  1.7e-03;
    0.4458  0.4809   636   -391   515    -7   590   102  49  6.6e-03  2.8e-03;
    0.4382  0.3417  1082   -329  1573  -335  1226   456  27  2.9e-03  2.1e-03;
    0.3822  0.3372   771   -696  1442   273   986  -283  32  4.8e-03  2.0e-03;
    0.2295  0.2987  1684    652   958  -198  1129   135  60  3.7e-03  1.9e-03;
    0.3310  0.2749  1874  -1127  2103   -51  1026   358  42  3.0e-03  1.6e-03;
    0.2766  0.2551  4476  -1596  1181   124   577   -14  68  3.8e-03  1.2e-03;
    0.2054  0.2544  1448   -429   485  -109  1106     3  66  4.0e-03  2.2e-03;
    0.5868  0.3360   325    -14  1324   260   444   133   1  4.9e-03  2.4e-03;
    0.3545  0.4518  2537   -814   587  -163   721   488  70  5.2e-03  1.7e-03;
    0.2618  0.2912  2368  -1013  1011  -194   946   319  62  4.1e-03  1.6e-03;
    0.3101  0.3163  3291  -1049  1598    35   888   336  64  2.7e-03  1.4e-03;
    0.3251  0.3325  2718   -831  1565   220  1102   474  62  2.6e-03  1.6e-03;
    0.2959  0.3009  2387   -686  1110   142   815   -97  66  3.1e-03  1.7e-03;
    0.3144  0.3355  1639   -303   853   -85   708   196  71  3.2e-03  2.1e-03;
    0.3425  0.5280  1566   -242   182   -58   630   360  80  7.4e-03  2.2e-03;
    0.3022  0.4725  1835    -14   180   102  1025   761  90  6.6e-03  2.1e-03;
    0.4873  0.3994   516   -279   911   310   555   -84  27  4.6e-03  2.7e-03;
    0.5428  0.3776   335    -49  1130   341   447    24   4  4.8e-03  2.6e-03;
    0.2670  0.4240  1723   -169   206   -13   683   291  84  6.4e-03  2.1e-03;
    0.3753  0.3867  1160   -596   917  -312  1058    17  51  4.3e-03  2.2e-03;
    0.5242  0.3140   241    -94   692   178   587    63  11  5.9e-03  3.3e-03;
    0.2730  0.3701  1836   -459   577   -63   464    -3  72  4.3e-03  2.0e-03;
    0.4650  0.3122   383   -143   788  -107   677   -95  18  4.8e-03  3.1e-03;
    0.4058  0.2909   613   -447   815   -73   176   198  39  5.5e-03  2.6e-03]

    # x0      y0    g11    g12   g22   g23    g33   g13   Θ      a       b
    const GW_ellipses=
   [0.4786  0.43273  674   -695  1176    34   773    -6  35  5.5e-03  2.2e-03;
    0.3909  0.43753 1492  -1042  1251  -104   577   252  48  4.9e-03  1.8e-03;
    0.4932  0.34623  118     30  1491   -46  1501   136 179  8.2e-03  2.3e-03;
    0.3079  0.35953 1938   -629  1054  -268   751   382  63  3.3e-03  1.9e-03;
    0.4458  0.48093  692   -393   540    -1   668   144  50  6.0e-03  2.8e-03;
    0.4382  0.34173  532    -75  1328   113   689    56   5  3.9e-03  2.4e-03;
    0.3822  0.33723 1185  -1092  2074    -1   816   -67  34  4.2e-03  1.7e-03;
    0.2295  0.29873 1628   -508  1089  -741  3695    61  59  3.2e-03  2.0e-03;
    0.3310  0.27493 1512   -859   936  -389   805    79  54  5.0e-03  1.9e-03;
    0.2766  0.25513 6793  -2841  1855  -593  1898  1591  65  3.7e-03  1.0e-03;
    0.2054  0.25443 1242   -249   482   -75   954   242  73  4.4e-03  2.4e-03;
    0.5868  0.33603   78     10  1025   184   454   -59 179 10.0e-03  2.8e-03;
    0.3545  0.45183 1342   -374   284   -40   529  -136  72  6.9e-03  2.3e-03;
    0.2618  0.29123 4691  -1527  1068  -127   295   211  70  3.9e-03  1.2e-03;
    0.3101  0.31633 3406  -1554  1641  -691  1930   739  60  3.3e-03  1.4e-03;
    0.3251  0.33253 2717  -1656  2146  -141   901   572  50  3.2e-03  1.4e-03;
    0.2959  0.30093 2365  -1215  1261  -499  1683   358  57  4.0e-03  1.6e-03;
    0.3144  0.33553 3426  -1759  1525   128   994   306  59  4.0e-03  1.3e-03;
    0.3425  0.52803  769   -181    85   -37   875   247  76 14.0e-03  3.1e-03;
    0.3022  0.47253 2192   -501   225   -52   717   487  77  8.6e-03  1.8e-03;
    0.4873  0.39943  223   -331   959   211  1034  -229  21  9.0e-03  2.7e-03;
    0.5428  0.37763  119    -65   573  -191   968  -108   8  8.4e-03  3.7e-03;
    0.2670  0.42403  965   -114   186    98  1546   366  82  6.8e-03  2.8e-03;
    0.3753  0.38673  857   -874  1396    19  1165   137  36  6.1e-03  2.0e-03;
    0.5242  0.31403  136    -58   559    -6   280   -88   8  7.8e-03  3.7e-03;
    0.2730  0.37013 1251    178   511   -71   875    70 103  4.1e-03  2.5e-03;
    0.4650  0.31223  177   -122   653   228   755   -69  14  7.3e-03  3.4e-03;
    0.4058  0.29093  441   -167   504  -199   794    94  40  5.1e-03  3.5e-03]


text:
using CliOptions

# An utility which executes a command line for many times
spec = CliOptionSpec(
    Option(Int, "-n", "--times"; default = 1,
           help = "Number of times to execute the command"),
    Positional("arguments"; multiple = true,
               help = "Command line arguments to execute"),
    HelpOption(),
    use_double_dash = true,
)
options = parse_args(spec)
for _ in 1:options.times
    cmd = Cmd(String[a for a in options.arguments])
    run(cmd)
end


text:
import AlphaZero.GI

using Crayons
using StaticArrays

const NUM_HOUSES_PER_PLAYER = 6

const NUM_INITIAL_SEEDS_PER_HOUSE = 3

const NUM_HOUSES = 2 * NUM_HOUSES_PER_PLAYER

const NUM_SEEDS =  NUM_HOUSES * NUM_INITIAL_SEEDS_PER_HOUSE

const Player = Int
const WHITE = 1
const BLACK = 2

other(p::Player) = 3 - p

struct Board
  stores :: SVector{2, UInt8}
  houses :: SMatrix{2, NUM_HOUSES_PER_PLAYER, UInt8, NUM_HOUSES}
end

const INITIAL_BOARD =
  let stores = SVector(0, 0)
  let houses = @SMatrix[
      NUM_INITIAL_SEEDS_PER_HOUSE
      for i in 1:2, j in 1:NUM_HOUSES_PER_PLAYER]
    Board(stores, houses)
  end end

const INITIAL_STATE = (board=INITIAL_BOARD, curplayer=WHITE)

struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  board :: Board
  curplayer :: Player
  finished :: Bool
end

function GI.init(::GameSpec)
  board = INITIAL_STATE.board
  curplayer = INITIAL_STATE.curplayer
  finished = false
  return GameEnv(board, curplayer, finished)
end

GI.spec(::GameEnv) = GameSpec()

GI.two_players(::GameSpec) = true

function GI.set_state!(g::GameEnv, state)
  g.board = state.board
  g.curplayer = state.curplayer
  if sum_houses(g.board, g.curplayer) == 0 || sum_houses(g.board, other(g.curplayer)) == 0
    g.finished = true
  end
end

const ACTIONS = collect(1:NUM_HOUSES_PER_PLAYER)
GI.actions(::GameSpec) = ACTIONS

#####
##### Position system to ease board manipulation
#####

struct HousePos
  player :: Int
  num :: Int
end

struct StorePos
  player :: Int
end

const Pos = Union{HousePos, StorePos}

function next_pos(pos::HousePos, player)
  if pos.num > 1
    return HousePos(pos.player, pos.num - 1)
  elseif pos.num == 1
    if pos.player == player
      return StorePos(player)
    else
      return HousePos(player, NUM_HOUSES_PER_PLAYER)
    end
  else
    @assert false
  end
end

function next_pos(pos::StorePos, player)
  @assert pos.player == player
  return HousePos(other(player), NUM_HOUSES_PER_PLAYER)
end

read_pos(b::Board, pos::HousePos) = b.houses[pos.player, pos.num]
read_pos(b::Board, pos::StorePos) = b.stores[pos.player]

function write_pos(b::Board, pos::HousePos, v)
  houses = setindex(b.houses, v, pos.player, pos.num)
  return Board(b.stores, houses)
end

function write_pos(b::Board, pos::StorePos, v)
  stores = setindex(b.stores, v, pos.player)
  return Board(stores, b.houses)
end

function opposite_pos(pos::HousePos)
  return HousePos(other(pos.player), NUM_HOUSES_PER_PLAYER - pos.num + 1)
end

#####
##### Defining game rules
#####

GI.actions_mask(g::GameEnv, player) = map(>(0), g.board.houses[player,:])

GI.actions_mask(g::GameEnv) = GI.actions_mask(g, g.curplayer)

function capture_last_and_opposite(b::Board, pos::HousePos)
  # @assert read_pos(b ,pos) == 1
  new_store = b.stores[pos.player] + read_pos(b, opposite_pos(pos)) + 1
  stores = setindex(b.stores, new_store, pos.player)
  b = write_pos(b, pos, 0)
  b = write_pos(b, opposite_pos(pos), 0)
  return Board(stores, b.houses)
end

sum_houses(b::Board, player) = sum(b.houses[player,:])

function capture_leftovers(b::Board, player)
  # @assert sum_houses(b, other(player)) == 0
  new_store = b.stores[player] + sum_houses(b, player)
  stores = setindex(b.stores, new_store, player)
  return Board(stores, zero(typeof(b.houses)))
end

# TODO tide up ending conditions
function GI.play!(g::GameEnv, a)
  pos = HousePos(g.curplayer, a)
  nseeds = read_pos(g.board, pos)
  @assert nseeds > 0
  g.board = write_pos(g.board, pos, 0)
  for i in 1:nseeds
    pos = next_pos(pos, g.curplayer)
    g.board = write_pos(g.board, pos, read_pos(g.board, pos) + 1)
  end
  # Check endgame
  if sum_houses(g.board, g.curplayer) == 0
    g.board = capture_leftovers(g.board, other(g.curplayer))
    g.finished = true
  elseif isa(pos, HousePos)
    # Capture opposite house if last seed was put in empty house on your side # MUST CHECK ENDGAME AFTER
    if read_pos(g.board, pos) == 1 && g.curplayer == pos.player
      g.board = capture_last_and_opposite(g.board, pos)
      # if captured last pieces of opposite player then game ends
      if sum_houses(g.board, other(g.curplayer)) == 0
        g.board = capture_leftovers(g.board, g.curplayer)
        g.finished = true
        return
      end
      if sum_houses(g.board, g.curplayer) == 0
        g.board = capture_leftovers(g.board, other(g.curplayer))
        g.finished = true
        return
      end
    end
    # Free turn if last seed was put in a store
    g.curplayer = other(g.curplayer)
  end
  return
end

GI.current_state(g::GameEnv) = (board=g.board, curplayer=g.curplayer)

GI.white_playing(g::GameEnv) = g.curplayer == WHITE

#####
##### Reward shaping
#####

function GI.game_terminated(g::GameEnv)
  return g.finished
end

function zero_one_reward(nw, nb)
  nw > nb && return 1.
  nw < nb && return -1.
  return 0.
end

linear_reward(nw, nb) = (nw - nb) / NUM_SEEDS

function GI.white_reward(g::GameEnv)
  if g.finished
    nw, nb = g.board.stores
    return zero_one_reward(nw, nb) # + linear_reward(nw, nb)
  else
    return 0.
  end
end

#####
##### Simple heuristic for minmax
#####

# This is an extremely naive heuristic and one can probably do better.
function GI.heuristic_value(g::GameEnv)
  nw, nb = g.board.stores
  v = nw - nb
  g.curplayer == BLACK && (v = -v)
  return Float64(v)
end

#####
##### ML interface
#####

function flip_colors(board)
  stores = @SVector [INITIAL_BOARD.stores[2],INITIAL_BOARD.stores[1]]
  houses = SMatrix{2, NUM_HOUSES_PER_PLAYER, UInt8, NUM_HOUSES}(
    [INITIAL_BOARD.houses[2,:]'; INITIAL_BOARD.houses[1,:]'])
  return Board(stores, houses)
end

function GI.vectorize_state(::GameSpec, state)
  board = state.curplayer == WHITE ? state.board : flip_colors(state.board)
  function cell(pos, chan)
    if chan == :nstones
      return read_pos(board, pos)
    elseif chan == :whouse
      return isa(pos, HousePos) && pos.player == WHITE
    elseif chan == :wstore
      return isa(pos, StorePos) && pos.player == WHITE
    elseif chan == :bhouse
      return isa(pos, HousePos) && pos.player == BLACK
    elseif chan == :bstore
      return isa(pos, StorePos) && pos.player == BLACK
    end
    @assert false
  end
  positions = [
    [HousePos(WHITE, i) for i in NUM_HOUSES_PER_PLAYER:-1:1];
    [StorePos(WHITE)];
    [HousePos(BLACK, i) for i in NUM_HOUSES_PER_PLAYER:-1:1];
    [StorePos(BLACK)]]
  return Float32[
    cell(p, c)
    for p in positions,
        y in 1:1,
        c in [:nstones, :whouse, :wstore, :bhouse, :bstore]]
end

#####
##### User interface
#####

GI.action_string(::GameSpec, a) = string(a)

function GI.parse_action(::GameSpec, str)
  try
    p = parse(Int, str)
    1 <= p <= NUM_HOUSES_PER_PLAYER ? p : nothing
  catch
    nothing
  end
end

#          1    2    3    4    5    6
#  +----+----+----+----+----+----+----+----+
#  |    |  3 |  3 |  3 |  3 |  3 |  3 |    |
#  |    +----+----+----+----+----+----+    |
#  |    |  3 |  3 |  3 |  3 |  3 |  3 |    |
#  +----+----+----+----+----+----+----+----+
#          6    5    4    3    2    1

function GI.render(g::GameEnv, with_position_names=true, botmargin=true)
  b = g.board
  gray(s...) = print(crayon"cyan", s..., crayon"reset")
  white(s...) = print(crayon"white", s..., crayon"reset")
  bold(s...) = print(crayon"yellow", s..., crayon"reset")
  blank_cell() = repeat(" ", 4)
  num_cell(n) = n == 0 ? blank_cell() : (n < 10 ? "  $n " : " $n ")
  function show_labels(labels)
    print(" "); print(blank_cell()); print(" ")
    for i in labels
      gray(num_cell(i)); print(" ")
    end
    print("\n")
  end
  function hline(len; newline=true)
    for i in 1:len
      gray("+----")
    end
    newline && gray("+\n")
  end
  function print_houses(houses, self)
    printing = self ? bold : white
    for v in houses
      gray("|"); printing(num_cell(v))
    end
  end
  function print_half(houses, self=false)
    gray("|"); print(blank_cell())
    print_houses(houses, self)
    gray("|"); print(blank_cell())
    gray("|\n")
  end
  function print_middle(stores)
    gray("|"); white(num_cell(stores[2]))
    hline(NUM_HOUSES_PER_PLAYER, newline=false)
    gray("+"); white(num_cell(stores[1]))
    gray("|\n")
  end
  show_labels(1:NUM_HOUSES_PER_PLAYER)
  hline(NUM_HOUSES_PER_PLAYER + 2)
  print_half(b.houses[2,:], g.curplayer == BLACK)
  print_middle(b.stores)
  print_half(reverse(b.houses[1,:]), g.curplayer == WHITE)
  hline(NUM_HOUSES_PER_PLAYER + 2)
  show_labels(reverse(1:NUM_HOUSES_PER_PLAYER))
  botmargin && print("\n")
end

function GI.read_state(::Type{GameEnv})
  try
    function read_houses(player)
      print("Player $(player) houses: ")
      houses = [parse(Int, s) for s in split(readline())]
      @assert all(>=(0), houses)
      @assert length(houses) == NUM_HOUSES_PER_PLAYER
      return houses
    end
    function read_store(player)
      print("Player $(player) store: ")
      store = parse(Int, readline())
      @assert store >= 0
      return store
    end
    # Read board
    h1 = read_houses(1)
    h2 = read_houses(2)
    s1 = read_store(1)
    #s2 = read_store(2)
    nsofar = sum(h1) + sum(h2) + s1
    @assert nsofar <= NUM_SEEDS
    s2 = NUM_SEEDS - nsofar
    board = Board([s1, s2], [h1'; h2'])
    # Read current player
    print("Current player (1/2): ")
    curplayer = parse(Int, readline())
    @assert 1 <= curplayer <= 2
    return GameEnv(board, curplayer == WHITE)
  catch
    return nothing
  end
end

# Example of creating game objects manually:
# Game(Board([10, 20], [[1 0 0 0 0 0]; [1 0 2 0 2 0]]), true)


text:
using GeometricTheoremProver
using OrderedCollections
using DynamicPolynomials
using Test
const GTP = GeometricTheoremProver

using GeometricTheoremProver: Segment, Parallelism, Perpendicularity, Congruence, Appartenence, Intersection, Midpoint, Parallelogram, Projection

include("test_parser.jl")
include("test_rittwu.jl")


text:
  using Random, Plots, LaTeXStrings, ProgressMeter
let

  # initialization
  alpha    = collect(0.1:0.1:1.14)
  results = zeros(length(alpha));
  energy   = 0.
  energy2  = 0.
  mcs = 10000000    # monte-carlo samplings

@showprogress 1 "COmputing MC simulation"  for a in 1:length(alpha)

    for i in 0:mcs
      x = rand() / sqrt(2) / alpha[a]
      local_energy = (alpha[a])*(alpha[a]) + x*x*(1 - ((alpha[a])^4))
      energy  += local_energy
      energy2 += local_energy*local_energy
    end

    # mean energy
    energy = energy / mcs

    # standard deviation
    sd  = sqrt(((energy2 /mcs ) - (energy/mcs)^2))

    results[a] = energy

  end
  gr(dpi=600)
  plot(xlabel=L"\mathrm{\alpha}", ylabel=L"\mathrm{E_{0}}", framestyle=:box)
  plot!(xtickfontsize=14,ytickfontsize=14,xguidefontsize=14,yguidefontsize=14,legendfontsize=14)
  scatter!(alpha, results, label="MC simulations with N = $mcs")
  plot!(size(400,400))
  savefig("monte_carlo_sp.png")

  println(results[10])

end


text:
module TestDataVectors
    using Base.Test
    using DataArrays

    # Base.getindex(d::DataVector, i::SingleIndex, j::SingleIndex)
    da = @data([1, 2, NA, 4])
    da[1, 1]
    da[2, 1]
    # da[1, 2]

    # Base.push!{T}(dv::DataVector{T}, v::NAtype)
    push!(da, NA)

    # Base.push!{S, T}(dv::DataVector{S}, v::T)
    push!(da, 6)

    # Base.pop!(dv::DataVector)
    pop!(da)

    # Base.unshift!{T}(dv::DataVector{T}, v::NAtype)
    unshift!(da, NA)

    # Base.unshift!{S, T}(dv::DataVector{S}, v::T)
    unshift!(da, -1)

    # Base.shift!{T}(dv::DataVector{T})
    pop!(da)

    # Base.map(f::Function, dv::DataVector)
    map(sin, da)

    # Base.push!{T,R}(pdv::PooledDataVector{T,R}, v::NAtype)
    pda = @pdata([1, 2, NA, 4])
    push!(pda, NA)

    # Base.push!{S,R,T}(pdv::PooledDataVector{S,R}, v::T)
    push!(pda, 6)

    # Base.pop!(pdv::PooledDataVector) = pdv.pool[pop!(pdv.refs)]
    pop!(pda)

    # Base.unshift!{T,R}(pdv::PooledDataVector{T,R}, v::NAtype)
    unshift!(pda, NA)

    # Base.unshift!{S,R,T}(pdv::PooledDataVector{S,R}, v::T)
    unshift!(pda, 6)

    # Base.shift!(pdv::PooledDataVector) = pdv.pool[shift!(pdv.refs)]
    shift!(pda)

    # Base.reverse(x::AbstractDataVector) = x[end:-1:1]
    reverse(da)
    reverse(pda)

    # padNA(dv::AbstractDataVector, front::Integer, back::Integer)
    padNA(da, 5, 5)
    padNA(pda, 5, 5)
end


text:
function device_range(ps_m::CanonicalModel,
                        range_data::Vector{NamedMinMax},
                        time_range::UnitRange{Int64},
                        cons_name::Symbol,
                        var_name::Symbol)

    ps_m.constraints[cons_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, [r[1] for r in range_data], time_range)

    for t in time_range, r in range_data

            ps_m.constraints[cons_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, r[2].min <= ps_m.variables[var_name][r[1], t] <= r[2].max)

    end

    return

end

function device_semicontinuousrange(ps_m::CanonicalModel,
                                    scrange_data::Vector{NamedMinMax},
                                    time_range::UnitRange{Int64}, cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

    ub_name = Symbol(cons_name,:_ub)
    lb_name = Symbol(cons_name,:_lb)

    #MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it. In the future this can be updated
    set_name = [r[1] for r in scrange_data]
    ps_m.constraints[ub_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[lb_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

    for t in time_range, r in scrange_data

            if r[2].min == 0.0

                ps_m.constraints[ub_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] <= r[2].max*ps_m.variables[binvar_name][r[1], t])
                ps_m.constraints[lb_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] >= 0.0)

            else

                ps_m.constraints[ub_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] <= r[2].max*ps_m.variables[binvar_name][r[1], t])
                ps_m.constraints[lb_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] >= r[2].min*ps_m.variables[binvar_name][r[1], t])

            end

    end

    return

end

function reserve_device_semicontinuousrange(ps_m::CanonicalModel,
                                            scrange_data::Vector{NamedMinMax},
                                            time_range::UnitRange{Int64},
                                            cons_name::Symbol,
                                            var_name::Symbol,
                                            binvar_name::Symbol)

    ub_name = Symbol(cons_name,:_ub)
    lb_name = Symbol(cons_name,:_lb)

    #MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it. In the future this can be updated
    set_name = [r[1] for r in scrange_data]
    ps_m.constraints[ub_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)
    ps_m.constraints[lb_name] = JuMP.Containers.DenseAxisArray{JuMP.ConstraintRef}(undef, set_name, time_range)

    for t in time_range, r in scrange_data

            if r[2].min == 0.0

                ps_m.constraints[ub_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] <= r[2].max*(1-ps_m.variables[binvar_name][r[1], t]))
                ps_m.constraints[lb_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] >= 0.0)

            else

                ps_m.constraints[ub_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] <= r[2].max*(1-ps_m.variables[binvar_name][r[1], t]))
                ps_m.constraints[lb_name][r[1], t] = JuMP.@constraint(ps_m.JuMPmodel, ps_m.variables[var_name][r[1], t] >= r[2].min*(1-ps_m.variables[binvar_name][r[1], t]))

            end

    end

    return

 end


text:
module RoME

using Reexport

@reexport using IncrementalInference
@reexport using TransformUtils
@reexport using Distributions
@reexport using KernelDensityEstimate
@reexport using ApproxManifoldProducts

using
  Dates,
  FileIO,
  Distributed,
  LinearAlgebra,
  Statistics,
  Rotations,
  CoordinateTransformations,
  JLD2,
  ProgressMeter,
  DocStringExtensions,
  DistributedFactorGraphs,
  TensorCast,
  ManifoldsBase

using StaticArrays

# to avoid name conflicts
import Manifolds
using Manifolds: hat, ProductGroup, ProductManifold, SpecialEuclidean, ProductRepr, SpecialOrthogonal, TranslationGroup, identity_element

import Manifolds: project, project!

export SpecialOrthogonal, SpecialEuclidean

# using Graphs,  # TODO determine how many parts still require Graphs still directly

import Base: +, \, convert
import TransformUtils: ⊖, ⊕, convert, ominus, veeQuaternion
import IncrementalInference: MB
import IncrementalInference: convert, getSample, reshapeVec2Mat, DFG
# not sure why this is gives import error
import DistributedFactorGraphs: compare
import DistributedFactorGraphs: getDimension, getManifold

# const AMP = ApproxManifoldProducts

const InstanceType{T} = Union{Type{<:T},T}

export
  IIF,  # Aliases for various packages
  KDE,
  TU,
  AMP,
  DFG,

  # RoME specific functions
  measureMeanDist,
  predictBodyBR,
  calcPosePointBearingRange,
  getLastPose,
  getLastPose2D,
  odomKDE,
  initFactorGraph!,
  addOdoFG!,
  addposeFG!,
  newLandm!,
  addBRFG!,
  addMMBRFG!,
  addAutoLandmBR!,
  projNewLandm!,
  malahanobisBR,
  \,
  RangeAzimuthElevation,

  # helper functions
  get2DSamples,
  get2DSampleMeans,
  getAll2DMeans,
  getAll2DPoses,
  get2DPoseSamples,
  get2DPoseMeans,
  get2DPoseMax,
  get2DLandmSamples,
  get2DLandmMeans,
  get2DLandmMax,

  # helper functions
  getLastLandm2D,
  getLastPose2D,
  getNextLbl,

  # RobotUtils
  getRangeKDEMax2D,
  getLastPoses,
  setSolvableOldPoses!,

  # some transform functions
  cart2pol,
  pol2cart,

  # Feature tracking code
  Feature,
  initTrackersFrom,
  propAllTrackers!,
  measUpdateTrackers!,
  assocMeasWFeats!,

  lsrBR,

  # Didson model
  evalPotential,
  LinearRangeBearingElevation,
  project!,
  project,
  # backprojectRandomized!,
  residual!,
  residualLRBE!,
  reuseLBRA,
  ominus,
  ominus!,
  +,
  evalPotential,
  getSample!,
  getSample,

  # Didson convenience function
  addLinearArrayConstraint,

  # camera model -- TODO --separate out
  CameraIntrinsic,
  CameraExtrinsic,
  CameraModelFull,
  project!,
  project,
  # keep
  cameraResidual!,

  # Point2D
  Point2,
  Point2Point2,
  PackedPoint2Point2,
  Point2Point2Range,
  PackedPoint2Point2Range,
  PriorPoint2,
  PackedPriorPoint2,
  Pose2Point2BearingRange,
  Pose2Point2BearingRangeMH,
  PackedPose2Point2BearingRange,
  PackedPose2Point2BearingRangeMH,
  Pose2Point2Bearing,
  PackedPose2Point2Bearing,
  Pose2Point2Range,
  PackedPose2Point2Range,
  PriorPoint2,
  PackedPriorPoint2,

  # Velocity in Point2 types
  DynPoint2,
  DynPoint2VelocityPrior,
  DynPoint2DynPoint2,
  VelPoint2VelPoint2,
  Point2Point2Velocity,
  PackedDynPoint2VelocityPrior,
  PackedVelPoint2VelPoint2,

  Point3Point3,
  PackedPoint3Point3,

  # likely to be deprecated
  solveLandm,
  solvePose2,
  solveSetSeps,
  addPose2Pose2!,
  compareDensity,

  # Pose2D
  Pose2,
  PriorPose2,
  PackedPriorPose2,
  PartialPriorYawPose2,
  PackedPartialPriorYawPose2,
  Pose2Pose2,
  PackedPose2Pose2,
  # Will be deprecated
  addPose2Pose2,

  # Polar types
  Polar,
  PolarPolar,
  PriorPolar,

  # MultipleFeatures2D constraint functions
  MultipleFeatures2D,
  getUvecScaleFeature2D,
  getUvecScaleBaseline2D,

  # Pose3, Three dimensional
  Pose3,
  Point3,
  PriorPoint3,

  # partial Pose3
  PriorPose3ZRP,
  Pose3Pose3XYYaw,
  PackedPose3Pose3XYYaw,
  PackedPriorPose3ZRP,
  # Various utilities
  passTypeThrough,

  # SLAM specific functions
  SLAMWrapper,

  # FG Analysis tools
  rangeErrMaxPoint2,
  rangeCompAllPoses,
  rangeCompAllPoses,

  # new robot navigation functionality
  triggerPose,
  GenericInSituSystem,
  InSituSystem,
  makeInSituSys,
  makeGenericInSituSys,
  advOdoByRules,
  poseTrigAndAdd!,
  poseTrigAndAdd!,
  processTreeTrackersUpdates!,
  addSoftEqualityPoint2D,
  vectoarr2,
  basicFactorGraphExample,
  predictVariableByFactor,

  # jld required Features Type
  LaserFeatures

include("SpecialDefinitions.jl")

#uses DFG v0.10.2 @defVariable for above
include("variables/Local_Manifold_Workaround.jl")
include("variables/VariableTypes.jl")

## More factor types
# RoME internal factors (FYI outside factors are easy, see Caesar documentation)
include("factors/Point2D.jl")
include("factors/Range2D.jl")
include("factors/Bearing2D.jl")
include("factors/BearingRange2D.jl")
include("factors/Polar.jl")
include("factors/PriorPose2.jl")
include("factors/PartialPriorPose2.jl")
include("factors/Pose2D.jl")
include("factors/Pose2Point2.jl")
include("factors/MutablePose2Pose2.jl")
include("factors/DynPoint2D.jl")
include("factors/VelPoint2D.jl")
include("factors/DynPose2D.jl")
include("factors/VelPose2D.jl")
include("factors/Point3D.jl")
include("factors/Point3Point3.jl")
include("factors/Pose3D.jl")
include("factors/Pose3Pose3.jl")
include("factors/PartialPose3.jl")
include("factors/MultipleFeaturesConstraint.jl")
include("factors/InertialPose3.jl")

# tools that come and go
include("TemporaryFunctionality.jl")

# additional tools
include("FactorGraphAnalysisTools.jl")

# tools related to robotics
include("BayesTracker.jl")
include("SensorModels.jl")
include("CameraModel.jl")
include("Slam.jl")
include("RobotUtils.jl")
include("SimulationUtils.jl")
include("OdometryUtils.jl")
include("RobotDataTypes.jl")
include("NavigationSystem.jl")

# generate canonical graphs
include("canonical/GenerateCommon.jl")
include("canonical/GenerateCircular.jl")
include("canonical/GenerateBox.jl")
include("canonical/GenerateHexagonal.jl")
include("canonical/GenerateHoneycomb.jl")
include("canonical/GenerateBeehive.jl")
include("canonical/GenerateHelix.jl")

# more utils requiring earlier functions
include("AdditionalUtils.jl")
include("g2oParser.jl")

# ScalarFields
include("services/ScalarFields.jl")

# things on their way out
include("Deprecated.jl")

# optional tools
using Requires

function __init__()
  # combining neural networks natively into the non-Gaussian  factor graph object
  @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
    @info "Loading RoME.jl tools related to Flux.jl."
    include("factors/flux/models/Pose2OdoNN_01.jl") # until a better way is found to deserialize
    include("factors/flux/MixtureFluxPose2Pose2.jl")
  end

  # Scalar field specifics

  @require ImageCore = "a09fc81d-aa75-5fe9-8630-4744c3626534" begin
    @require ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" include("services/RequiresImages.jl")
  end
  # Images="916415d5-f1e6-5110-898d-aaa5f9f070e0"

  # @require Interpolations="a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" begin
  #   include("services/ScalarFieldsInterpolations.jl")
  # end
end

# manifold conversions required during transformation

end


text:
module SVG

export Style, Svg, Polyline, Circle, Line, write, translate, bounds

import Base.==

dp(n; digits=2) = round(n; digits)

abstract type SvgObject end

struct Style
	fill::String
	strokecolor::String
	strokewidth::Float64
	Style(;fill="none", strokecolor="black", strokewidth=1) = new(fill, strokecolor, strokewidth)
end

struct Polyline <: SvgObject
	xs
	ys
	style::Style
	Polyline(; style::Style=Style()) = new(Vector{Float64}(), Vector{Float64}(), style)
	Polyline(npoints::Int; style::Style=Style()) = new(Vector{Float64}(undef, npoints), Vector{Float64}(undef, npoints), style)
	Polyline(xs, ys; fx=identity, fy=identity, style::Style=Style()) = new(map(fx, xs), map(fy, ys), style)
	Polyline(xys; fx=identity, fy=identity, style::Style=Style()) = new(map(xy->fx(xy[1]), xys), map(xy->fy(xy[2]), xys), style)
end

==(p1::Polyline, p2::Polyline) = p1.xs == p2.xs && p1.ys == p2.ys && p1.style == p2.style

struct Line <: SvgObject
	x1
	y1
	x2
	y2
	style::Style
	Line(x1, y1, x2, y2; style::Style=Style()) = new(x1, y1, x2, y2, style)
end

struct Circle <: SvgObject
	x
	y
	r
	s::Style
	Circle(x, y, r; style::Style=Style()) = new(x, y, r, style)
end

struct Svg
	objects::Vector{SvgObject}
	Svg() = new(Vector{SvgObject}())
	Svg(objs) = new(objs)
end

function bounds(svg::Svg)
	xmin, xmax, ymin, ymax = Inf, -Inf, Inf, -Inf
	for b in map(bounds, svg.objects)
		if b.xmin < xmin
			xmin = b.xmin
		end
		if b.xmax > xmax
			xmax = b.xmax
		end
		if b.ymin < ymin
			ymin = b.ymin
		end
		if b.ymax > ymax
			ymax = b.ymax
		end
	end
	xmin, ymin = max(0,xmin), max(0,ymin)
	(;xmin, xmax, ymin, ymax)
end

Base.write(fn::String, svg::Svg, width, height; viewbox="", inhtml=false, digits=2, objwrite_fn=write_objs) = open(fn, "w+") do io write(io, svg, width, height; viewbox, inhtml, digits, objwrite_fn) end

write_objs(io::IO, svg::Svg) = foreach(o->write(io, o), svg.objects)

function Base.write(io::IO, svg::Svg, width, height; viewbox="", inhtml=false, digits=2, objwrite_fn=write_objs)
	if inhtml
		println(io, "
")
	end
	if viewbox == ""
		viewbox = ""
	else
		viewbox = " viewBox=\"$viewbox\""
	end
	println(io, """<svg width="$width" height="$height"$viewbox>""")
	objwrite_fn(io, svg)
	println(io, "</svg>")
	if inhtml
		println(io, "
")
	end
end

Base.write(io::IO, s::Style) = print(io, "style=\"", "fill:", s.fill, ";stroke:", s.strokecolor, ";stroke-width:", dp(s.strokewidth), "\"")

function Base.write(io::IO, p::Polyline)
	print(io, "<polyline points=\"")
	foreach(i->print(io, p.xs[i], ",", p.ys[i], " "), 1:length(p.xs))
	print(io, "\" ")
	write(io, p.style)
	println(io, " />")
end

bounds(p::Polyline) = (xmin=minimum(p.xs), xmax=maximum(p.xs), ymin=minimum(p.ys), ymax=maximum(p.ys))

function Base.write(io::IO, L::Line)
	print(io, "<line x1=\"", L.x1, "\" y1=\"", L.y1, "\" x2=\"", L.x2, "\" y2=\"", L.y2, "\" ")
	write(io, L.style)
	println(io, " />")
end

bounds(l::Line) = (xmin=min(l.x1,l.x2), xmax=max(l.x1,l.x2), ymin=min(l.y1,l.y2), ymax=max(l.y1,l.y2))

function Base.write(io::IO, c::Circle)
	print(io, "<circle cx=\"", c.x, "\" cy=\"", c.y, "\" r=\"", c.r, "\" ")
	write(io, c.s)
	println(io, " />")
end

points2circles(xs, ys; r=1) = map(i->Circle(xs[i], ys[i], r), 1:length(xs))

bounds(c::Circle) = (xmin=c.x-c.r, xmax=c.x+c.r, ymin=c.y-c.r, ymax=c.y+c.r)

#helper
translate(points::Vector, delta) = map(p->p+delta, points)

###
end

max_stars_repo_name: lawless-m/jsvg

max_issues_repo_name: lawless-m/jsvg

max_forks_repo_name: lawless-m/jsvg


text:
## Exercise 5-4
## If you are given three sticks, you may or may not be able to arrange them in a triangle. For example, if one of the sticks is 12 inches long and the other two are one inch long, you will not be able to get the short sticks to meet in the middle. For any three lengths, there is a simple test to see if it is possible to form a triangle:

## TIP
## If any of the three lengths is greater than the sum of the other two, then you cannot form a triangle. Otherwise, you can. (If the sum of two lengths equals the third, they form what is called a “degenerate” triangle.)

## 1. Write a function named istriangle that takes three integers as arguments, and that prints either “Yes” or “No”, depending on whether you can or cannot form a triangle from sticks with the given lengths.
println("Ans 1: ")

function istriangle(a, b, c)
   return (((a + b) > c) && ((b + c) > a) && ((c + a) > b)) ? "Yes" : "No"
end

println(istriangle(13, 1, 1))
println(istriangle(4, 3, 5))

## 2. Write a function that prompts the user to input three stick lengths, converts them to integers, and uses istriangle to check whether sticks with the given lengths can form a triangle.
println("Ans 2: ")

function check_is_triangle()
    print("Enter the first side: ")
    a = parse(Int64, readline())
    print("Enter the second side: ")
    b = parse(Int64, readline())
    print("Enter the third side: ")
    c = parse(Int64, readline())

    println(istriangle(a, b, c))
end

check_is_triangle()

println("End.")


text:
using MPI  # needed to generate docs for parallel HDF5 API
using HDF5
using Documenter

makedocs(;
    modules=[HDF5],
    authors="Mustafa Mohamad <mus-m@outlook.com> and contributors",
    repo="https://github.com/JuliaIO/HDF5.jl/blob/{commit}{path}#L{line}",
    sitename="HDF5.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaIO.github.io/HDF5.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Low-level library bindings" => "api_bindings.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaIO/HDF5.jl",
)


text:
function generate_integrand_uv(kernel, testref, trialref, testel, trialel)

    function k3(u,v)
        out = @SMatrix zeros(3,3)

        x = neighborhood(testel,u)
        y = neighborhood(trialel,v)

        kernelval = kernel(x,y)
        f = testref(x)
        g = trialref(y)

        jx = jacobian(x)
        jy = jacobian(y)
        ds = jx*jy

        return  SMatrix{3,3}([dot(f[i][1], kernelval*g[j][1])*ds for i=1:3, j=1:3])
    end

    return k3
end


text:
#=
   6 digit numbers

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   Page 19f, Problem 2.5 Problems concerning random numbers

   Given the 6 digits numbers:

   a) Problem 1
      find the probability that at least one of the digits 0..9
      appears exactly twice.
      Answer: 2943/4000 ~ 0.7358

   b) Problem 2
      find the probability that at least two of the digits 0..9
      appears exactly once.
      Answer: 1179/1250 ~ 0.9432

=#

using Turing # , StatsPlots, DataFrames
include("jl_utils.jl")

# Note that we collect the interesting data to post process it.
# This is because we cannot add non random values in the chain.
k_occ_coll1 = []
# occ_coll = []
@model function six_digits(n, k, m)
        # d = TArray{Int}(undef, n)
        d = Vector{Int}(undef, n)
        for i in 1:n
            d[i] ~ DiscreteUniform(0,9)
        end

        # Number of occurrences of each digit in digits
        occ = make_hash(d)

        # push!(occ_coll, occ)
        # Number of digits that occurs exactly k times
        k_occ = sum(values(occ).|>x->x==k)>=m

        push!(k_occ_coll1, k_occ)
end

# 1) Find the probability that at least one of the digits 0..9
#    appears exactly twice.
println("1) Find the probability that at least one of the digits 0..9\nappears exactly twice.")
model = six_digits(6,2,1)
num_chains = 4
num_samples = 10_000
chains = sample(model, MH(), MCMCThreads(), 1000, num_chains)

# It seems that MH() is the only sampler that can handle this
# chains = sample(model, MH(), 10000) #
# chains = sample(model, IS(), MCMCThreads(), num_samples, num_chains)
# chains = sample(model, PG(20), MCMCThreads(), 1000, num_chains)
# chains = sample(model, SMC(1000), MCMCThreads(), 1000, num_chains)

# NUTS, HMC, and HMCDA throws this error:
# chains = sample(model, NUTS(), 1000) # Error
# chains = sample(model, HMC(0.1, 5), 1000) # Same error as NUTS
# chains = sample(model, HMCDA(0.15, 0.65), 1000) # Same errors as NUTS

# display(chains)

# display(group(chains,:d))

# println("occ_coll:$occ_coll")
println("mean: $(sum(k_occ_coll1)/length(k_occ_coll1))")
display(sort(make_hash(k_occ_coll1)))

# Note: Since we are updating an external array we have to
# build a new model here...

# 2) Find the probability that at least two of the digits 0..9
#    appears exactly once.
k_occ_coll2 = []
@model six_digits_2(n, k, m) = begin
        # d2 = TArray{Int}(undef, n)
        d2 = Vector{Int}(undef, n)
        for i in 1:n
            d2[i] ~ DiscreteUniform(0,9)
        end

        # Number of occurrences of each digit in digits
        occ2 = make_hash(d2)

        # Number of digits that occurs exactly k times
        k_occ2 = sum(values(occ2).|>x->x==k)>=m
        push!(k_occ_coll2, k_occ2)
end

println("\n2) Find the probability that at least two of the digits 0..9\nappears exactly once.")
model2 = six_digits_2(6,1,2)

num_chains = 4
num_samples = 10_000
chains = sample(model2, MH(), MCMCThreads(), num_samples, num_chains)

println("mean: $(sum(k_occ_coll2)/length(k_occ_coll2))")
display(sort(make_hash(k_occ_coll2)))

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@testset "dist RF" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    tree2 = ParseNewick("((E:5.0,(B:5.0,A:5.0)C:9.0)F:5.0,D:5.0)G:5.0;")
    @test RF(tree, tree2) == 2
end

@testset "dist RF_weighted" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    tree2 = ParseNewick("((E:5.0,(B:5.0,A:5.0)C:9.0)F:5.0,D:5.0)G:5.0;")
    @test RF_weighted(tree, tree2) == 2/12
end


text:
using Compat
import Compat.String

mutable struct RawEEG{S<:Real, P<:String}
    data::AbstractMatrix{S}
    eventTable::Any
    sampRate::Real
    nChannels::Int
    nSamples::Int
    chanLabels::AbstractVector{P}
end

function toRawEEG(data::AbstractMatrix{S}, eventTable, sampRate::Real, ; chanLabels=["" for i=1:size(data)[1]]) where {S<:Real}
    thisRawEEG = RawEEG(data, eventTable, sampRate, size(data)[1], size(data)[2], chanLabels)

    return thisRawEEG
end


text:
module RConnector

using Sockets
using Logging
import REPL
import Main.Tables: istable, columnaccess,
                    columns, AbstractColumns,
                    getcolumn, columnnames

const TYPE_ID_NOTHING = 0x00
const TYPE_ID_FLOAT64 = 0x01
const TYPE_ID_COMPLEX = 0x02
const TYPE_ID_RAW = 0x03
const TYPE_ID_INT = 0x04
const TYPE_ID_BOOL = 0x05
const TYPE_ID_STRING = 0x06
const TYPE_ID_LIST = 0x07
const TYPE_ID_EXPRESSION = 0xee
const TYPE_ID_CALLBACK = 0xcb
const TYPE_ID_NAMED_FUNCTION = 0xfc
const TYPE_ID_OBJECT_REFERENCE = 0xce
const TYPE_ID_SYMBOL = 0x5b

const OBJECT_CLASS_ID_ARRAY = 0xaa
const OBJECT_CLASS_ID_SIMPLE_ARRAY = 0x5a
const OBJECT_CLASS_ID_ANONYMOUS_FUNCTION = 0xaf
const OBJECT_CLASS_ID_STRUCT = 0x5c
const OBJECT_CLASS_ID_NO_INFO = 0x00

const CALL_INDICATOR = 0x01
const RESULT_INDICATOR = 0x00
const STDOUT_INDICATOR = 0x50
const STDERR_INDICATOR = 0x5e
const FAIL_INDICATOR = 0xff
const BYEBYE = 0xbb

const NO_ATTRIBUTES = Base.ImmutableDict{String, Any}()

const SEND_AS_RAW_TYPES = Union{UInt64, Int128, UInt128, Ptr}
const SEND_AS_INT32 = Union{Int8, Int16, UInt16, Char}
const SEND_AS_DOUBLE = Union{Float16, Float32, UInt32}
const SEND_AS_COMPLEX = Complex{T} where T <: Union{
      Int8, Int16, Int32, Int64, Float16, Float32}

const R_NA_INTEGER = Int32(-2147483648)
const R_NA_REAL = ntoh(reinterpret(Float64,
      UInt8[0x7f, 0xf0, 0x00, 0x00, 0x00, 0x00, 0x07, 0xa2])[1])
const R_NA_COMPLEX = Complex{Float64}(R_NA_REAL, R_NA_REAL)

include("communicating.jl")
include("sharing.jl")
include("handling_undefined.jl")
include("reading.jl")
include("handling_dataframes.jl")
include("evaluating.jl")
include("writing.jl")
include("exporting.jl")
include("accessing_mutating.jl")

end


text:
using MatrixCompletion

import Distributions

@testset "$(format("SGD Optimizer: Bernoulli Loss [Small][Forgiving][Native]"))" begin
  let
    input_size = 500
    y = rand(Distributions.Bernoulli(0.6), input_size) * 1.0
    mle_x = MatrixCompletion.Losses.sgd_train(Loss{Bernoulli}(),
                                              fx         = rand(input_size),
                                              y          = y,
                                              c          = zeros(input_size),
                                              ρ          = 0,
                                              α          = 0.2,
                                              ρ₁         = 0.9,
                                              ρ₂         = 0.999,
                                              batch_size = 500,
                                              epoch      = 10);
    prediction = predict(Bernoulli(),forward_map(Bernoulli(),mle_x))
    tc = provide(Diagnostics{Poisson()}(),
                 input_data=prediction, reference=y)
    display(tc)
  end
end