text:
function gale_shapley(proposing_group::AbstractVector{Agent}, other_group::AbstractVector{Agent}, rankings::RankedPreferences)
  rankings = copy(rankings)
  free_proposing_group = Dict{Agent, Bool}()
  matchings = Dict{Agent, Agent}() #keys correspond to agents in the other group!
  for i in 1:length(proposing_group)
    free_proposing_group[proposing_group[i]] = true
  end

  while length(free_proposing_group) > 0
    proposing_man = first(free_proposing_group)[1]
    first_woman = rankings.ranking[proposing_man][1]
    if haskey(matchings, first_woman)
      current_match = matchings[first_woman]
      if prefers(rankings, first_woman, proposing_man, current_match)
        matchings[first_woman] = proposing_man
        free_proposing_group[current_match] = true
        delete!(free_proposing_group, proposing_man)
      end
    else
      matchings[first_woman] = proposing_man
      delete!(free_proposing_group, proposing_man)
    end
    deleteat!(rankings.ranking[proposing_man], 1)
  end
  Match(matchings)
end

# group_1 is the proposing group
function extended_gale_shapley(proposing_group::AbstractVector{Agent}, other_group::AbstractVector{Agent})
  proposing_group = copy(proposing_group)
  other_group = copy(other_group)
  free_proposing_group = Dict{Agent, Bool}()
  free_other_group = Dict{Agent, Bool}()
  matchings = Dict{Agent, Agent}() #keys are from the other group!
  for i in 1:length(proposing_group)
    free_proposing_group[proposing_group[i]] = true
    free_other_group[other_group[i]] = true
    for j in 1:length(proposing_group)
      matchings[proposing_group[i]] = other_group[j]
    end
  end

  matchings_2 = Dict{Agent, Agent}()
  while length(free_proposing_group) > 0
    proposing_man = first(free_proposing_group)[1]
    first_woman = proposing_man.ranked_list[1]
    if haskey(matchings, first_woman)
      delete!(matchings, first_woman)
    end
    matchings[first_woman] = proposing_man
    first_woman
  end

end


text:
function warmup()
    t1 = now()
    ds = Dataset(x1 = rand(1:2, 10000), x2 = rand(1:2, 10000), x3 = rand(1:2, 10000) * 1.1, x4 = rand([Date(1), Date(2)], 10000), x5 = string.(rand(1:2), 10000),
                 x6 = PooledArray(Characters{5, UInt8}.(rand(1:2, 10000))), x7 = Int32.(rand(1:2, 10000)), x8 = Float32.(rand(1:2, 10000) * 1.1), x9 = PooledArray(Characters{3, UInt8}.(rand(1:2, 10000))), x10 = PooledArray(Characters{12, UInt8}.(rand(1:2, 10000))))
    for i in 1:ncol(ds)
        sortperm(ds, i)
    end
    for i in 1:ncol(ds)
        sortperm(ds, i, alg = QuickSort)
    end
    for i in 1:ncol(ds)
        sortperm(ds, 1:i)
    end
    for i in 1:ncol(ds)-1
        sortperm(ds, [ncol(ds), i])
    end
    for i in 1:ncol(ds)
        groupby(ds, i)
        gatherby(ds, i)
    end
    for i in 1:ncol(ds)
        groupby(ds, 1:i)
        gatherby(ds, 1:i)
    end
    byrow(ds, all, :, by = isequal(1))
    byrow(ds, sum)
    for op in (+, *, -, /)
        byrow(ds, op, 1:2)
    end
    combine(groupby(ds, [6,1]), (1,2)=>cor)
    combine(groupby(ds,1), Ref([1,2,3,7,8]) .=> [sum, mean, length, maximum, minimum, var, std, median, median!, sort])
    combine(groupby(ds,1), r"x1$" .=> [sum, mean, length, maximum, minimum, var, std, median, sort])
    ds2 = ds[1:2, [1,3]]
    combine(gatherby(ds,1), Ref([1,2,3,7,8]) .=> [median, sort])
    combine(gatherby(ds,1), Ref([1,2,3,7,8]) .=> [sum, mean, length, maximum, minimum, var, std])
    combine(gatherby(ds,1), r"x1$" .=> [sum, mean, length, maximum, minimum, var, std])
    IMD._gather_groups_hugeds_multicols(ds, 1:6, Val(Int32), threads = true)

    ds2 = ds[1:2, [1,3,7]]
    innerjoin(ds, ds2, on = [:x1, :x3, :x7])
    leftjoin(ds, ds2, on = [:x1, :x3, :x7])
    innerjoin(ds, ds2, on = [:x1, :x3, :x7], method = :hash)
    leftjoin(ds, ds2, on = [:x1, :x3, :x7], method = :hash)
    leftjoin(ds, ds2, on = [:x1, :x3, :x7], accelerate = true)
    ds3 = unique(ds, 1)[:, [1,3,5,7]]
    leftjoin(ds, ds3, on = :x1, makeunique = true)
    innerjoin(ds, ds3, on = :x1, makeunique = true)
    leftjoin(ds, ds3, on = :x1, makeunique = true, method = :hash)
    innerjoin(ds, ds3, on = :x1, makeunique = true, method = :hash)
    transpose(ds, 1:ncol(ds))
    transpose(groupby(ds,1:8), [2,3])
    show(IOBuffer(), MIME("text/html"), ds);
    # views
    ds = Dataset("A"=> ["A", "A" ,"A", "A", "B", "B","G"],
        "B"=> ["C", "D", "E", "B", "F", "G","N"]
        )
    select!(ds, 2,1)
    rename!(ds, [:child, :parent])

    function ff(ds, i)
        newds = leftjoin(ds, unique(dropmissing(ds,ncol(ds), view=true)[!, reverse(ncol(ds):-1:ncol(ds)-1)], [1,2] , view=true), on = [i+1=>1], makeunique = true)
        rename!(x->replace(x,"1"=>"grand"), newds)
        newds
    end

    ff(ff(ff(ds, 1),2), 3)

    ds = Dataset("A"=> ["A", "A" ,"A", "A", "B", "B","G"],
        "B"=> ["C", "D", "E", "B", "F", "G","N"]
        )
    select!(ds, 2,1)
    rename!(ds, [:child, :parent])

    function ff_hash(ds, i)
        newds = leftjoin(ds, unique(dropmissing(ds,ncol(ds), view=true)[!, reverse(ncol(ds):-1:ncol(ds)-1)], [1,2] , view=true), on = [i+1=>1], makeunique = true, method = :hash)
        rename!(x->replace(x,"1"=>"grand"), newds)
        newds
    end

    ff_hash(ff_hash(ff_hash(ds, 1),2), 3)
    @assert ff(ff(ff(ds, 1),2), 3) == ff_hash(ff_hash(ff_hash(ds, 1),2), 3)

    # some from test
    store = Dataset([[Date("2019-10-05"), Date("2019-10-04"), Date("2019-10-02"), Date("2020-01-01"), Date("2019-10-01"), Date("2019-10-02"), Date("2019-10-05"), Date("2019-10-04"), Date("2019-10-03"), Date("2019-10-03")],
        ["A", "A", "B", "A", "B", "A", "B", "B", "A", "B"]], [:date, :store])
    roster = Dataset([Union{Missing, String}["A", "A", "B", "A", "B", "A", "B", "B"],
             Union{Missing, Int64}[4, 1, 8, 2, 5, 3, 6, 7],
             Union{Missing, Date}[Date("2019-10-04"), Date("2019-09-30"), Date("2019-10-04"), Date("2019-10-02"), Date("2019-09-30"), Date("2019-10-03"), Date("2019-10-02"), Date("2019-10-03")],
             Union{Missing, Date}[Date("2019-10-06"), Date("2019-10-04"), Date("2019-10-06"), Date("2019-10-04"), Date("2019-10-04"), Date("2019-10-05"), Date("2019-10-04"), Date("2019-10-05")]],
             ["store", "employee_ID", "start_date", "end_date"])
    inn_r1 =  innerjoin(store, roster, on = [:date => (:start_date, nothing)], makeunique = true, stable = true)
    inn_r1_v =  innerjoin(store, view(roster, :, :), on = [:date => (:start_date, nothing)], makeunique = true, stable = true)
    inn_r1_a =  innerjoin(store, roster, on = [:date => (:start_date, nothing)], makeunique = true, stable = true, accelerate = true)
    inn_r1_v_a =  innerjoin(store, view(roster, :, :), on = [:date => (:start_date, nothing)], makeunique = true, stable = true, accelerate = true)

    inn_r1 =  innerjoin(store, roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true)
    inn_r1_a =  innerjoin(store, roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true, accelerate = true)
    inn_r1 =  innerjoin(store, roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true)
    inn_r1_a =  innerjoin(store, roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true, accelerate = true)

    inn_r1 =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true)
    inn_r1_a =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true, accelerate = true)
    inn_r1 =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true)
    inn_r1_a =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true, accelerate = true)

    inn_r1 =  innerjoin(store, roster, on = [:date => (:start_date, nothing)], makeunique = true, stable = true, method = :hash)
    inn_r1_v =  innerjoin(store, view(roster, :, :), on = [:date => (:start_date, nothing)], makeunique = true, stable = true, method = :hash)
    inn_r1_a =  innerjoin(store, roster, on = [:date => (:start_date, nothing)], makeunique = true, stable = true, accelerate = true, method = :hash)
    inn_r1_v_a =  innerjoin(store, view(roster, :, :), on = [:date => (:start_date, nothing)], makeunique = true, stable = true, accelerate = true, method = :hash)

    inn_r1 =  innerjoin(store, roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true, method = :hash)
    inn_r1_a =  innerjoin(store, roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true, accelerate = true, method = :hash)
    inn_r1 =  innerjoin(store, roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true, method = :hash)
    inn_r1_a =  innerjoin(store, roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true, accelerate = true, method = :hash)

    inn_r1 =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true, method = :hash)
    inn_r1_a =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (nothing, :start_date)], stable = true, accelerate = true, method = :hash)
    inn_r1 =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true, method = :hash)
    inn_r1_a =  innerjoin(view(store, :, :), roster, on = [:store => :store, :date => (:end_date, :start_date)], stable = true, accelerate = true, method = :hash)

    ds = Dataset(foo = ["one", "one", "one", "two", "two","two"],
                      bar = ['A', 'B', 'C', 'A', 'B', 'C'],
                      baz = [1, 2, 3, 4, 5, 6],
                      zoo = ['x', 'y', 'z', 'q', 'w', 't'])
    ds2 = transpose(groupby(ds, :foo, stable = true), :baz, id = :bar)
    ds = Dataset(g1 = ["g2", "g1", "g1", "g2", "g1", "g3"], c1 = ["c1", "c1", "c1", "c1", "c1", "c1"],
               id1 = ["id2_g2", "id1_g1", "id2_g1", "id1_g2", "id3_g1", "id4_g3"], id2 = [2, 3, 1, 5, 4, 6],
               val1 = [1, 2, 3, 4, 5, 6], val2 = [1.1, 2.1, 0.5, 1.3, 1.0, 2])
   ds2 = transpose(ds, [:val1, :val2])
   ds3 = transpose(ds, [:val1, :val2], id = :id2)
   ds4 = transpose(ds, [:val1, :val2], id = [:id1, :id2])
   ds5 = transpose(groupby(ds, 1, stable = true), [:val1, :val2])
   ds6 = transpose(groupby(ds, 1, stable = true), [:val1, :val2], id = :id1)
   ds7 = transpose(groupby(ds, 1, stable = true), [:val1, :val2], id = [:id1, :id2])
   ds8 = transpose(groupby(ds, 1:2, stable = true), [:val1, :val2], id = :id2)
   ds9 = transpose(groupby(ds, 1:2, stable = true), [:val1], id = [:id2, :id1])

   ds = Dataset([Union{Missing, Int64}[1, 1, 1, 2, 2, 2],
                Union{Missing, String}["foo", "bar", "monty", "foo", "bar", "monty"],
                Union{Missing, String}["a", "b", "c", "d", "e", "f"],
                Union{Missing, Int64}[1, 2, 3, 4, 5, 6]], [:g, :key, :foo, :bar])
    dst = transpose(groupby(ds, :g), (:foo, :bar), id = :key, variable_name = "_variables_")

    ds = Dataset(g = [1, 1, 1, 2, 2],
                        x1_int = [0, 0, 1, missing, 2],
                        x2_int = [3, 2, 1, 3, -2],
                        x1_float = [1.2, missing, -1.0, 2.3, 10],
                        x2_float = [missing, missing, 3.0, missing, missing],
                        x3_float = [missing, missing, -1.4, 3.0, -100.0])
    isequal(byrow(ds, argmin, r"float"), [:x1_float, missing, :x3_float, :x1_float, :x3_float])
    ds1 = Dataset(a = Union{String, Missing}["a", "b", "a", "b", "a", "b"],
                    b = Vector{Union{Int, Missing}}(1:6),
                    c = Union{Int, Missing}[1:3;1:3])
    ds = vcat(ds1, ds1)
    findall(duplicates(ds)) == collect(7:12)
    fmt(x)=1
    setformat!(ds, :a=>fmt)
    findall(duplicates(ds, :a, mapformats = true)) == 2:12
    unique(ds) == ds1
    unique(ds, 2:3) == ds1

    t2 = now()
    Dataset(x1 = "Finished warmup in", x2 = t2-t1)
end

ext: jl


text:
"""
	MM1K(λ, μ, k)

Tạo mô hình `M/M/1/K`.
"""
struct MM1K{T} <: AbstractMMCK
	λ::Union{T, Real}
	μ::Union{T, Real}
	k::Union{T, Integer}
	ρ::Union{T, Real}
	function MM1K(λ, μ, k)
		T = Union{typeof(μ), typeof(λ)}
		new{T}(λ, μ, k, λ/μ)
	end
end

function pn(m::MM1K, n::Int)
	ρ = m.ρ
	k = m.k
	if ρ == 1
		1 / (k + 1)
	else
		ρ^n * (1 - ρ) / (1 - ρ^(k+1))
	end
end

function L(m::MM1K)
	ρ = m.ρ
	k = m.k
	if ρ == 1
		k / 2
	else
		l = k + 1
		ρ * (1 + k*ρ^l - l*ρ^k) / (1 - ρ) / (1 - ρ^l)
	end
end

function Lq(m::MM1K)
	ρ = m.ρ
	k = m.k
	if ρ == 1
		k * (k - 1) / 2 / (k + 1)
	else
		L(m) - ρ * (1 - ρ^k) / (1 - ρ^(k+1))
	end
end

function W(m::MM1K)
	L(m) / λe(m)
end

function Wq(m::MM1K)
	W(m) - 1/m.μ
end


text:
"""
    TupleTabular{N}(i₁ => data₁, i₂ => data₂, ...)

Constructs an immutable `Tabular` data structure which stores data as a tuple of
`Pair`s between indices and elements (for `N == 1`) or nested containers (for
`N > 1`). Efficient for small numbers of indices (more than 14 is a bad idea).

This type is particularly useful for storing different data types (e.g. a
`Table` with columns of different element types). If the indices have distinct
types, then the compiler can deduce the return type from `getindex`, etc, at
compile-time. To ensure optimial compilation, ensure that `==` is a pure
function for the index values.
"""
struct TupleTabular{N, Data <: Tuple{Vararg{Pair}}} <: Tabular{N}
    data::Data

    function TupleTabular{N, Data}(d::Data) where {N, Data <: Tuple{Vararg{Pair}}}
        if N !== 1
            same_indices(last.(d))
        end
        new{N, Data}(d)
    end
end

@inline TupleTabular{N}(data::Data) where {N, Data <: Tuple{Vararg{Pair}}} = TupleTabular{N, Data}(data)
@inline TupleTabular{N}(data::Pair...) where {N} = TupleTabular{N}(data)

# Aliases for common sizes
const TupleTable{Data <: Tuple{Vararg{Pair}}} = TupleTabular{2, Data}
const TupleSeries{Data <: Tuple{Vararg{Pair}}} = TupleTabular{1, Data}

# indices
@inline indices(t::TupleTabular{1}) = (first.(t.data),)
@inline indices(t::TupleTabular) = (indices(t.data[1].second)..., first.(t.data))

# Get the tuple element
@inline get_tuple_value(data::Tuple{Vararg{Pair}}, i) = _get_tuple_value(i, data...)
@inline function _get_tuple_value(i, d::Pair)
    if i == d.first
        return d.second
    else
        throw(IndexError("Can't find index $i"))
    end
end
@inline function _get_tuple_value(i, d::Pair, ds::Pair...)
    if i == d.first
        return d.second
    else
        return _get_tuple_value(i, ds...)
    end
end

# Get the tuple pair
@inline get_tuple_pair(data::Tuple{Vararg{Pair}}, i) = _get_tuple_pair(i, data...)
@inline function _get_tuple_pair(i, d::Pair)
    if i == d.first
        return d
    else
        throw(IndexError("Can't find index $i"))
    end
end
@inline function _get_tuple_pair(i, d::Pair, ds::Pair...)
    if i == d.first
        return d
    else
        return _get_tuple_pair(i, ds...)
    end
end

# getindex
@propagate_inbounds function getindex(t::TupleSeries, i)
    get_tuple_value(t.data, i)
end

@propagate_inbounds function getindex(t::TupleSeries, ::Colon)
    t
end

@propagate_inbounds function getindex(t::TupleSeries, inds::Tuple)
    TupleSeries(_map(i -> get_tuple_pair(t.data, i), inds))
end

@propagate_inbounds function getindex(t::TupleTable, other_inds, i)
    get_tuple_value(t.data, i)[other_inds]
end

@propagate_inbounds function getindex(t::TupleTable, other_inds, ::Colon)
    data = _map(kv -> Pair(kv.first, kv.second[other_inds]), t.data)
    if aretype(Series, _map(last, data))
        TupleTable(data)
    else
        TupleSeries(data)
    end
end

@propagate_inbounds function getindex(t::TupleTable, other_inds, this_inds::Tuple)
    data = _map(i -> (kv = get_tuple_pair(t.data, i); Pair(kv.first, kv.second[other_inds])), this_inds)
    if aretype(Series, _map(last, data))
        TupleTable(data)
    else
        TupleSeries(data)
    end
end

# ===========
#  setindex!
# ===========

@propagate_inbounds function setindex!(s::TupleSeries, v, i)
    error("Cannot perform `setindex!` on an immutable object") # TODO choose an exception type... though this matches inbuilt Julia behavior...
end

@propagate_inbounds function setindex!(t::TupleTable, v, i1, i2)
    get_tuple_value(t.data, i2)[i1] = v
end

@propagate_inbounds function setindex!(t::TupleTable, v, i1, ::Colon)
    # TODO support v being a table or series.
    # TODO make faster for hetrogenously typed tuples.
    for kv ∈ t.data
        kv.second[i1] = v
    end
end

@propagate_inbounds function setindex!(t::TupleTable, v, i1, i2::Tuple)
    # TODO support v being a table or series.
    # TODO make faster for hetrogenously typed tuples.
    for i ∈ i2
        get_tuple_value(t.data, i)[i1] = v
    end
end

# ======
#  view
# ======

# TODO

#=
@propagate_inbounds function getindex(t::TupleTabular{N}, inds::Vararg{Any, N}) where {N}
    (other_inds, this_ind) = pop(inds)
    get_tuple_element(t.data, this_ind)[other_inds...]
end

# setindex!
@propagate_inbounds function setindex!(t::TupleTabular{1}, value, i)
    error("TupleSeries is immutable")
end

@propagate_inbounds function setindex!(t::TupleTabular{N}, value, inds::Vararg{Any, N}) where {N}
    (other_inds, this_ind) = pop(inds)
    t.dict[this_ind][other_inds...] = value
end
=#

# These optimizations are not safe...
function same_indices(iter::NTuple{N, T}) where {N, T <: TupleSeries}
    return true
end

function same_indices(iter::AbstractArray{T}) where {T <: TupleSeries}
    if isleaftype(T)
        return true
    else
        return false # should invoke more abstract method?
    end
end

function same_indices(iter::Associative{<:Any, T}) where {T <: TupleSeries}
    if isleaftype(T)
        return true
    else
        return false # should invoke more abstract method?
    end
end


text:
using DataFrames, StructuralEquationModels, Symbolics,
    LinearAlgebra, SparseArrays, Optim, LineSearches,
    BenchmarkTools, CSV, Statistics

cd("cfa")

include("functions.jl")

config = DataFrame(CSV.File("config.csv"))
config2 = copy(config)
config.backend .= "Optim.jl"
config2.backend .= "NLopt.jl"

config = [config; config2]

config = filter(row -> (row.backend == "Optim.jl") & (row.meanstructure == 0), config)

data_vec = read_files("data", get_data_paths(config))
par_vec = read_files("parest", get_data_paths(config))
# start_vec = read_files("start", get_data_paths(config))

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

models = gen_models(config, data_vec)

# fits = get_fits(models)

# correct = compare_estimates(fits, par_vec, config)

# config.correct = correct

############################################
# some benchmarks
############################################

using BenchmarkTools

semdiff = SemDiffOptim(
    LBFGS(
        ;linesearch = BackTracking(order=3),
        alphaguess = InitialHagerZhang()
        ),
    Optim.Options(
        ;f_tol = 1e-10,
        x_tol = 1.5e-8)
    )

ram_matrices_small = gen_CFA_RAM(3, 5)

ram_matrices_big = gen_CFA_RAM(5, 40)

model_ml_small = Sem(
    ram_matrices = ram_matrices_small,
    data = Matrix(data_vec[2]),
    imply = RAM,
    diff = semdiff
)

model_ml_big = Sem(
    ram_matrices = ram_matrices_big,
    data = Matrix(data_vec[6]),
    imply = RAM,
    diff = semdiff
)

############################################
# some benchmarks
############################################

bm_nonsymbolic_small = @benchmark sem_fit(model_ml_small)

bm_nonsymbolic_big = @benchmark sem_fit(model_ml_big)

# MKL

using MKL

bm_nonsymbolic_small_mkl = @benchmark sem_fit(model_ml_small)

bm_nonsymbolic_big_mkl = @benchmark sem_fit(model_ml_big)

# BFGS

algo = BFGS(;linesearch = BackTracking(order=3), alphaguess = InitialHagerZhang())

model_ml_small = Sem(
    ram_matrices = ram_matrices_small,
    data = Matrix(data_vec[2]),
    imply = RAM,
    diff = semdiff,
    algorithm = algo
)

model_ml_big = Sem(
    ram_matrices = ram_matrices_big,
    data = Matrix(data_vec[6]),
    imply = RAM,
    diff = semdiff,
    algorithm = algo
)

bm_algo_small = @benchmark sem_fit(model_ml_small)

bm_algo_big = @benchmark sem_fit(model_ml_big)

using ProfileView

function profile_test(model, n)
    for _ = 1:n
        sem_fit(model)
    end
end

ProfileView.@profview profile_test(model_ml_big, 10)


text:
module Wordlegames

using AbstractTrees
using DataFrames
using Random
using Tables

using AbstractTrees: print_tree

include("utilities.jl")
include("GamePool.jl")
include("trees.jl")

export GameNode,
    GamePool,
    GuessScore,
    GuessType,
    MinimizeExpected,
    MaximizeEntropy,
    Random,
    Tables,
    bincounts!,
    entropy2,
    expectedpoolsize,
    optimalguess,
    playgame!,
    print_tree,
    reset!,
    rowtable,
    scorecolumn!,
    scoreupdate!,
    showgame!,
    tiles,
    tree,
    updateguess!

end


text:
module PRS

include("GeneticMap.jl")
include("PlinkReader.jl")
include("LDMatrix.jl")
include("GaussianMixtureModel.jl")
include("utils.jl")
include("ldpred.jl")
include("lassosum.jl")

end # module


text:
using DelimitedFiles

# Generates a non-trivial random cone program
# min c'x
#  st Ax + s = b
#          s in K
# K is a convex cone (product of given cones)
function random_cone_program(dims, cone_prod)
    m,n = dims

    # Break into ortho parts in K and in K*
    z = randn(m)
    s_star = project_onto_cone(z, cone_prod)
    y_star = z - s_star
    A = sparse(randn(dims))
    x_star = randn(n)
    b = A*x_star + s_star
    c = -A' * y_star

    # TODO: return as a dict instead of tuple?
    return Dict(:A => A, :b => b, :c => c, :x_star => x_star,
                :y_star => y_star, :s_star => s_star)
    # return (A, b, c), (x_star, y_star, s_star)
end

function l1_minimization_program(dims; solve=false)
    m,n = dims
    m >= n || throw(ArgumentError("m < n"))
    A = randn(dims)
    b = randn(m)
    @assert rank(A) == n

    c_ = [zeros(n); ones(m)]
    A_ = -[
        -A                          Matrix{Float64}(I, m, m);
        A                           Matrix{Float64}(I, m, m);
        Matrix{Float64}(I, n, n)    zeros(n, m);
        ones(n)'                    zeros(m)'
    ]
    b_ = [b; -b; zeros(n); -1.0]
    # -Ax + b == s ∈ K
    cones = [MOI.Nonnegatives(2m+n), MOI.Zeros(1)]

    if solve
        model = Model()
        set_optimizer(model, optimizer_with_attributes(SCS.Optimizer,
            "eps" => 1e-10, "verbose" => 0))
        @variable(model, x[1:10])
        @variable(model, t[1:20])
        @objective(model, Min, sum(t))
        @constraint(model, 0 .<= A*x - b + t)
        @constraint(model, 0 .<= t - A*x + b)
        @constraint(model, x .>= 0)
        @constraint(model, sum(x) == 1.0)
        optimize!(model)
        return A_, b_, c_, cones, model
    else
        return A_, b_, c_, cones
    end

end

#
# A, b, c, cone_prod, model = l1_minimization_program((20, 10))
# m,n = 20,10
# xs = [value(x) for x in all_variables(model)][1:n]
# ts = value.(all_variables(model)[n+1:end])
# opt_val = sum(ts)
# sum(abs.(A[1:m,1:n]*xs - b[1:m]))
# all(A*[xs;ts]  - b .>= -1e10)
# (A*[xs;ts]  - b)[end] <= 1e10
# dot(c,[xs;ts])
#
# solve_opt_problem(A, b, c, cone_prod, nothing, SCS.Optimizer)
# # solve_opt_problem(A, b, c, [MOI.Nonnegatives(size(A)[1])], nothing, SCS.Optimizer)
# m,n = size(A)
# using Convex
# x = Variable(n)
# s = Variable(m)
# p = minimize(dot(x,c), A*x + s == b, s >= 0)
# solve!(p, () -> SCS.Optimizer())
#
# model = Model()
# optimizer_constructor = optimizer_with_attributes(SCS.Optimizer, "eps" => 1e-10, "verbose" => 1)
# set_optimizer(model, optimizer_constructor)
# @variable(model, x[1:size(A)[2]])
# @variable(model, s[1:size(A)[1]])
# @objective(model, Min, dot(c,x))
# @constraint(model, b - A*x .>= s)
# @constraint(model, s[end:end] in MOI.Zeros(1))
# @constraint(model, s[1:end-1] .>= 0)
# optimize!(model)
# optimizer_constructor
# MOI.get(owner_model(x[1]), MOI.VariablePrimal(1), x[1])
# typeof(backend(model)) <: SCS.Optimizer
# MOI.optimize!(backend(model))
#
# solution = SCS_solve(m, n, A, b, c, f=0, l=);


text:
using LTC
using BenchmarkTools
using Zygote

function test()

  function generate_data()
      in_features = 2
      out_features = 1
      N = 48
      data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
      data_x = [reshape([Float32(data_x[1][i]),Float32(data_x[2][i])],2,1) for i in 1:N]# |> f32
      data_y = [reshape([Float32(y)],1) for y in sin.(range(0,stop=6π,length=N))]# |> f32

      data_x = [repeat(x,1,20) for x in data_x]
      data_y = [repeat(x,1,20) for x in data_y]

      data_x, data_y
  end

  function lossf(x,y)
    Flux.reset!(ltc)
    ŷ = ltc.(x)
    sum(sum([(ŷ[i] .- y[i]) .^ 2 for i in 1:length(y)]))/length(y)#, ŷ
  end

  data_x,data_y = generate_data()

  ltc = NCP(Wiring(2,1))

  θ = Flux.params(ltc)

  @time gs = Zygote.gradient(θ) do
    lossf(data_x,data_y)#[1]
  end
  @time gs = Zygote.gradient(θ) do
    lossf(data_x,data_y)#[1]
  end
  @time gs = Zygote.gradient(θ) do
    lossf(data_x,data_y)#[1]
  end
end
@time test()


text:
module Simplify

# Base types
include("types.jl")

# Rule and context logic
include("rule.jl")
include("context.jl")

# Matching and unification
include("properties.jl")
include("match.jl")

# Default algebraic setup
include("rules.jl")
include("algebra.jl")

# Utility functions
include("utils.jl")

end # module


text:
include(joinpath(@__DIR__, "en", "standard_dictionary_numbers_extended.jl"))
include(joinpath(@__DIR__, "en", "large_standard_dictionary_numbers_extended.jl"))
include(joinpath(@__DIR__, "en", "ordinal_dictionaries.jl"))

# convert a value < 100 to English.
function small_convert_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    scale_numbers = _scale_modern # define scale type
    if number < 20
        word = _small_numbers[number + 1]

        return word
    end

    v = 0
    while v < length(_tens)
        d_cap = _tens[v + 1]
        d_number = BigInt(20 + 10 * v)

        if d_number + 10 > number
            if mod(number, 10) ≠ 0
                word = d_cap * "-" * _small_numbers[mod(number, 10) + 1]

                return word
            end

            return d_cap
        end
        v += 1
    end
end

# convert a value < 1000 to english, special cased because it is the level that excludes
# the < 100 special case.  The rest are more general.  This also allows you to get
# strings in the form of "forty-five hundred" if called directly.
function large_convert_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    scale_numbers = _scale_modern # define scale type
    word = string() # initiate empty string
    divisor = div(number, 100)
    modulus = mod(number, 100)

    if divisor > 0
        word = _small_numbers[divisor + 1] * " hundred"
        if modulus > 0
            word = word * " "
        end
    end

    if british
        if ! iszero(divisor) && ! iszero(modulus)
            word = word * "and "
        end
    end

    if modulus > 0
        word = word * small_convert_en(modulus, british=british, dict=dict)
    end

    return word
end

function spelled_out_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    scale_numbers = _scale_modern # default to :modern
    if isequal(dict, :british)
        scale_numbers = _scale_traditional_british
    elseif isequal(dict, :european)
        scale_numbers = _scale_traditional_european
    elseif isequal(dict, :modern)
    else
        throw(error("unrecognized dict value: $dict"))
    end

    number = big(number)
    isnegative = false
    if number < 0
        isnegative = true
    end

    number = abs(number)
    if number > limit - 1
        throw(error("""SpelledOut.jl does not support numbers larger than $(limit_str * " - 1").  Sorry about that!"""))
    end

    if number < 100
        word = small_convert_en(number, british=british, dict=dict)

        if isnegative
            word = "negative " * word
        end

        return word
    end

    if number < 1000
        word = large_convert_en(number, british=british, dict=dict)

        if isnegative
            word = "negative " * word
        end

        return word
    end

    v = 0
    while v ≤ length(scale_numbers)
        d_idx = v
        d_number = BigInt(round(big(1000)^v))

        if d_number > number
            modulus = BigInt(big(1000)^(d_idx - 1))
            l, r = divrem(number, modulus)
            word = large_convert_en(l, british=british, dict=dict) * " " * scale_numbers[d_idx - 1]

            if r > 0
                word = word * ", " * spelled_out_en(r, british=british, dict=dict)
            end

            if isnegative
                word = "negative " * word
            end

            return word
        end

        v += 1
    end
end

# Need to print ordinal numbers for the irrational printing
function spell_ordinal_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    s = spelled_out_en(number, british = british, dict = dict)

    lastword = split(s)[end]
    redolast = split(lastword, "-")[end]

    if redolast != lastword
        lastsplit = "-"
        word = redolast
    else
        lastsplit = " "
        word = lastword
    end

    firstpart = reverse(split(reverse(s), lastsplit, limit = 2)[end])
    firstpart = (firstpart == word) ? string() : firstpart * lastsplit

    if haskey(irregular, word)
        word = irregular[word]
    elseif word[end] == 'y'
        word = word[1:end-1] * ysuffix
    else
        word = word * suffix
    end

    return firstpart * word
end

# This method is an internal method used for spelling out floats
function decimal_convert_en(number::AbstractString; british::Bool = false, dict::Symbol = :modern)
    # decimal, whole = modf(number)
    # whole = round(BigInt, whole)
    whole, decimal = split(number, ".")
    word = spelled_out_en(parse(BigInt, whole), british=british, dict=dict) * string(" point")
    # word = spelled_out_en(whole, british=british, dict=dict) * string(" point")

    for i in decimal
        word = word * " " * _small_number_dictionary[i]
    end

    return word
end

function spelled_out_en(number::AbstractFloat; british::Bool = false, dict::Symbol = :modern)
    str_number = format(number)
    if occursin('.', str_number)
        _decimal, _ = modf(Dec64(number))
        _length_of_presicion = length(string(_decimal)) - 2 # (ndigits(_whole) + 1)
        number = format(number, precision = _length_of_presicion)
    else
        # It is an integer is scientific notation, treat normally without decimal precision considerations
        # E.g., 1e10 should be parsed as an integer (as should 1.0e10)
        number = parse(BigInt, str_number)
    end

    if isa(number, AbstractString)
        # if the number is a string then it is a decimal, which is formatted precisely
        # for correct precision with decimal places
        return decimal_convert_en(number, british = british, dict = dict)
    elseif isinteger(number)
        # otherwise, it is an integer
        return spelled_out_en(number, british = british, dict = dict)
    else
        throw(error("Cannot parse type $(typeof(number)).  Please make an issue."))
    end

    # should never get here
    return nothing
end

# Spell out complex numbers
function spelled_out_en(number::Complex; british::Bool = false, dict::Symbol = :modern)
    return spelled_out_en(real(number), british = british, dict = dict) * " and " * spelled_out_en(imag(number), british = british, dict=dict) * " imaginaries"
end

function spelled_out_en(number::Rational; british::Bool = false, dict::Symbol = :modern)
		_num, _den = number.num, number.den

        # return the number itself if the denomimator is one
		isone(_den) && return spelled_out_en(_num, british = british, dict = dict)

		word = spelled_out_en(_num, british = british, dict = dict) * " " * spell_ordinal_en(_den, british = british, dict = dict)

        # account for pluralisation
		return isone(_num) ? word : word * "s"
end

function spelled_out_en(number::AbstractIrrational; british::Bool = false, dict::Symbol = :modern)
    throw(error("Please round the input number, as we support floating point printing but cannot support irrational printing."))
end

# Fallback method if we do not know how to handle the input
function spelled_out_en(number; british::Bool = false, dict::Symbol = :modern)
    throw(error("Cannot parse type $(typeof(number)).  Please make an issue."))
end


text:
module CleverDicts

# The following two functions are overloaded for `MOI.VariableIndex` here
# it is the original use-case for `CleverDict`, and it would be type-piracy for
# solvers using `CleverDicts` to implement it themselves.

import MathOptInterface
import OrderedCollections

"""
    index_to_key(::Type{K}, index::Int)

Create a new key associated with the integer value `index`.
"""
function index_to_key end

"""
    key_to_index(key::K)

Map `key` to an integer valued index, assuming that there have been no
deletions.
"""
function key_to_index end

function index_to_key(::Type{MathOptInterface.VariableIndex}, index::Int64)
    return MathOptInterface.VariableIndex(index)
end

function index_to_key(
    ::Type{MathOptInterface.ConstraintIndex{F,S}},
    index::Int64,
) where {F,S}
    return MathOptInterface.ConstraintIndex{F,S}(index)
end

key_to_index(key::MathOptInterface.Index) = key.value

# Now, on with `CleverDicts`.

"""
    CleverDict{K, V}

A smart storage type for managing sequential objects with non-decreasing integer
indices. Note that the integers are `Int64`s.

Provided no keys are deleted, the backing storage is a `Vector{V}`. Once a key
has been deleted, the backing storage switches to an `OrderedDict{K, V}`.

Use the `add_item` to enforce adding item in sequence. Once an object is added
out of order `add_item` does not work anymore and the storage is switched to
`OrderedDict{K, V}` if it is not one already.

The i'th ordered element can be obtained with `c[LinearIndex(i)]`.

Note that querying a `LinearIndex` immediately after deleting a key via
`delete!` is very slow. (It requires a rebuild of an ordered list of variables.)

Store an item `val` using `add_item(c::CleverDict, val)`. `add_item` returns a
key corresponding to the stored item.

It is possible to initialize the CleverDict as `CleverDict{K, V}(n)` so that
`n` elements are pre-allocated and can be efficiently added even out of order
as long as their key hashes are between 1 and `n`.

Overload the functions `index_to_key` and `key_to_index` to enable mappings
between the integer index of the vector and the dictionary key.

## Example

```julia
struct MyKey
    x::Int64
end
index_to_key(::Type{MyKey}, i::Int64) = MyKey(i)
key_to_index(key::MyKey) = key.x
```
"""
mutable struct CleverDict{K,V,F<:Function,I<:Function} <: AbstractDict{K,V}
    last_index::Int64
    hash::F
    inverse_hash::I
    is_dense::Bool
    vector::Vector{V}
    dict::OrderedCollections.OrderedDict{K,V}
    function CleverDict{K,V}(
        hash::F = key_to_index,
        inverse_hash::I = index_to_key,
    ) where {K,V,F<:Function,I<:Function}
        return new{K,V,F,I}(
            0,
            hash,
            inverse_hash,
            true,
            K[],
            OrderedCollections.OrderedDict{K,V}(),
        )
    end
end

_is_dense(c::CleverDict) = c.is_dense

function _inverse_hash(c::CleverDict{K}, index::Integer) where {K}
    return c.inverse_hash(Int64(index))::K
end

function _inverse_hash(
    ::CleverDict{K,V,F,typeof(index_to_key)},
    index::Integer,
) where {K,V,F<:Function}
    return index_to_key(K, Int64(index))::K
end

"""
    add_item(c::CleverDict{K, V}, val::Val)::K where {K, V}

Set `val` in the next available key, and return that key.
"""
function add_item(c::CleverDict{K,V}, val::V)::K where {K,V}
    if c.last_index == -1
        error(
            "Keys were added out of order. `add_item` requires that keys are always added in order.",
        )
    end
    # adding a key in order
    key = _inverse_hash(c, c.last_index + 1)
    c[key] = val
    return key
end

function Base.empty!(c::CleverDict)
    empty!(c.vector)
    empty!(c.dict)
    c.last_index = 0
    c.is_dense = true
    return
end

function Base.length(c::CleverDict)
    return _is_dense(c) ? length(c.vector) : length(c.dict)
end

function Base.haskey(c::CleverDict{K}, key::K) where {K}
    if _is_dense(c)
        return 1 <= c.hash(key)::Int64 <= length(c.vector)
    end
    return haskey(c.dict, key)
end

function Base.keys(c::CleverDict{K}) where {K}
    if _is_dense(c)
        return K[_inverse_hash(c, i) for i in 1:length(c.vector)]
    end
    return collect(keys(c.dict))
end

function Base.get(c::CleverDict, key, default)
    if _is_dense(c)
        if !haskey(c, key)
            return default
        end
        return c.vector[c.hash(key)::Int64]
    end
    return get(c.dict, key, default)
end

function Base.getindex(c::CleverDict{K,V}, key::K) where {K,V}
    if _is_dense(c)
        if !haskey(c, key)
            throw(KeyError(key))
        end
        return c.vector[c.hash(key)::Int64]
    end
    return c.dict[key]
end

function Base.setindex!(c::CleverDict{K,V}, value::V, key::K)::V where {K,V}
    h = c.hash(key)::Int64
    if c.last_index != -1
        if h == c.last_index + 1
            c.last_index = h
        elseif h <= 0 || h > c.last_index
            c.last_index = -1
        end
    end
    if 1 <= h <= length(c.vector) && _is_dense(c)
        c.vector[h] = value
    elseif h == length(c.vector) + 1 && _is_dense(c)
        push!(c.vector, value)
    else
        if _is_dense(c)
            _rehash(c)
        end
        c.dict[key] = value
        # If there is a vector (e.g., because it has been rebuilt for
        # `LinearIndex`), clear it.
        if !isempty(c.vector)
            empty!(c.vector)
        end
    end
    return value
end

function _rehash(c::CleverDict{K}) where {K}
    sizehint!(c.dict, length(c.vector))
    @assert _is_dense(c)
    # Since c is currently dense, iterator protocol from CleverDict is used.
    for (k, v) in c
        c.dict[k] = v
    end
    empty!(c.vector)
    c.is_dense = false
    return
end

function Base.delete!(c::CleverDict{K}, k::K) where {K}
    if _is_dense(c)
        _rehash(c)
    end
    delete!(c.dict, k)
    if !isempty(c.vector)
        empty!(c.vector)
    end
    return c
end

struct LinearIndex
    i::Int64
end

function Base.getindex(c::CleverDict{K,V}, index::LinearIndex)::V where {K,V}
    if !(1 <= index.i <= length(c))
        throw(KeyError(index))
    end
    # Get the `index` linear element. If `c.vector` is currently `nothing`
    # (i.e., there has been a deletion), rebuild `c.vector`. This is a
    # trade-off: We could ensure `c.vector` is always updated, but this requires
    # a `splice!` in `delete!`, making deletions costly. However, it makes this
    # `getindex` operation trival because we would never have to rebuild the
    # vector.
    # The current implemented approach offers quick deletions, but an expensive
    # rebuild the first time you query a `LinearIndex` after a deletion or a new
    # key is added. Once the rebuild is done, there are quick queries until the
    # next deletion or addition. Thus, the worst-case is a user repeatedly
    # deleting a key and then querying a LinearIndex (e.g., getting the MOI
    # objective function).
    if !_is_dense(c) && length(c.dict) != length(c.vector)
        c.vector = Vector{V}(undef, length(c))
        for (i, val) in enumerate(values(c.dict))
            c.vector[i] = val
        end
    end
    return c.vector[index.i]::V
end

function Base.isempty(c::CleverDict)
    return _is_dense(c) ? isempty(c.vector) : isempty(c.dict)
end

Base.haskey(::CleverDict, key) = false

function Base.iterate(c::CleverDict{K,V}) where {K,V}
    if _is_dense(c)
        itr = iterate(c.vector)
        if itr === nothing
            return
        end
        return _inverse_hash(c, 1) => itr[1]::V, itr[2]::Int
    end
    itr = iterate(c.dict)
    if itr === nothing
        return
    end
    return convert(Pair{K,V}, itr[1]), itr[2]
end

function Base.iterate(c::CleverDict{K,V}, s::Int) where {K,V}
    if _is_dense(c)
        itr = iterate(c.vector, s)
        if itr === nothing
            return
        end
        return _inverse_hash(c, s) => itr[1]::V, itr[2]::Int
    end
    itr = iterate(c.dict, s)
    if itr === nothing
        return
    end
    return convert(Pair{K,V}, itr[1]), itr[2]
end

# we do not have to implement values and keys functions because they can rely
# on `Base`s fallback that uses the iterator protocol

function Base.sizehint!(c::CleverDict{K,V}, n) where {K,V}
    if _is_dense(c)
        sizehint!(c.vector, n)
    else
        sizehint!(c.dict, n)
    end
    return
end

function Base.resize!(c::CleverDict{K,V}, n) where {K,V}
    if _is_dense(c)
        if n < length(c.vector)
            error(
                "CleverDict cannot be resized to a size smaller than the current",
            )
        end
        resize!(c.vector, n)
    else
        sizehint!(c.dict, n)
    end
    return
end

Base.values(d::CleverDict) = _is_dense(d) ? d.vector : values(d.dict)

# TODO `map!(f, values(dict::AbstractDict))` requires Julia 1.2 or later,
#      use `map_values` once we drop Julia 1.1 and earlier.
function map_values!(f::Function, d::CleverDict)
    if _is_dense(d)
        map!(f, d.vector, d.vector)
    else
        for (k, v) in d.dict
            d.dict[k] = f(v)
        end
    end
    return
end

function Base.convert(
    ::Type{CleverDict{K,V,typeof(key_to_index),typeof(index_to_key)}},
    d::CleverDict{K,V,typeof(key_to_index),typeof(index_to_key)},
) where {K,V}
    return d
end

function Base.convert(
    ::Type{CleverDict{K,V,typeof(key_to_index),typeof(index_to_key)}},
    src::AbstractDict{K,V},
) where {K,V}
    dest = CleverDict{K,V}()
    for key in sort(collect(keys(src)), by = dest.hash)
        dest[key] = src[key]
    end
    return dest
end

end


text:
"""
Provides some useful functions for dealing with tree models and tree approximation using nested distance
Needs more ideas!

"""

function summary(tree0::Tree{A,B,C,D}) where {A,B,C,D}
    println()
    println("Tree with $(maximum(stage(tree0))+1) stages and branching structure $(tree0.name[6:end]) .")
    println("The tree has $(length(tree0.state)) nodes")
end

"""
	@author: Alois Pichler
	Created September 2018
		provides ontoSimplex: function to project onto the simplex
"""

function ontoSimplex!(probabilities::Vector{Float64})   # projects the probabilities onto the simplex
    for i= 1:length(probabilities)
        if probabilities[i] < 0.0   # make them nonnegative
            probabilities[i]= 0.0
        end
    end
	summ= sum(probabilities)
	if summ <= 0.0
		probabilities= fill!(probabilities, 1.0/ length(probabilities))
		summ= sum(probabilities)
		@info "ontoSimplex: created probabilities." maxlog=1
	end
	if summ != 1
		probabilities./= summ
		tmpi= rand(1:length(probabilities)) 	# modify a random index
		summ= -sum(probabilities[1:tmpi-1])+ 1.0- sum(probabilities[tmpi+1:end])
		probabilities[tmpi]= max(0, summ)
		@info "ontoSimplex: modified probabilities." maxlog=1
	end
    nothing     # modifies probabilities
end