text:
n, p = 100, 5
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p; seed=52)

@testset "linreg" begin
    lr  = LinearRegression(fit_intercept=false)
    lr1 = LinearRegression()
    β_ref = X \ y
    @test β_ref == fit(lr, X, y)

    # fit_intercept
    β_ref = X1 \ y1
    @test β_ref == fit(lr1, X, y1)

    # == iterative solvers
    β_cg = fit(lr1, X, y1; solver=CG())
    @test norm(β_cg - β_ref) / norm(β_ref) ≤ 1e-12
end

@testset "ridgereg" begin
    λ   = 1.0
    rr  = RidgeRegression(lambda=λ, fit_intercept=false)
    rr1 = RidgeRegression(λ, penalize_intercept=true)

    β_ref  = (X'X + λ*I) \ (X'y)
    β_ref1 = (X1'X1 + λ*I) \ (X1'y1)
    @test β_ref ≈ fit(rr, X, y)
    @test β_ref1 ≈ fit(rr1, X, y1)

    β_cg  = fit(rr, X, y; solver=CG())
    β_cg1 = fit(rr1, X, y1; solver=CG())

    @test norm(β_cg - β_ref) / norm(β_ref)   ≤ 1e-12
    @test norm(β_cg1 - β_ref1) / norm(β_ref) ≤ 1e-12
end

n, p = 500, 100
((X, y, θ), (X1, y1, θ1)) = generate_continuous(n, p;  seed=512, sparse=0.1)

@testset "lasso" begin
    # no intercept
    λ  = 50
    lr = LassoRegression(λ; fit_intercept=false)
    J  = objective(lr, X, y)
    θ_ls    = X \ y
    θ_ista  = fit(lr, X, y, solver=ISTA())
    θ_fista = fit(lr, X, y)
    @test isapprox(J(θ_ls),    680.26525, rtol=1e-5)
    @test isapprox(J(θ_fista), 631.58622, rtol=1e-5) # <- ref value
    @test isapprox(J(θ_ista),  631.58622, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 100  # not sparse
    @test nnz(θ_fista) == 12 # sparse
    @test nnz(θ_ista)  == 12

    # with intercept
    lr1 = LassoRegression(λ, penalize_intercept=true)
    J   = objective(lr1, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(lr1, X, y1)
    θ_ista  = fit(lr1, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    353.85709, rtol=1e-5)
    @test isapprox(J(θ_fista), 312.93487, rtol=1e-5) # <- ref values
    @test isapprox(J(θ_ista),  312.93487, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 9   # sparse
    @test nnz(θ_ista)  == 9

    # with intercept and not penalizing intercept
    lr1 = LassoRegression(λ)
    J   = objective(lr1, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(lr1, X, y1)
    θ_ista  = fit(lr1, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    353.71225, rtol=1e-5)
    @test isapprox(J(θ_fista), 312.93073, rtol=1e-5) # <- ref values
    @test isapprox(J(θ_ista),  312.93073, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 10  # sparse
    @test nnz(θ_ista)  == 10

    if DO_COMPARISONS
        lr_sk = SKLEARN_LM.Lasso(alpha=λ/n)
        lr_sk.fit(X, y1)
        θ1_sk = vcat(lr_sk.coef_[:], lr_sk.intercept_)
        @test isapprox(J(θ1_sk), 312.93072, rtol=1e-5)
        @test nnz(θ1_sk) == 10
    end
end

@testset "elnet" begin
    ρ = 0.3
    α = 0.1
    λ = α * (1 - ρ) * n
    γ = α * ρ * n
    enr = ElasticNetRegression(λ, γ; penalize_intercept=true)
    J   = objective(enr, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(enr, X, y1)
    θ_ista  = fit(enr, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    220.38693, rtol=1e-5)
    @test isapprox(J(θ_fista), 199.79860, rtol=1e-5)
    @test isapprox(J(θ_ista),  199.79860, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 10   # sparse
    @test nnz(θ_ista)  == 10

    # not penalizing intercept
    enr = ElasticNetRegression(λ, γ)
    J   = objective(enr, X, y1)
    θ_ls    = X1 \ y1
    θ_fista = fit(enr, X, y1)
    θ_ista  = fit(enr, X, y1, solver=ISTA())
    @test isapprox(J(θ_ls),    220.34333, rtol=1e-5)
    @test isapprox(J(θ_fista), 199.78497, rtol=1e-5)
    @test isapprox(J(θ_ista),  199.78497, rtol=1e-5)
    # sparsity
    @test nnz(θ_ls)    == 101 # not sparse
    @test nnz(θ_fista) == 11   # sparse
    @test nnz(θ_ista)  == 11

    if DO_COMPARISONS
        enr_sk = SKLEARN_LM.ElasticNet(alpha=α, l1_ratio=ρ)
        enr_sk.fit(X, y1)
        θ_sk = vcat(enr_sk.coef_[:], enr_sk.intercept_)
        @test isapprox(J(θ_sk), 199.78496, rtol=1e-5)
        @test nnz(θ_sk) == 11
    end
end


text:
struct Widget
    type
    body
end

function display(w::Widget)
    HTML(w.body)
end

function widget_graph(
        v2he::String,
        he2v::String,
        div_id::Int;
        v_meta::Union{AbstractVector, Nothing}=nothing,
        he_meta::Union{AbstractVector, Nothing}=nothing,
        width::Int=500,
        height::Int=500,
        radius::Real=10,
        node_radii::Union{AbstractVector{<:Real}, Nothing}=nothing,
        node_color::String="#999",
        node_colors::Union{AbstractVector{String}, Nothing}=nothing,
        node_stroke::Union{String, Nothing} = nothing,
        node_strokes::Union{AbstractVector{String}, Nothing}=nothing,
        stroke_width::Real=0,
        stroke_widths::Union{AbstractVector{<:Real}, Nothing}=nothing,
        node_opacity::Real=1,
        node_opacities::Union{AbstractVector{<:Real}, Nothing}=nothing,
        stroke_opacity::Real=1,
        stroke_opacities::Union{AbstractVector{<:Real}, Nothing}=nothing,
        with_node_labels::Bool=false,
        node_labels::Union{AbstractVector{String}, Nothing}=nothing,
        with_node_metadata_hover::Bool=false,
        with_node_weight::Bool=false,
        he_colors::Union{AbstractVector{String}, Nothing}=nothing,
        with_he_labels::Bool=false,
        he_labels::Union{AbstractVector{String}, Nothing}=nothing,
        with_he_metadata_hover::Bool=false
        )

    w = Widget("graph", """

            """)
    return w
end

ext: jl


text:
module AxisIndicesGraphs

using AxisIndices, LightGraphs, SimpleWeightedGraphs
using Base: to_index, @propagate_inbounds

struct AxisIndicesGraph{T,G<:AbstractGraph{T},AI} <: AbstractGraph{T}
    graph::G
    axes::AI
end

Base.parent(g::AxisIndicesGraph) = getfield(g, :graph)

# Although it seems weird to overload this for graphs it makes sense in the context
# of how an "axis" is defined (e.g., keys + indices)
Base.axes(g::AxisIndicesGraph) = getfield(g, :axes)

# Shouldn't these be performed on the type
#eltype(g::AbstractMetaGraph) = eltype(g.graph)
#edgetype(g::AbstractMetaGraph) = edgetype(g.graph)

###
### Vertices
###
LightGraphs.nv(g::AxisIndicesGraph) = nv(parent(g))
LightGraphs.vertices(g::AxisIndicesGraph) = vertices(parent(g))

LightGraphs.has_vertex(g::AxisIndicesGraph, x...) = has_vertex(parent(g), x...)
@inline LightGraphs.has_edge(g::AxisIndicesGraph, x...) = has_edge(parent(g), x...)

@propagate_inbounds function LightGraphs.add_vertex!(g::AxisIndicesGraph, u, v,)
    return add_vertex!(parent(g), to_index(axes(g), u), to_index(axes(g), v))
end

@propagate_inbounds function LightGraphs.rem_vertex!(g::AxisIndicesGraph, u, v)
    return rem_vertex!(parent(g), to_index(axes(g), u), to_index(axes(g), v))
end

###
### Edges
###
LightGraphs.ne(g::AxisIndicesGraph) = ne(parent(g))
LightGraphs.edges(g::AxisIndicesGraph) = edges(parent(g))

@propagate_inbounds function LightGraphs.add_edge!(g::AxisIndicesGraph, u, v)
    return add_edge!(parent(g), to_index(axes(g), u), to_index(axes(g), v))
end
@propagate_inbounds function LightGraphs.add_edge!(g::AxisIndicesGraph, u, v, val)
    return add_edge!(parent(g), to_index(axes(g), u), to_index(axes(g), v), val)
end

@propagate_inbounds function LightGraphs.rem_edge!(g::AxisIndicesGraph, u, v)
    return rem_edge!(parent(g), to_index(axes(g), u), to_index(axes(g), v))
end

function LightGraphs.inneighbors(g::AxisIndicesGraph, v)
    return inneighbors(parent(g), to_index(axes(g), v))
end
function LightGraphs.outneighbors(g::AxisIndicesGraph, v)
    return outneighbors(parent(g), to_index(axes(g), v))
end

Base.Base.issubset(g::AxisIndicesGraph, h::AbstractGraph) = issubset(parent(g), h)

LightGraphs.is_directed(::Type{<:AxisIndicesGraph{T,G}}) where {T,G} = is_directed(G)

SimpleWeightedGraphs.weighttype(::Type{<:AxisIndicesGraph{T,G}}) where {T,G} = weighttype(G)

# TODO these should be AxisIndicesArrays
SimpleWeightedGraphs.degree_matrix(g::AxisIndicesGraph, args...) = degree_matrix(parent(g), args...)
function LightGraphs.adjacency_matrix(g::AxisIndicesGraph, args...)
    return adjacency_matrix(parent(g), args...)
end

function LightGraphs.laplacian_matrix(g::AxisIndicesGraph, args...)
    return laplacian_matrix(parent(g), args...)
end

function LightGraphs.pagerank(g::AxisIndicesGraph, α=0.85, n::Integer=100, ϵ=1.0e-6)
    return pagerank(parent(g), α, n, ϵ)
end

function SimpleWeightedGraphs.get_weight(g::AxisIndicesGraph, u, v)
    return get_weight(parent(g), to_index(axes(g), u), to_index(axes(g), v))
end

LightGraphs.connected_components(g::AxisIndicesGraph) = connected_components(parent(g))

function LightGraphs.induced_subgraph(g::AxisIndicesGraph, vlist)
    return induced_subgraph(parent(g), to_index(axes(g), vlist))
end

###
### Base
###
Base.zero(g::AxisIndicesGraph) = AxisIndicesGraph(zero(parent(g)), axes(g))

Base.copy(g::AxisIndicesGraph) = AxisIndicesGraph(copy(parent(g)), copy(axes(g)))

end


text:
function m_i(multi_index::S) where {S <: Integer}

  if multi_index == 1
	  m_node_number = one(S)
  else
	  m_node_number = 2^(multi_index-1)+1
   end

  return m_node_number

end

function m_i(multi_index::Array{S,1}) where {S <: Integer}

  m_node_number = Array{S,1}(undef,length(multi_index))

  @inbounds for i in eachindex(multi_index)
    m_node_number[i] = m_i(multi_index[i])
  end

  return m_node_number

end

function m_i(multi_index::Array{S,2}) where {S <: Integer}

  m_node_number = Array{S,2}(undef,size(multi_index))

  @inbounds for i in eachindex(multi_index)
    m_node_number[i] = m_i(multi_index[i])
  end

  return m_node_number

end


text:
# For very large inputs, it is practical to preprocess data as much as possible.
# So in distributed memory mode, we expect inputs to:
#   - have normalized all user and item ids to be integers starting from 1
#   - have user-item and item-user (the transpose) rating data pre-created so that we don't have to transpose a large matrix every time
#   - have filtered out all empty users and items from the data
type DistInputs <: Inputs
    Rfile::FileSpec
    RTfile::FileSpec
    nusers::Int
    nitems::Int
    R::Nullable{InputRatings}
    RT::Nullable{InputRatings}

    function DistInputs(ratings_file::FileSpec, transposed_ratings_file::FileSpec)
        new(ratings_file, transposed_ratings_file, 0, 0, nothing, nothing)
    end
end

function clear(inp::DistInputs)
    inp.R = nothing
    inp.RT = nothing
end

localize!(inp::DistInputs) = nothing
share!(inp::DistInputs) = nothing

function ensure_loaded(inp::DistInputs)
    if isnull(inp.R)
        # R is the user x item matrix (users in rows and items in columns)
        R = read_input(inp.Rfile)
        RT = read_input(inp.RTfile)
        inp.nusers, inp.nitems = size(R)

        inp.R = R
        inp.RT = RT
    end
    nothing
end

item_idmap(inp::DistInputs) = Int64[]
user_idmap(inp::DistInputs) = Int64[]


text:
using Pkg
Pkg.update()
Pkg.develop(path = ".")
Pkg.add("DashBase")
Pkg.add("PlotlyBase")
Pkg.add("HTTP")
Pkg.instantiate()
Pkg.build("Dash")
Pkg.build("HTTP")
Pkg.test("Dash", coverage=true)


text:
#### query.
# # switch to exact search.
# function searchmoleculelist(list_names::Vector{String},
#                 target_strings::Vector{String})
#
#     M = length(target_strings)
#     N = length(list_names)
#
#     inds = collect(1:N)
#
#     target_exist_flags = falses(M)
#
#
#     keep_flags = falses(N)
#     for i = 1:N
#         db_string = Unicode.normalize(list_names[i], casefold=true)
#
#         for m = 1:M
#
#             search_string = Unicode.normalize(target_strings[m], casefold=true)
#
#             if db_string == search_string
#                 keep_flags[i] = true
#                 target_exist_flags[m] = true
#             end
#         end
#     end
#
#     return inds[keep_flags], target_exist_flags
# end

# single entry.
function searchmoleculelist(list_names::Vector{String},
                target_string::String)

    N = length(list_names)

    inds = collect(1:N)

    target_exist_flag = false

    search_string = Unicode.normalize(target_string, casefold=true)

    keep_flags = falses(N)
    for i = 1:N

        db_string = Unicode.normalize(list_names[i], casefold=true)

        if search_string == db_string
            keep_flags[i] = true
            target_exist_flag = true
        end

    end

    return inds[keep_flags], target_exist_flag
end

########## molar mass look up.

function getmolarmasssolvent(name)::Float64
    #

    if name == "D2O"
        return 20.028 #g/mol, D2O.
    end

    return NaN
end

# in grams/mol.
function getmolarmass(name)::Float64
    #

    if name == "DSS"
        return 218.32177
    end

    if name == "TSP"
        return 163.94
    end

    if name == "TMSP"
        return 146.26
    end

    if name == "TMS"
        return 88.22
    end

    return NaN
end

function loadGISSMOmolecule(base_path, file_name)

   load_path = joinpath(base_path, "$(file_name).jld")

   dic = JLD.load(load_path)

   css_sys = dic["css_sys"]
   J_IDs_sys = dic["J_IDs_sys"]
   J_vals_sys = dic["J_vals_sys"]
   cs_singlets = dic["cs_singlets"]

   J_lb_sys = dic["J_lb_sys"]
   J_ub_sys = dic["J_ub_sys"]
   cs_lb_sys = dic["cs_lb_sys"]
   cs_ub_sys = dic["cs_ub_sys"]
   cs_LUT = dic["cs_LUT"]

   ref_mM = dic["ref_mM"]
   solute_mM = dic["solute_mM"]
   solubility_gL_solute = dic["solubility_gL_solute"]
   solubility_gL_ref = dic["solubility_gL_ref"]
   molar_mass_ref = dic["molar_mass_ref"]
   molar_mass_solute = dic["molar_mass_solute"]
   pH = dic["pH"]
   temperature = dic["temperature"]
   ref_molecule_name = dic["ref_molecule_name"]

   return css_sys, J_IDs_sys, J_vals_sys, cs_singlets,
            J_lb_sys, J_ub_sys, cs_lb_sys, cs_ub_sys, cs_LUT,
            ref_mM, solute_mM, solubility_gL_solute,
            solubility_gL_ref, molar_mass_ref,
            molar_mass_solute, pH, temperature, ref_molecule_name
end


text:
using Isaac
using Isaac.TestSets
using Isaac.TestSets: precond

println("Testing nsolimod for index-1 saddles")
@testset "nsolidmod-index1" begin

#[1] Muller Test
xe = []
for init in (:near, :far)
   println("Testing nsolimod, Muller, init = $(init)")
   V = MullerPotential()
   x0, v0 = ic_dimer(V, init)
   E, dE = objective(V)
   x, ndE = nsolimod(dE, x0, 1; V0=v0, maxnumdE = 100, maxstep = 0.5)
   println("   num_dE = ", ndE)
   if init == :near
      xe = copy(x)
   else
      println("  |x - xe| = ", norm(x - xe, Inf))
      @test norm(x - xe, Inf) < 1e-8
   end
end

# [2] Vacancy Test
V = LJVacancy2D(bc = :clamped)
E, dE = objective(V)
x0, v0 = ic_dimer(V, :near)
P = precond(V, x0)
Iprep = (P, x) -> P
Pprep = (P, x) -> Isaac.TestSets.precond(V, x)
xe = []

for (init, precond, precon_prep) in
      [(:near, I, Iprep), (:near, P, Pprep), (:far, I, Iprep), (:far, P, Pprep)]

   x0, v0 = ic_dimer(V, init)
   pre = !(precond == I)

   println("Testing nsolimod, 2D LJ Vacancy, init = $init, precond = $pre")
   x, ndE = nsolimod(dE, x0, 1; V0=v0, P = precond, precon_prep=precon_prep)
   println("   num_dE = ", ndE)
   @test norm(dE(x), Inf) < 1e-5
   println("   |dE| = ", norm(dE(x), Inf))
   if xe == []
      # to make sure all 4 methods converge to the same saddle,
      # store the first one ...
      xe = x
      H2 = TestSets.hessian(V, x)
      σ = sort(eigvals(Symmetric(H2)))
      @test σ[1] * σ[2] < 0
      println("   min-eigs = ", σ[1:2])
   else
      # ... and then compare the others
      @test norm(x - xe, Inf) < 1e-6
      println("   |x - xe| = ", norm(x - xe, Inf))
   end
end

end  # @testset "nsolidstab-index1"


text:
function initializeACTR(;parms...)
    chunks = Chunk[Chunk(;isa=:bafoon,animal=:dog,name=:Sigma),
        Chunk(;isa=:mammal,animal=:cat,name=:Butters),
        Chunk(;isa=:mammal,animal=:rat,name=:Joy),
        Chunk(;isa=:mammal,animal=:rat,name=:Bernice)]
    memory = Declarative(;memory=chunks)
    imaginal = Imaginal(buffer=typeof(chunks)(undef,1))
    actr = ACTR(;declarative=memory, imaginal=imaginal, parms...)
    return actr,memory,chunks
end


text:
@testset "Avg cut pruning" begin
    algo = AvgCutPruningAlgo(2)
    pruner = CutPruner{2, Int}(algo, :≤)
    # check redundancy between new and old cuts
    pruner.excheck = true

    addposition!(pruner, [1, 0])
    @test 1:1 == addcuts!(pruner, [1 0], [1], [true])
    @test 2:2 == addcuts!(pruner, [0 1], [1], [false])
    addusage!(pruner, [1, 0])
    @test [2, 0, 0] == addcuts!(pruner, [1 1; -1 -1; 0 1], [1, 1, 2], [true, false, true])
    @test pruner.A == [1 0; 1 1]
    @test pruner.b == [1, 1]
    @test pruner.ids == [1, 3]
    @test [0, 0] == addcuts!(pruner, [1 0; 1 1], [1, 1], [true, true])
    @test pruner.A == [1 0; 1 1]
    @test pruner.b == [1, 1]
    @test pruner.ids == [1, 3]
    addusage!(pruner, [1, 0])
    @test [2, 0] == addcuts!(pruner, [2 0; 0 2], [0, 0], [false, false])
end


text:
module LazyFieldArraysTests

using Test
using Gridap.Arrays
using Gridap.Fields
using Gridap.TensorValues
using FillArrays
using Gridap.Helpers

# using BenchmarkTools
using Test

np = 4
p1 = Point(1.0,2.0)
p2 = Point(2.0,1.0)
p3 = Point(4.0,3.0)
p4 = Point(6.0,1.0)

p = p1

np = 4
p1 = Point(1.0,2.0)
p2 = Point(2.0,1.0)
p3 = Point(4.0,3.0)
p4 = Point(6.0,1.0)

p = p1
np = 4
x = [p1,p2,p3,p4]

f1 = GenericField(x->1*x)
f2 = GenericField(x->2*x)
f3 = GenericField(x->3*x)
f4 = GenericField(x->4*x)
f5 = GenericField(x->5*x)

h(x) = 2.0*x
f = GenericField(h)

q(x) = Point(sqrt.(x.data))
g = GenericField(q)

nf = 2
# basis = [f1, f2, f3, f4, f5]
basis = fill(f,nf)

field = g

na = 3

x_a = fill(x,na)
basis_a = fill(basis,na)
field_a = fill(field,na)

v = rand(nf)
v_a = fill(v,na)

vm = rand(nf,nf)
vm_a = fill(vm,na)

# Evaluate field

res_field_a = lazy_map(evaluate,field_a,x_a)

c_r = array_cache(res_field_a)
# @btime getindex!($c_r,$res_field_a,1);

# Evaluate basis

res_basis_a = lazy_map(evaluate,basis_a,x_a)

c_r = array_cache(res_basis_a)
# @btime getindex!($c_r,$res_basis_a,1);

cb = return_cache(basis,x)
# @btime evaluate!(cb,basis,x)

# Transpose basis

tbasis_a = lazy_map(transpose,basis_a)
@test all([ (tbasis_a[i] == transpose(basis_a[i])) for i in 1:na])

res_tbasis_a = lazy_map(evaluate,tbasis_a,x_a)
for j in 1:na
  @test all([res_tbasis_a[j][i,1,:] == res_basis_a[j][i,:] for i in 1:np])
end

c_r = array_cache(res_tbasis_a)
# @btime getindex!($c_r,$res_tbasis_a,1);

# Broadcast operation basis basis

op = +
brbasis_a = lazy_map(Broadcasting(Operation(op)),basis_a,basis_a)
res_brbasis_a = lazy_map(evaluate,brbasis_a,x_a)

c_r = array_cache(res_brbasis_a)
# @btime getindex!($c_r,$res_brbasis_a,$1);

# Broadcast operation basis field

op = +
brbasis_a = lazy_map(Broadcasting(Operation(op)),basis_a,field_a)

res_brbasisfield_a = lazy_map(evaluate,brbasis_a,x_a)
for j in 1:na
  for i in 1:np
    @test res_brbasisfield_a[j][i,:] == res_basis_a[j][i,:] .+ res_field_a[j][i]
  end
end

c_r = array_cache(res_brbasisfield_a)
# @btime getindex!($c_r,$res_brbasisfield_a,1);

# Linear Combination basis values (vector values)

vxbasis_a = lazy_map(linear_combination,v_a,basis_a)
res_vxbasis_a = lazy_map(evaluate,vxbasis_a,x_a)

for j in 1:na
  for i in 1:np
    res_vxbasis_a[j][i] == transpose(res_basis_a[j][i,:])*v_a[j]
  end
end

c_r = array_cache(res_vxbasis_a)
# @btime getindex!(c_r,res_vxbasis_a,1)

# Linear Combination basis values (matrix values)

vmxbasis_a = lazy_map(linear_combination,vm_a,basis_a)
res_vmxbasis_a = lazy_map(evaluate,vmxbasis_a,x_a)

for j in 1:na
  for i in 1:np
    res_vmxbasis_a[j][i] == transpose(res_basis_a[j][i,:])*vm_a[j]
  end
end

c_r = array_cache(res_vmxbasis_a)
# @btime getindex!(c_r,res_vmxbasis_a,1)

# basis*transpose(basis)

basisxtbasis_a = lazy_map(Broadcasting(Operation(⋅)),basis_a,tbasis_a)
res_basisxtbasis_a = lazy_map(evaluate,basisxtbasis_a,x_a)
size(res_basisxtbasis_a[1])
for j in 1:na
  for i in 1:np
    @test res_basisxtbasis_a[j][i,:,:] == broadcast(⋅,res_basis_a[j][i,:],res_tbasis_a[j][i,:,:])
  end
end

c_r = array_cache(res_basisxtbasis_a)
# @btime getindex!($c_r,$res_basisxtbasis_a,1);

# composition basis(field)

op = ∘
compfield_a = lazy_map(op,field_a,field_a)

res_field_a = lazy_map(evaluate,field_a,x_a)

res_compfield_a = lazy_map(evaluate,compfield_a,x_a)

@test all([
  res_compfield_a[i]  ==   lazy_map(evaluate,field_a,res_field_a)[i]
  for i in 1:na])

c_r = array_cache(res_compfield_a)
# @btime getindex!($c_r,$res_compfield_a,1);

# composition basis(basis, field)

op = ∘
compbasisfield_a = lazy_map(Broadcasting(op),basis_a,field_a)

res_compbasisfield_a = lazy_map(evaluate,compbasisfield_a,x_a)
@test all([ res_compbasisfield_a[i] == lazy_map(evaluate,basis_a,res_field_a)[i] for i in 1:na])

c_r = array_cache(res_compbasisfield_a)
# @btime getindex!($c_r,$res_compbasisfield_a,1);

## transpose(basis)*basis
#
#tbasisxbasis_a = lazy_map(*,tbasis_a,basis_a)
#res_tbasisxbasis_a = lazy_map(evaluate,tbasisxbasis_a,x_a)
#for j in 1:na
#  for i in 1:np
#    res_tbasisxbasis_a[j][i] == res_tbasis_a[j][i,1,:] ⋅ res_basis_a[j][i,:]
#  end
#end
#
#c_r = array_cache(res_tbasisxbasis_a)
## @btime getindex!($c_r,$res_tbasisxbasis_a,1);

# integration (field)

w = rand(np)
w_a = fill(w,na)
jac = ConstantField(TensorValue(4.0,0.0,0.0,4.0))
jac_a = fill(jac,na)
res_jac_a = lazy_map(evaluate,jac_a,x_a)

field_a
integrate_a = lazy_map(integrate,field_a,x_a,w_a,jac_a)

# integration (basis)

# w_a = v_a
jac = ConstantField(TensorValue(4.0,0.0,0.0,4.0))
jac_a = fill(jac,na)
res_jac_a = lazy_map(evaluate,jac_a,x_a)

integrate_a = lazy_map(integrate,basis_a,x_a,w_a,jac_a)

# Gradient

∇field_a = lazy_map(gradient,field_a)
res_∇field_a = lazy_map(evaluate,∇field_a,x_a)
@test res_∇field_a[1][1] == TensorValue(0.5, 0.0, 0.0, 0.35355339059327373)

∇basis_a = lazy_map(Broadcasting(∇),basis_a)
res_∇basis_a = lazy_map(evaluate,∇basis_a,x_a)
@test all(res_∇basis_a[1] .== TensorValue(2.0, 0.0, 0.0, 2.0))

end # module


text:
using .Networks

using Dates
using JuMP
using Printf
using Parameters

@with_kw mutable struct EODAssessmentConfigs
    cut_prod_coeff = 100 #must be greater than LoL coeff in objective
    inj_prod_coeff = 10 #must be greater than LoL coeff in objective
    loss_of_load_coeff = 1 #LoL
    out_path = nothing
    problem_name = "EODAssessment"
end

@with_kw struct EODAssessment <: AbstractRunnable
    configs::EODAssessmentConfigs = EODAssessmentConfigs()
end

function run(runnable::EODAssessment, ech::Dates.DateTime, TS::Vector{Dates.DateTime}, context::AbstractContext)
    problem_name_l = @sprintf("eod_assessment_%s", ech)
    runnable.configs.out_path = context.out_dir
    runnable.configs.problem_name = problem_name_l

    # @assert( runnable.configs.loss_of_load_penalty >=
    #         biggest_min_imposition(get_tso_actions(context)) * runnable.configs.cut_prod_penalty)
    @assert( runnable.configs.loss_of_load_coeff
                <= runnable.configs.inj_prod_coeff
                    <= runnable.configs.cut_prod_coeff )

    #FIXME : only create vars for started units

    return eod_assessment(get_network(context),
                        TS,
                        get_assessment_uncertainties(context),
                        get_tso_actions(context),
                        runnable.configs
                    )
end

function eod_assessment(network,
                        TS,
                        assessment_uncertainties,
                        tso_actions,
                        configs::EODAssessmentConfigs
                    )
    model_container_l = formulate_eod_assessment(network, TS, assessment_uncertainties, tso_actions, configs)
    solve!(model_container_l, configs.problem_name, configs.out_path)

    for var in all_variables(model_container_l.model)
        println(name(var), " = ", value(var))
    end

    return model_container_l
end

function biggest_min_imposition(tso_actions::TSOActions)
    return maximum(x->x[2][1], get_impositions(tso_actions))
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
    LaneChangeChoice
A choice of whether to change lanes, and what direction to do it in
"""
const DIR_RIGHT = -1
const DIR_MIDDLE =  0
const DIR_LEFT =  1
struct LaneChangeChoice
    dir::Int # -1, 0, 1
end
Base.show(io::IO, a::LaneChangeChoice) = @printf(io, "LaneChangeChoice(%d)", a.dir)
Base.length(::Type{LaneChangeChoice}) = 1
Base.convert(::Type{LaneChangeChoice}, v::Vector{Float64}) = LaneChangeChoice(convert(Int, v[1]))
function Base.copyto!(v::Vector{Float64}, a::LaneChangeChoice)
    v[1] = a.dir
    v
end
function get_lane_offset(a::LaneChangeChoice, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    if a.dir == DIR_MIDDLE
        rec[pastframe][vehicle_index].state.posF.t
    elseif a.dir == DIR_LEFT
        convert(Float64, get(LANEOFFSETLEFT, rec, roadway, vehicle_index, pastframe))
    else
        @assert(a.dir == DIR_RIGHT)
        convert(Float64, get(LANEOFFSETRIGHT, rec, roadway, vehicle_index, pastframe))
    end
end

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

abstract type LaneChangeModel end
get_name(::LaneChangeModel) = "???"
set_desired_speed!(::LaneChangeModel, v_des::Float64) = model # # do nothing by default
reset_hidden_state!(model::LaneChangeModel) = model # do nothing by default
observe!(model::LaneChangeModel, scene::Scene, roadway::Roadway, egoid::Int) = model  # do nothing by default
Base.rand(model::LaneChangeModel) = error("rand not implemented for model $model")

max_stars_repo_name: intelligent-control-lab/Autoenv

max_issues_repo_name: intelligent-control-lab/Autoenv

max_forks_repo_name: intelligent-control-lab/Autoenv


text:
#----------------------------------------------------------------------------------------------#
#                            cloze.jl -- Moodle cloze question type                            #
#----------------------------------------------------------------------------------------------#


text:
module EthylDynamic
  # Function to call the GUI
  export EthylDynamicGUI

  using Gtk, Gtk.ShortNames, JLD, Suppressor, CSV, Mustache, Dates
  import DataFrames

  # Path to CSS Gtk-Style dataFile
  global style_file = joinpath(dirname(Base.source_path()), "style2020.css")

  # Logo
  if Sys.iswindows()
    global image_path = joinpath(dirname(Base.source_path()), "media\\mainlogo.png")
  end

  if Sys.islinux()
    global image_path = joinpath(dirname(Base.source_path()), "media//mainlogo.png")
  end

  # General Settings
  global pathDataBase = joinpath(dirname(Base.source_path()), "database.csv")

  # Load default database
  global database = CSV.read(pathDataBase)

  # Include the main file .fl
  include("mainGUI_EthylDynamic.jl")
end


text:
mutable struct Graph
    nodes::Dict
    relationships::Dict
end

Graph() = Graph(Dict(), Dict())

abstract type Entity end

mutable struct Node <: Entity
    graph::Graph
    id::Integer
    labels::Set
    properties::Dict
end

function Node(graph::Graph, id::Integer)
    if !haskey(graph.nodes, id)
        graph.nodes[id] = Node(graph, id, Set(), Dict())
    end
    return graph.nodes[id]
end

function Base.:(==)(x::Node, y::Node)
    return x.id == y.id && x.labels == y.labels && x.properties == y.properties
end

function Base.hash(x::Node)
    return Base.hash(x.properties, Base.hash(x.id))
end

function Base.keys(node::Node)
    Base.keys(node.properties)
end

function Base.values(node::Node)
    Base.values(node.properties)
end

mutable struct Relationship <: Entity
    graph::Graph
    id::Integer
    type::String
    properties::Dict
    start_node::Union{Node, Nothing}
    end_node::Union{Node, Nothing}
end

function Relationship(graph::Graph, id::Integer, type::String)
    if !haskey(graph.relationships, id)
        graph.relationships[id] = Relationship(graph, id, type, Dict(), nothing, nothing)
    end
    return graph.relationships[id]
end

function Base.hash(x::Relationship)
    return Base.hash(x.properties, Base.hash(x.id))
end

function Base.keys(r::Relationship)
    Base.keys(r.properties)
end

function Base.values(r::Relationship)
    Base.values(r.properties)
end

mutable struct Path
    nodes::Vector
    relationships::Vector
end

function Path(start_node::Node, relationships::Vector{Relationship})
    nodes = [start_node]
    for (i, relationship) in enumerate(relationships)
        if relationship.start_node == nodes[length(nodes)]
            push!(nodes, relationship.end_node)
        elseif relationship.end_node == nodes[length(nodes)]
            push!(nodes, relationship.start_node)
        else
            throw(ErrorException("Relationship $(i) does not connect to the last node"))
        end
    end
    return Path(nodes, relationships)
end

function Base.length(p::Path)
    return length(p.relationships)
end

function Base.:(==)(x::Path, y::Path)
    return x.nodes[1] == y.nodes[1] && x.relationships == y.relationships
end

function Base.hash(x::Path)
    return Base.hash(x.nodes[1], Base.hash(x.relationships))
end

function start_node(path::Path)
    return path.nodes[1]
end

function end_node(path::Path)
    return path.nodes[length(path.nodes)]
end

function put_node(graph::Graph, n_id::Integer, labels::Set, properties::Dict=Dict(); kwproperties...)
    for (k, v) in kwproperties
        if v == nothing
            continue
        end
        properties[String(k)] = v
    end

    node = Node(graph, n_id)
    union!(node.labels, labels)
    merge!(node.properties, properties)
    return node
end

put_node(graph::Graph, n_id::Integer, labels::Vector, properties::Dict=Dict(); kwproperties...) = put_node(graph, n_id, Set(labels), properties; kwproperties...)

function put_relationship(graph::Graph, r_id::Integer, start_node::Node, end_node::Node, r_type::String, properties::Dict=Dict(); kwproperties...)
    r = put_unbound_relationship(graph, r_id, r_type, properties::Dict; kwproperties...)
    r.start_node = start_node
    r.end_node = end_node
    return r
end

function put_unbound_relationship(graph::Graph, r_id::Integer, r_type::String, properties::Dict=Dict(); kwproperties...)
    for (k, v) in kwproperties
        if v == nothing
            continue
        end
        properties[String(k)] = v
    end

    relationship = Relationship(graph, r_id, r_type)
    merge!(relationship.properties, properties)
    return relationship
end

function hydrate_path(nodes::Vector, relationships::Vector, sequence::Vector)
    if length(nodes) == 0
        throw(ErrorException("Empty nodes array"))
    elseif length(sequence) % 2 != 0
        throw(ErrorException("Invalid sequence length"))
    end

    last_node = nodes[1]
    entities = Relationship[]
    for (i, rel_index) in enumerate(sequence[1:2:length(sequence)])
        if rel_index == 0
            throw(ErrorException("Invalid rel_index"))
        end
        n = sequence[2*i]
        if n < 0
            n = length(nodes) + n
        end
        next_node = nodes[n+1] #index
        if rel_index > 0
            r = relationships[rel_index] #index
            r.start_node = last_node
            r.end_node = next_node
            push!(entities, r)
        else
            r = relationships[-rel_index] #index
            r.start_node = next_node
            r.end_node = last_node
            push!(entities, r)
        end
        last_node = next_node
    end
    return Path(nodes[1], entities)
end

function graph_hydration_functions(graph::Graph)
    return Dict(
        UInt8('N')=>(psh, values)->put_node(graph, values...),
        UInt8('R')=>(psh, values)->put_relationship(graph, values[1], Node(graph, values[2]), Node(graph, values[3]), values[4], values[5]),
        UInt8('r')=>(psh, values)->put_unbound_relationship(graph, values...),
        UInt8('P')=>(psh, values)->hydrate_path(values...)
    )
end


text:
""" Conjugate gradient algorithm using a line search supporting stopping
    Generic version receiving a conjugacy formula
    kwargs may specify options for the line search.
"""
function CG_generic(nlp       :: AbstractNLPModel;
                    x         :: Vector{T}=copy(nlp.meta.x0),
                    stp       :: NLPStopping = NLPStopping(nlp,
                                                           NLPAtX(nlp.meta.x0)),
                    scaling   :: Bool = true,
                    LS_algo   :: Function = bracket{T},
                    LS_logger :: AbstractLogger = Logging.NullLogger(),
                    CG_formula :: Function = formula_HZ,
                    kwargs...      # eventually options for the line search
                    ) where T

    @info log_header([:iter, :f, :dual, :step, :slope], [Int, T, T, T, T],
                     hdr_override=Dict(:f=>"f(x)", :dual=>"‖∇f‖", :slope=>"IterLS"))

    n = nlp.meta.nvar

    f = obj(nlp,x)
    ∇f = grad(nlp, x)

    xt = similar(x)
    ∇ft = similar(∇f)

    τ₀ = 0.05
    τ₁ = 0.16

#    ϕ, ϕstp = prepare_LS(stp, x, ∇f, τ₀, f, ∇f)

    OK = update_and_start!(stp, x = x, fx = f, gx = ∇f)
    #update_and_stop!(stp,  x = x, fx = f, gx = ∇f)
    ϕ, ϕstp = prepare_LS(stp, x, ∇f, τ₀, f, ∇f)
    @info log_row(Any[0, f, norm(∇f)])

    β = 0.0
    d = zeros(n)
    scale = 1.0

    while !OK
        d = - ∇f + β*d

        slope = ∇f⋅d
        if slope > 0.0  # restart with negative gradient
            d = - ∇f
          slope =  ∇f⋅d
        end

        # Simple line search call
        # returns  t, xt = x + t*d,  ft=f(xt), ∇ft = ∇f(xt)

        t, xt, ft, ∇ft = linesearch(ϕ, ϕstp, x, d*scale, f, ∇f, τ₀, τ₁, logger = LS_logger, algo = LS_algo; kwargs...)
        fail_sub_pb = ~ϕstp.meta.optimal

        if fail_sub_pb
            OK = true
            stp.meta.fail_sub_pb = true
        else
            #move on
            s = xt - x
            y = ∇ft - ∇f
            β = 0.0
            #if (∇ft⋅∇f) < 0.75 * (∇ft⋅∇ft)   # Powell restart
            β = CG_formula(∇f,∇ft,s,d)
            #end
            if scaling
                scale = (y⋅s) / (y⋅y)
            end
            if scale <= 0.0
                scale = 1.0
            end
            x .= xt
            f = ft
            ∇f .= ∇ft
            OK = update_and_stop!(stp, x = x, gx = ∇f, fx = f)
        end

        @info log_row(Any[stp.meta.nb_of_stop, f, norm(∇f), t, ϕ.counters.neval_obj])
    end

    if !stp.meta.optimal
        @warn status(stp,list=true)
        @debug x
    end

    return stp
end

max_stars_repo_name: vepiteski/LSDescent

max_issues_repo_name: vepiteski/LSDescent

max_forks_repo_name: vepiteski/LSDescent


text:
using Pkg, Random
using Dates
using PyPlot, Logging
import Base.show

struct Node
    id::Int16
    x::Float32
    y::Float32
end

function show(io::IO, n::Node)
    print(string("Node:", n.id))
end

import Base.copy

mutable struct Firefly
    path::Array{Node}
    cost::Float32
end

function show(io::IO, n::Firefly)
    print(string("Firefly:", n.cost, " "))
end

copy(s::Firefly) = Firefly(s.path, s.cost)

function arc_distance(f1::Firefly, f2::Firefly)::Int
    """Calculates the arc distance between two firefly paths.
    Example:
        Firefly i -> 1 2 4 6 3 5 7
        Firefly j -> 1 2 4 3 5 7 6
        1-2, 2-4, 3-5, 5-7 arcs in firefly i are also presented
        in firefly j. But 4-6, 6-3 arcs are missing in
        firefly j. so the arc distance between firefly i
        and firefly j is 2 here.
    """
    missmatched_pair_counts = 0
    f1_pairs = zip(f1.path[1:end-1], f1.path[2:end])
    f2_pairs = zip(f2.path[1:end-1], f2.path[2:end])

    for f1_pair in f1_pairs
        if f1_pair ∉ f2_pairs missmatched_pair_counts += 1 end
    end

    return missmatched_pair_counts
end

function attractiveness(source_node::Firefly, destination_node::Firefly, λ::Float32)::Float32
    """Firefly’s attractiveness is proportional to the light
    intensity seen by adjacent fireflies, we can now define the
    attractiveness β of a firefly depending on distance between
    two fireflies.
    """
    β = source_node.cost
    r, _ = hamming_distance(source_node, destination_node)
    return β * exp(-λ ^ (r^2))
end

function create_distance_matrix(nodes::Array{Node})
    """Creates and returns the matrix of distances between given nodes.
    Uses euclidean_distance function to calculate distance.
    """
    node_num = length(nodes)
    result_matrix = zeros(node_num, node_num)

    for n1 in nodes
        for n2 in nodes
            result_matrix[n1.id,n2.id] = haversine_distance(n1, n2)
        end
    end
    return result_matrix
end

function generate_mutated_fireflies(firefly::Firefly, distance_matrix, r::Int, k::Int)::Array{Firefly}
    """Generates mutated fireflies using inversion_mutation function.
    Returns the new generated firefly array.
    """
    return [inversion_mutation(copy(firefly), distance_matrix, r) for _ in 1:k]
end

function hamming_distance(firefly_1::Firefly, firefly_2::Firefly)
    """Calculates the hamming distance between two fireflies.
    """
    dist = 0
    dist_info = []
    for (i, n) in enumerate(firefly_1.path)
        if n.id != firefly_2.path[i].id
            dist += 1
            push!(dist_info, true)
        else
            push!(dist_info, false)
        end
    end
    return dist, dist_info
end

function haversine_distance(n1::Node, n2::Node)::Float32
    #converting decimal degree to radians
    lat1=n1.y
    long1=n1.x
    lat2=n2.y
    long2=n2.x
    lat1,long1,lat2,long2=map(deg2rad,[lat1,long1,lat2,long2])

    # Using Haversine formula for distance calculation
    dlat=lat2-lat1
    dlong=long2-long1
    a=sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlong/2)^2
    c=2 * asin(sqrt(a))
    r=6371 #radius of earth in km
    return c*r
end

function init_firefly_paths(fireflies::Array{Firefly}, distance_matrix)
    """Initializes random paths for the given firefly array.
    """
    println("inside init_firefly_paths")
    for f in fireflies
        shuffle!(f.path)
        f.cost = path_cost(f, distance_matrix)
    end
end

function inversion_mutation(firefly::Firefly, distance_matrix, r::Int)::Firefly
    """Reverses the random part of given firefly path with
    given difference.
    """
    length_of_mutation = rand(2:r);
    max_len = length(firefly.path)-length_of_mutation
    index1 = rand(1:max_len)
    index2 = index1 + length_of_mutation
    mutated_path = reverse(firefly.path, index1, index2)
    new_f = copy(firefly)
    new_f.path = mutated_path
    new_f.cost = path_cost(new_f, distance_matrix)
    return new_f
end

function light_intensity(source_node::Firefly, destination_node::Firefly, λ::Float32)::Float32
    """Calculates the light intensity for the source firefly.
    """
    I_0 = source_node.cost
    r, _ = hamming_distance(source_node, destination_node)
    return I_0 * ℯ ^ (-λ * (r^2))
end

function move_firefly(f1::Firefly, f2::Firefly, r::Int)
    """Moving f1 to f2 less than r swap operations.
    """
    number_of_swaps = Int(ceil(rand() * (r-2) ))
    d, d_info = hamming_distance(f1, f2)

    while number_of_swaps > 0
        d, d_info = hamming_distance(f1, f2)
        random_index = rand([i for i in 1:length(d_info) if d_info[i]])
        value_to_copy = f2.path[random_index]
        index_to_move = findfirst(x -> x==value_to_copy, f1.path)

        if number_of_swaps == 1 && f1.path[index_to_move] == f2.path[random_index] && f1.path[random_index] == f2.path[index_to_move]
            break
        end

        f1.path[random_index], f1.path[index_to_move] = f1.path[index_to_move], f1.path[random_index]
        # println("swapped! ", f1.path[random_index], " -- ", f1.path[index_to_move])
        if f1.path[index_to_move] == f2.path[index_to_move] number_of_swaps -= 1 end
        number_of_swaps -= 1
    end
end

function path_cost(firefly::Firefly, distance_matrix)::Float32
    """Calculates the total path cost from given distance matrix.
    """
    path = firefly.path
    path_len = length(path) - 1
    p_cost = distance_matrix[path[end].id, path[1].id]
    for i in range(1, length=path_len)
        p_cost +=  distance_matrix[path[i].id, path[i+1].id]
    end

    return p_cost
end

function read_nodes(cluster_info)::Array{Node}
    """Reading nodes from giving file path.
    """
    nodes = Node[]
    for coord in cluster_info
            id = parse(Int16, coord[1])
            x_coord = parse(Float32, coord[2])
            y_coord = parse(Float32, coord[3])
            push!(nodes, Node(id, x_coord, y_coord))
    end
    return nodes
end

function plot_firefly(best_x, best_y, x, y, x_cost, best_cost, file_name, tsp_file_name)

    plot(best_x, best_y, "-")
    plot(x, y, "ro", markersize=2.0)
    xlabel("X Dimension")
    ylabel("Y Dimension")
    title("Firefly Cost: $best_cost", fontsize=10)
    suptitle("File: $file_name", fontweight="bold")
    fig_name = "graphs/best_firefly_path_$file_name-cost_$best_cost.png"
    savefig(fig_name)
    #println("Saved figure:",fig_name)
    clf()
    plot(x_cost)
    grid(true)
    ylabel("Path Cost")
    xlabel("Iteration Number")
    title("Firefly Iteration Graphic", fontsize=12)
    fig_name = "graphs/iteration_$tsp_file_name-cost_$best_cost.png"
    savefig(fig_name)
    clf()
end

function main(cluster_data, name)
    final_route=[]

    for ncluster in 1:length(cluster_data)

        file_name=name[ncluster]
        nodes = read_nodes(cluster_data[ncluster])
        distance_matrix = create_distance_matrix(nodes)
        # Checking if logs direcitory is available!
        # If is not create one!
        if !isdir("logs") mkdir("logs"); end

        tsp_file_name = file_name
        io_log = open("logs/$tsp_file_name.txt", "w+")

        # Setting Logger!
        logger = SimpleLogger(io_log)
        global_logger(logger)

        # Setting hyperparameters!
        LIGHT_ABSORPTION_COEFF = 0.2
        ATTRACTION_COEFF = 1
        ITERATION_NUMBER = 100
        POPULATION_NUMBER = 150
        NUMBER_OF_MUTATION = 10

        # Main firefly algorithm
        bests = []
        pop = [Firefly(copy(nodes), -1.0) for _ in 1:POPULATION_NUMBER]
        init_firefly_paths(pop, distance_matrix)
        println(length(pop), " Sized population created!")
        for t in 1:ITERATION_NUMBER
            new_pop = []
            for f1 in pop
                for f2 in pop
                    if f1.cost < f2.cost
                        mutated_f1s = generate_mutated_fireflies(f2, distance_matrix, arc_distance(f1, f2), NUMBER_OF_MUTATION)
                        new_pop = [new_pop ; mutated_f1s]
                        push!(new_pop, f2)
                    end
                end
            end
            [push!(pop, n) for n in new_pop]
            sorted_pop = sort!(pop, by=p->p.cost)
            while length(pop) != POPULATION_NUMBER
                pop!(pop)
            end
            push!(bests, sorted_pop[1])
            current_best_cost = copy(bests[end].cost)
            now = string(Dates.now())
            println("Step:",t, " - Best ", bests[end])
            @info("Step:$t, Population Size:$POPULATION_NUMBER, Mutation Number:$NUMBER_OF_MUTATION, Date:$now",
                  current_best_cost)
        end

        close(io_log)

        # Checking if graphs direcitory is available for saving graphs!
        # If is not create one!
        if !isdir("graphs") mkdir("graphs"); end

        # Plotting and printing the best solutions.
        best = sort(pop, by=p->p.cost)[1]
        best_x = [n.x for n in best.path]
        best_y = [n.y for n in best.path]
        best_id= [n.id for n in best.path]
        best_cost = best.cost
        x = [n.x for n in nodes]
        y = [n.y for n in nodes]
        x_cost = [b.cost for b in bests]
        plot_firefly(best_x, best_y, x, y, x_cost, best_cost, file_name, tsp_file_name)
        z=[best_id, best_x, best_y]
        push!(final_route, z)
    end
    return final_route
end


text:
#!/usr/bin/env julia

using Test

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions
using .SFB.SeparableArrays

## This might be quicker instead of loading the full SFB module:
#include("../src/SeparableArrays.jl")
#using .SeparableArrays

@testset "SeparableArrays" begin
    @testset "SeparableArrays matrix functionality" begin
        phi = rand(5)
        mask = rand(3)
        A = phi * mask'  # our reference
        S = SeparableArray(phi, mask; name1=:phi, name2=:mask)
        for i in eachindex(S)
            @test S[i] == A[i]
        end
        for i=1:length(S)
            @test S[i] == A[i]
        end
        for i=1:size(S,1), j=1:size(S,2)
            @test S[i,j] == A[i,j]
        end
        @test all(S .== A)
        @test S[5,3,1] == A[5,3,1]
        @test_throws BoundsError S[5,3,2]
        @test_throws BoundsError S[6,1]
        @test_throws BoundsError S[2,4]
        @test_throws BoundsError S[0,1]
        @test_throws BoundsError S[1,0]
        @test_throws BoundsError S[6,0]
        @test ndims(S) == ndims(A)
        @test S.arr1 == phi
        @test S.arr2 == mask
        @test S.phi == phi
        @test S.mask == mask
        @test propertynames(S) == (:phi, :mask)
        S2 = @SeparableArray phi mask
        @test S2 == S
        @test propertynames(S2) == (:phi, :mask)
        @test S2.phi == S.phi
        @test S2.mask == S.mask

        # squaring
        S3 = S .^ 2
        S4 = S .* S
        #S5 = [S[i,j]*S[i,j] for i=1:size(S,1),j=1:size(S,2)]
        S5 = exponentiate(S, 2)
        @show S3 S4 S5
        @show (S4 .- S3) ./ eps.(S3)
        @show (S5 .- S3) ./ eps.(S3)
        @test all(abs.(S4 .- S3) .<= 10*eps.(S4))
        @test all(abs.(S5 .- S3) .<= 10*eps.(S5))
        @test all(abs.(S5 .- S4) .<= 10*eps.(S5))
        #@test typeof(S3) <: SeparableArray
        #@test typeof(S4) <: SeparableArray
        @test typeof(S5) <: SeparableArray
    end

    @testset "SeparableArrays ndims>2 functionality" begin
        phi = rand(2,3)
        mask = rand(4,5)
        A = [phi[i,j]*mask[k,l]
             for i=1:size(phi,1), j=1:size(phi,2), k=1:size(mask,1), l=1:size(mask,2)]
        S = SeparableArray(phi, mask; name1=:phi, name2=:mask)

        pass = true
        for i in eachindex(S)
            S[i] == A[i] || (pass = false)
            #@show i S[i] A[i]
            #@assert pass
        end
        @test pass

        pass = true
        for i=1:length(S)
            S[i] == A[i] || (pass = false)
        end
        @test pass

        pass = true
        for i=1:size(S,1), j=1:size(S,2), k=1:size(S,3), l=1:size(S,4)
            S[i,j,k,l] == A[i,j,k,l] || (pass = false)
        end
        @test pass

        @test all(S .== A)

        @test_throws BoundsError S[5,3,2,3]
        @test_throws BoundsError S[1,0,1,1]
        @test_throws BoundsError S[6,0,1,1]
        @test_throws BoundsError S[0,1,1,1]
        @test_throws BoundsError S[0,6,1,1]
        @test_throws BoundsError S[1,1,1,0]
        @test_throws BoundsError S[1,1,6,0]
        @test_throws BoundsError S[1,1,0,1]
        @test_throws BoundsError S[1,1,0,6]
        @test ndims(S) == ndims(A)
        @test S.arr1 == phi
        @test S.arr2 == mask
        @test S.phi == phi
        @test S.mask == mask
        @test propertynames(S) == (:phi, :mask)
        S2 = @SeparableArray phi mask
        @test S2 == S
        @test propertynames(S2) == (:phi, :mask)
        @test S2.phi == S.phi
        @test S2.mask == S.mask
    end

    @testset "SeparableArrays (1,1) timing and allocations" begin
        phi = rand(1000)
        mask = rand(3400)
        # our reference
        A = phi * mask'
        S = @SeparableArray phi mask

        S[1,1]  # compile
        @time S[1,1]

        S[:,110]
        S[110,:]
        @time S[:,110]
        @time S[110,:]
        @time phi[110]*mask
        @time phi[110]*mask
    end

    @testset "SeparableArrays (2,2) timing and allocations" begin
        phi = rand(100,120)
        mask = rand(340,110)
        # our reference
        A = [phi[i] * mask[j,k] for i=1:length(phi), j=1:size(mask,1), k=1:size(mask,2)]
        S = @SeparableArray phi mask

        S[1,1,1,1]  # compile
        @time S[1,1,1,1]

        S[1,:,2,1]
        S[1,2,2,:]
        @time S[1,:,2,1]
        @time S[1,2,2,:]
    end
end

# vim: set sw=4 et sts=4 :


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# This code was ported from a MATLAB implementation available at
# http://www.albertostrumia.it/Fractals/FractalMatlab/Julia/juliaSH.m.
# Original MATLAB implementation is copyright (c) the author.

using DocOpt
if "--demo" in ARGS
    using Winston
end

include("ndgrid.jl") # provides the meshgrid function

function juliaset(demo::Bool)
    col = 128             # color depth
    m = 1000              # image size
    cx = 0                # center X
    cy = 0                # center Y
    l = 1.5               # span

    if demo
        iters = 50
        zoomAmount = 0.9
    else
        iters = 10
        zoomAmount = 0.6
    end

    # The c constant.
    c = -.745429 + .11308*im

    for zoom = 1:iters

        # `x` and `y` are two 1000-element arrays representing the x
        # and y axes: [-1.5, -1.497, ..., 0, ..., 1.497, 1.5] on the
        # first iteration of this loop.
        x = linspace(cx-l, cx+l, m)
        y = linspace(cy-l, cy+l, m)

        # `X` and `Y` are two arrays containing, respectively, the x-
        # and y-coordinates of each point on a 1000x1000 grid.
        (X, Y) = meshgrid(x, y)

        # Let `Z` represent the complex plane: a 1000x1000 array of
        # numbers each with a real and a complex part.
        Z = X + im*Y

        # Iterate the Julia set computation (squaring each element of
        # Z and adding c) for `col` steps.
        W = iterate(col, Z, c)

        # Mask out the NaN values (overflow).
        minval = minimum(W)
        W[isnan(W)] = minval - minval/10

        if demo
            # Plot the result as a "heat map"-style plot.  The values
            # of the elements of W determine the "color" of the
            # corresponding point.
            plot = Winston.imagesc(W)
            Winston.display(plot)
        end

        # Zoom into the next frame, shrinking the distance that `x`
        # and `y` will cover.
        l = l * zoomAmount
    end

end

function iterate(col, Z, c)
    for k = 1:col
        Z = Z.*Z .+ c
    end
    return exp(-abs(Z))
end

function main()
    doc = """juliaset.jl

Julia set fractal visualization.

Usage:
  juliaset.jl -h | --help
  juliaset.jl [--demo]

Options:
  -h --help  Show this screen.
  --demo     Use settings that look good for the purposes of a demo, and plot the results using Winston.
"""
    arguments = docopt(doc)
    demo = arguments["--demo"]

    tic()
    iterate(0, Array{Complex{Float64}}(0,2), Complex(0.0))
    time = toq()
    println("SELFPRIMED ", time)

    tic()
    juliaset(demo)
    time = toq()
    println("SELFTIMED ", time)
end

main()


text:
## spherical.jl : implementation of spherical covariance function

## Spherical ##
struct Spherical{T} <: IsotropicCovarianceStructure{T}
    λ::T
    σ::T
    p::T

    function Spherical{T}(λ::T, σ::T, p::T) where T
        λ > 0 || throw(DomainError(λ, "correlation length λ of spherical covariance cannot be negative or zero"))
        σ > 0 || throw(DomainError(σ, "marginal standard deviation σ of spherical covariance cannot be negative or zero"))
        p >= 1 || throw(DomainError(p, "in p-norm, p must be greater than or equal to 1"))
        isinf(p) && throw(DomainError(p, "in p-norm, p cannot be infinity"))

        new{T}(λ, σ, p)
    end
end

"""
    Spherical(λ, [σ = 1], [p = 2])

Spherical covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm, defined as

``C(x, y) = \\begin{cases} σ \\left(1 - \\displaystyle\\frac{3}{2}\\frac{ρ}{λ} + \\frac{1}{2}\\left(\\frac{ρ}{λ}\\right)^3\\right) & \\text{for }ρ≤λ\\\\0 & \\text{for }ρ>λ\\end{cases}``

with ``ρ = ||x - y||_p``.

# Examples
```jldoctest
julia> Spherical(0.1)
spherical (λ=0.1, σ=1.0, p=2.0)

julia> Spherical(1.0, σ=2)
spherical (λ=1.0, σ=2.0, p=2.0)

```
See also: [`Exponential`](@ref), [`Linear`](@ref), [`Whittle`](@ref), [`Gaussian`](@ref), [`SquaredExponential`](@ref), [`Matern`](@ref)
"""
Spherical(λ::Real; σ::Real=1.0, p::Real=2) = Spherical{promote_type(typeof(λ),typeof(σ),typeof(p))}(promote(λ, σ, p)...)

# evaluate spherical covariance
apply(s::Spherical, x::Real) = max(zero(x), 1 - 3/2 * x/s.λ + 1/2 * (x/s.λ)^3)

# short name
shortname(::Spherical) = "spherical"


text:
@testset "archmCopula.jl" begin
    @test length(claytonCopula(2,3,2)) > 0
    @test length(claytonCopula(2,4,2)) > 0
end


text:
using Test
using DFTK

include("testcases.jl")

if mpi_nprocs() == 1  # not easy to distribute
@testset "High-symmetry kpath construction for silicon" begin
    testcase = silicon
    Ecut = 2

    ref_kcoords = [
        [0.000000000000, 0.000000000000, 0.000000000000],
        [0.038461538462, 0.000000000000, 0.038461538462],
        [0.076923076923, 0.000000000000, 0.076923076923],
        [0.115384615385, 0.000000000000, 0.115384615385],
        [0.153846153846, 0.000000000000, 0.153846153846],
        [0.192307692308, 0.000000000000, 0.192307692308],
        [0.230769230769, 0.000000000000, 0.230769230769],
        [0.269230769231, 0.000000000000, 0.269230769231],
        [0.307692307692, 0.000000000000, 0.307692307692],
        [0.346153846154, 0.000000000000, 0.346153846154],
        [0.384615384615, 0.000000000000, 0.384615384615],
        [0.423076923077, 0.000000000000, 0.423076923077],
        [0.461538461538, 0.000000000000, 0.461538461538],
        [0.500000000000, 0.000000000000, 0.500000000000],
        [0.500000000000, 0.000000000000, 0.500000000000],
        [0.531250000000, 0.062500000000, 0.531250000000],
        [0.562500000000, 0.125000000000, 0.562500000000],
        [0.593750000000, 0.187500000000, 0.593750000000],
        [0.625000000000, 0.250000000000, 0.625000000000],
        [0.375000000000, 0.375000000000, 0.750000000000],
        [0.348214285714, 0.348214285714, 0.696428571429],
        [0.321428571429, 0.321428571429, 0.642857142857],
        [0.294642857143, 0.294642857143, 0.589285714286],
        [0.267857142857, 0.267857142857, 0.535714285714],
        [0.241071428571, 0.241071428571, 0.482142857143],
        [0.214285714286, 0.214285714286, 0.428571428571],
        [0.187500000000, 0.187500000000, 0.375000000000],
        [0.160714285714, 0.160714285714, 0.321428571429],
        [0.133928571429, 0.133928571429, 0.267857142857],
        [0.107142857143, 0.107142857143, 0.214285714286],
        [0.080357142857, 0.080357142857, 0.160714285714],
        [0.053571428571, 0.053571428571, 0.107142857143],
        [0.026785714286, 0.026785714286, 0.053571428571],
        [0.000000000000, 0.000000000000, 0.000000000000],
        [0.000000000000, 0.000000000000, 0.000000000000],
        [0.041666666667, 0.041666666667, 0.041666666667],
        [0.083333333333, 0.083333333333, 0.083333333333],
        [0.125000000000, 0.125000000000, 0.125000000000],
        [0.166666666667, 0.166666666667, 0.166666666667],
        [0.208333333333, 0.208333333333, 0.208333333333],
        [0.250000000000, 0.250000000000, 0.250000000000],
        [0.291666666667, 0.291666666667, 0.291666666667],
        [0.333333333333, 0.333333333333, 0.333333333333],
        [0.375000000000, 0.375000000000, 0.375000000000],
        [0.416666666667, 0.416666666667, 0.416666666667],
        [0.458333333333, 0.458333333333, 0.458333333333],
        [0.500000000000, 0.500000000000, 0.500000000000],
        [0.500000000000, 0.500000000000, 0.500000000000],
        [0.500000000000, 0.472222222222, 0.527777777778],
        [0.500000000000, 0.444444444444, 0.555555555556],
        [0.500000000000, 0.416666666667, 0.583333333333],
        [0.500000000000, 0.388888888889, 0.611111111111],
        [0.500000000000, 0.361111111111, 0.638888888889],
        [0.500000000000, 0.333333333333, 0.666666666667],
        [0.500000000000, 0.305555555556, 0.694444444444],
        [0.500000000000, 0.277777777778, 0.722222222222],
        [0.500000000000, 0.250000000000, 0.750000000000],
        [0.500000000000, 0.250000000000, 0.750000000000],
        [0.500000000000, 0.208333333333, 0.708333333333],
        [0.500000000000, 0.166666666667, 0.666666666667],
        [0.500000000000, 0.125000000000, 0.625000000000],
        [0.500000000000, 0.083333333333, 0.583333333333],
        [0.500000000000, 0.041666666667, 0.541666666667],
        [0.500000000000, 0.000000000000, 0.500000000000],
    ]

    ref_klabels = Dict(
        "U"=>[0.625, 0.25, 0.625],
        "W"=>[0.5, 0.25, 0.75],
        "X"=>[0.5, 0.0, 0.5],
        "Γ"=>[0.0, 0.0, 0.0],
        "L"=>[0.5, 0.5, 0.5],
        "K"=>[0.375, 0.375, 0.75]
    )

    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    model = model_DFT(silicon.lattice, [spec => testcase.positions], [:lda_xc_teter93])
    kcoords, klabels, kpath = high_symmetry_kpath(model; kline_density=22.7)

    @test length(ref_kcoords) == length(kcoords)
    for ik in 1:length(ref_kcoords)
        @test ref_kcoords[ik] ≈ kcoords[ik] atol=1e-11
    end

    @test length(klabels) == length(ref_klabels)
    for key in keys(ref_klabels)
        @test klabels[key] ≈ ref_klabels[key] atol=1e-15
    end

    @test kpath[1] == ["Γ", "X", "U"]
    @test kpath[2] == ["K", "Γ", "L", "W", "X"]
end

@testset "High-symmetry kpath construction for 1D system" begin
    lattice = diagm([8.0, 0, 0])
    model = Model(lattice, n_electrons=1, terms=[Kinetic()])
    kcoords, klabels, kpath = high_symmetry_kpath(model; kline_density=20)

    @test length(kcoords) == 17
    @test kcoords[1]  ≈ [-1/2, 0, 0]
    @test kcoords[9]  ≈ [   0, 0, 0]
    @test kcoords[17] ≈ [ 1/2, 0, 0]
    @test length(kpath) == 1
end

@testset "Compute bands for silicon" begin
    testcase = silicon
    Ecut = 7
    n_bands = 8

    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    model = model_DFT(silicon.lattice, [spec => testcase.positions], :lda_xc_teter93)
    basis = PlaneWaveBasis(model, Ecut, testcase.kcoords, testcase.ksymops)

    # Build Hamiltonian just from SAD guess
    ρ0 = guess_density(basis, [spec => testcase.positions])
    ham = Hamiltonian(basis; ρ=ρ0)

    # Check that plain diagonalization and compute_bands agree
    eigres = diagonalize_all_kblocks(lobpcg_hyper, ham, n_bands + 3, n_conv_check=n_bands,
                                     tol=1e-5)

    band_data = compute_bands(basis, [k.coordinate for k in basis.kpoints]; ρ=ρ0, n_bands)
    for ik in 1:length(basis.kpoints)
        @test eigres.λ[ik][1:n_bands] ≈ band_data.λ[ik] atol=1e-5
    end
end

@testset "prepare_band_data" begin
    testcase = silicon
    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    model = model_DFT(silicon.lattice, [spec => testcase.positions], :lda_xc_teter93)

    # k coordinates simulating two band branches, Γ => X => W and U => X
    kcoords = [
        [0.000, 0.000, 0.000],
        [0.250, 0.000, 0.250],
        [0.500, 0.000, 0.500],
        #
        [0.500, 0.000, 0.500],
        [0.500, 0.125, 0.625],
        [0.500, 0.250, 0.750],
        #
        [0.625, 0.250, 0.625],
        [0.575, 0.150, 0.575],
        [0.500, 0.000, 0.500],
    ]
    ksymops = [[DFTK.identity_symop()] for _ in 1:length(kcoords)]
    basis = PlaneWaveBasis(model, 5, kcoords, ksymops)
    klabels = Dict("Γ" => [0, 0, 0], "X" => [0.5, 0.0, 0.5],
                   "W" => [0.5, 0.25, 0.75], "U" => [0.625, 0.25, 0.625])

    # Setup some dummy data
    λ = [10ik .+ collect(1:4) for ik = 1:length(kcoords)]  # Simulate 4 computed bands
    λerror = [λ[ik]./100 for ik = 1:length(kcoords)]       # ... and 4 errors

    ret = DFTK.prepare_band_data((basis=basis, λ=λ, λerror=λerror), klabels=klabels)

    @test ret.n_spin   == 1
    @test ret.n_kcoord == 9
    @test ret.n_bands  == 4

    @test ret.branches[1].kindices == [1, 2, 3]
    @test ret.branches[2].kindices == [4, 5, 6]
    @test ret.branches[3].kindices == [7, 8, 9]

    @test ret.branches[1].klabels == ("Γ", "X")
    @test ret.branches[2].klabels == ("X", "W")
    @test ret.branches[3].klabels == ("U", "X")

    for iband in 1:4
        @test ret.branches[1].λ[:, iband, 1] == [10ik .+ iband for ik in 1:3]
        @test ret.branches[2].λ[:, iband, 1] == [10ik .+ iband for ik in 4:6]
        @test ret.branches[3].λ[:, iband, 1] == [10ik .+ iband for ik in 7:9]

        for ibr in 1:3
            @test ret.branches[ibr].λerror[:, iband, 1] == ret.branches[ibr].λ[:, iband, 1] ./ 100
        end
    end

    B = model.recip_lattice
    ref_kdist = zeros(3, 3)  # row idx is k-point, col idx is branch,
    ikpt = 1
    for ibr in 1:3
        ibr != 1 && (ref_kdist[1, ibr] = ref_kdist[end, ibr-1])
        ikpt += 1
        for ik in 2:3
            ref_kdist[ik, ibr] = (
                ref_kdist[ik-1, ibr] + norm(B * (kcoords[ikpt-1] - kcoords[ikpt]))
            )
            ikpt += 1
        end
    end
    for ibr in 1:3
        @test ret.branches[ibr].kdistances == ref_kdist[:, ibr]
    end

    @test ret.ticks.labels == ["Γ", "X", "W | U", "X"]
    @test ret.ticks.distances == [0.0, ref_kdist[end, 1], ref_kdist[end, 2], ref_kdist[end, 3]]
end

@testset "is_metal" begin
    testcase = silicon
    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    model = model_LDA(silicon.lattice, [spec => testcase.positions])

    basis = PlaneWaveBasis(model, 5, testcase.kcoords, testcase.ksymops)
    λ = [[1, 2, 3, 4], [1, 1.5, 3.5, 4.2], [1, 1.1, 3.2, 4.3], [1, 2, 3.3, 4.1]]

    @test !DFTK.is_metal((λ=λ, basis=basis), 2.5)
    @test DFTK.is_metal((λ=λ, basis=basis), 3.2)
end
end


text:
include("friend_sets_struct.jl")
include("enemy_sets_struct.jl")
include("item_use_struct.jl")
include("irt_struct.jl")

mutable struct Settings
    n_items::Int64
    n_fs::Int64
    bank::DataFrames.DataFrame
    irt::IRT
    theta_bounds::Vector{Vector{Float64}}
    forced0::Vector{Vector{Bool}}
    n_groups::Int64
    T::Int64
    Tg::Vector{Int64}
    fs::FriendSets # friend Sets
    es::EnemySets # enemy Sets
    iu::ItemUse
    ol_max::Matrix{Float64}
    to_apply::Vector{Bool} #[iu_max?, iu_min?, ol?]
    Settings(
        n_items,
        n_fs,
        bank,
        irt,
        theta_bounds,
        forced0,
        n_groups,
        T,
        Tg,
        fs,
        es,
        iu,
        ol_max,
        to_apply,
    ) = new(
        n_items,
        n_fs,
        bank,
        irt,
        theta_bounds,
        forced0,
        n_groups,
        T,
        Tg,
        fs,
        es,
        iu,
        ol_max,
        to_apply,
    ) # no pattern mode
    Settings() = new(
        0,
        0,
        DataFrames.DataFrame(),
        IRT(),
        [[-6.0, 6.0]],
        Vector{Vector{Bool}}(undef, 0),
        1,
        1,
        [1],
        FriendSets(),
        EnemySets(),
        ItemUse(),
        zeros(Int64, 0, 0),
        [false, false, false],
    )
end


text:
import Statistics, DataFrames, StatsBase

function fix_variable_type(x::Vector)
  if all(isa.(x, Number))
    x = Vector{Number}(x)
  elseif all(isa.(x, String))
    if length(unique(x)) < 7
      x = DataFrames.categorical(x)
    else
      x = Vector{String}(x)
    end
  else
    @warn "Undetected variable type."
  end
  return x
end

function report(x::Vector{String}; n_strings::Int=4, kwargs...)
    # Sort by frequency
    sorted =first.(sort(collect(StatsBase.countmap(x)), by = i -> i[2], rev=true))
    if length(sorted) > n_strings
        sorted = "(" * join(sorted[1:n_strings], ", ") * ", ...)"
    else
        sorted = "(" * join(sorted, ", ", " and ") * ")"
    end

    text = "$(length(unique(x))) different entries " * sorted
    return Report(text=text)
end

function report(x::Vector{<:Any}; missing_percentage::Bool=true, kwargs...)
    n_missings = sum(ismissing.(x))
    if n_missings == 0
        x = fix_variable_type(x)
        text = report(x; kwargs...)
        return Report(text=text)
    end
    if missing_percentage == true
        n_missings = "$(round(n_missings / length(x), digits=2))%"
    end
    x = collect(skipmissing(x))
    text = report(x; kwargs...).text
    text *= ", missing = $n_missings"
    return Report(text=text)
end

function report(x::DataFrames.CategoricalVector; levels_percentage::Bool=true, kwargs...)
    text = string.(zeros(length(DataFrames.levels(x))))
    for (i, level) in enumerate(DataFrames.levels(x))
        if levels_percentage == true
            n = ", $(round(sum(x .== level)/length(x)*100, digits=2))%"
        else
            n = ", n = $(sum(x .== level))"
        end
        text[i] = level * n
    end
    text = join(text, "; ")

    return Report(text=text)
end

function report(x::Vector{<:Real}; median::Bool=false, dispersion::Bool=true, range::Bool=true, kwargs...)

    # Centrality
    if median == false
        var_centrality = "Mean = $(round(Statistics.mean(x), digits=2))"
    else
        var_centrality = "Median = $(round(Statistics.median(x), digits=2))"
    end

    # Dispersion
    if dispersion == true
        if median == true
            var_dispersion = " ± $(round(StatsBase.mad(x, normalize=true), digits=2))"
        else
            var_dispersion = " ± $(round(Statistics.std(x), digits=2))"
        end
    else
        var_dispersion = ""
    end

    # Range
    if range == true
        var_range = " [$(round(minimum(x), digits=2)), $(round(maximum(x), digits=2))]"
    else
        var_range = ""
    end

    text = var_centrality * var_dispersion * var_range
    text = format_text(text)

    return Report(text=text)
end