text:
# MetidaBase.jl

struct MetidaTable{T <: NamedTuple}
    table::T
end
function metida_table(table::NamedTuple)
    MetidaTable(table)
end
function metida_table(args...; kwargs...)
    MetidaTable(metida_table_(args...; kwargs...))
end
function metida_table_(args...; names = nothing)
    if length(args) > 1
        e1 = length(args[1])
        i = 2
        @inbounds for i = 2:length(args)
            length(args[i]) == e1 || error("Length not equal")
        end
    end
    if isnothing(names)
        names = Tuple(Symbol.(:x , Symbol.(collect(1:length(args)))))
    else
        if length(args) != length(names) error("Length args and names not equal") end
        if !(typeof(names) <: Tuple)
            if !(typeof(names) <: AbstractVector{Symbol})
                names = Tuple(Symbol.(names))
            else
                names = Tuple(names)
            end
        end
    end
    NamedTuple{names}(args)
end

table(t::MetidaTable) = getfield(t, :table)
################################################################################
# TABLES
################################################################################
Tables.istable(t::MetidaTable) = true

Tables.columnaccess(t::MetidaTable) = true

Tables.columns(t::MetidaTable) = t

Tables.getcolumn(t::MetidaTable, i::Int) = getfield(t, :table)[i]

Tables.getcolumn(t::MetidaTable, nm::Symbol) = getfield(t, :table)[nm]

Tables.getcolumn(t::MetidaTable, ::Type{T}, col::Int, nm::Symbol) where {T} = t[:, col]

Tables.columnnames(t::MetidaTable) = names(t)

Tables.rowaccess(::Type{<:MetidaTable}) = true
# just return itself, which means MatrixTable must iterate `Tables.AbstractRow`-compatible objects
Tables.rows(t::MetidaTable) = t

# a custom row type; acts as a "view" into a row of an AbstractMatrix
struct MetidaTableRow{T} <: Tables.AbstractRow
    row::Int
    source::MetidaTable{T}
end

Base.iterate(t::MetidaTable, st=1) = st > length(t) ? nothing : (MetidaTableRow(st, t), st + 1)

Tables.getcolumn(t::MetidaTableRow, ::Type, col::Int, nm::Symbol) = getfield(t, :source)[getfield(t, :row), col]

Tables.getcolumn(t::MetidaTableRow, i::Int) = getfield(t, :source)[getfield(t, :row), i]

Tables.getcolumn(t::MetidaTableRow, nm::Symbol) = getfield(t, :source)[getfield(t, :row), nm]

Tables.columnnames(t::MetidaTableRow) = names(getfield(t, :source))

Tables.schema(t::MetidaTable) = Tables.Schema(names(t), eltype(y) for y in t.table)

################################################################################
# BASE
################################################################################

Base.names(t::MetidaTable) = collect(keys(t.table))

function Base.getindex(t::MetidaTable, col::Colon, ind::T) where T <: Union{Symbol, Int}
    Tables.getcolumn(t, ind)
end
function Base.getindex(t::MetidaTable, row::Int, ind::T) where T <: Union{Symbol, Int}
    Tables.getcolumn(t, ind)[row]
end

function Base.setindex!(t::MetidaTable, val, row::Int, ind::T) where T <: Union{Symbol, Int}
    Tables.getcolumn(t, ind)[row] = val
end

function Base.append!(t::MetidaTable, t2::MetidaTable)
    if !(names(t) ⊆ names(t2)) error("Names for t not in t2") end
    for n in names(t)
        append!(t[:, n], t2[:, n])
    end
    t
end

function Base.pushfirst!(t::MetidaTable, row::AbstractVector)
    if length(row) != length(keys(t.table)) error("Size not equal") end
    i = 1
    for i = 1:length(row)
        pushfirst!(t.table[i], row[i])
    end
    t
end
function Base.pushfirst!(t::MetidaTable, row::NamedTuple)
    kt = keys(t.table)
    kr = keys(row)
    if !issetequal(kt, kr) error("Size not equal") end
    for i in kt
        pushfirst!(t.table[i], row[i])
    end
    t
end

Base.length(t::MetidaTable) = length(first(t.table))

function Base.size(t::MetidaTable, i::Int)
    if i == 1
        return length(first(t.table))
    elseif i == 2
        return length(t.table)
    else
        error("Wrong dimention!")
    end
end

function Base.show(io::IO, table::MetidaTable)
    pretty_table(io, table; tf = PrettyTables.tf_compact)
end

# All
################################################################################
# DATASET
################################################################################
struct DataSet{T <: AbstractData} <: AbstractDataSet{AbstractData}
    ds::Vector{T}
end

function getdata(d::DataSet)
    d.ds
end
@inline function getindormiss(d::Dict{K}, i::K) where K
    ind::Int = ht_keyindex(d, i)
    if ind > 0 return d.vals[ind]  end
    missing
end
################################################################################
# BASE
################################################################################

function Base.getindex(d::DataSet, ind::Int)
    d.ds[ind]
end
@inline function getresultindex_safe(rd::T, ind::Symbol) where T <: AbstractResultData
    getindormiss(rd.result, ind)
end
@inline function getresultindex_unsafe(rd::T, ind::Symbol) where T <: AbstractResultData
    rd.result[ind]
end

function Base.getindex(d::DataSet{T}, col::Int, ind) where T <: AbstractResultData
    getresultindex_safe(d[col], ind)
end
function Base.getindex(d::DataSet{T}, col::Colon, ind) where T <: AbstractResultData
    @inbounds for i in Base.OneTo(length(d))
        if Base.ht_keyindex(d.ds[i].result, ind) < 1 return getresultindex_safe.(d.ds, ind) end
    end
    getresultindex_unsafe.(d.ds, ind)
end

Base.first(d::DataSet) = first(getdata(d))

function Base.length(d::DataSet)
    length(getdata(d))
end

function Base.iterate(d::DataSet)
    return Base.iterate(getdata(d))
end

function Base.iterate(d::DataSet, i::Int)
    return Base.iterate(getdata(d), i)
end

function Base.map(f, d::DataSet)
    DataSet(map(f, getdata(d)))
end
################################################################################
# BASE.SORT
################################################################################
function islessdict(a::Dict{A1,A2}, b::Dict{B1,B2}, k::Union{AbstractVector, Set}) where A1 where A2 where B1 where B2
    l = length(k)
    av = Vector{Union{Missing, A2}}(undef, l)
    bv = Vector{Union{Missing, B2}}(undef, l)
    @inbounds for i = 1:l
        av[i] = getindormiss(a, k[i])
        bv[i] = getindormiss(b, k[i])
    end
    isless(av, bv)
end
function islessdict(a::Dict, b::Dict, k)
    isless(getindormiss(a, k), getindormiss(b, k))
end
function Base.sort!(d::DataSet{T}, k; alg::Base.Algorithm = QuickSort, lt=nothing, by=nothing, rev::Bool=false, order::Base.Ordering = Base.Forward) where T <: Union{AbstractIdData, AbstractIDResult}
    if isnothing(by) by = x -> getid(x) end
    if isnothing(lt) lt = (x, y) -> islessdict(x, y, k) end
    sort!(d.ds;  alg = alg, lt = lt, by = by, rev = rev, order = order)
    d
end

################################################################################
# filter
################################################################################
function Base.filter(f::Function, d::DataSet)
    ds   =  getdata(d)
    inds  = findall(f, ds)
    DataSet(ds[inds])
end
function Base.filter!(f::Function, d::DataSet)
    filter!(f, getdata(d))
    d
end
################################################################################
# SELF
################################################################################

getid_safe(idd::AbstractIdData, ind) = getindormiss(idd.id, ind)

getid_unsafe(idd::AbstractIdData, ind) = idd.id[ind]

getid_safe(asr::AbstractIDResult, ind) = getindormiss(asr.data.id, ind)

getid_unsafe(asr::AbstractIDResult, ind) = asr.data.id[ind]

getid(idd::AbstractIdData, ind) = getid_safe(idd, ind)

getid(asr::AbstractIDResult, ind) = getid_safe(asr, ind)

getid(idd::AbstractIdData) = idd.id

getid(asr::AbstractIDResult) = asr.data.id

function getid(d::DataSet{T}, col::Int, ind) where T <: Union{AbstractIdData, AbstractIDResult}
    getid(d[col], ind)
end
function getid(d::DataSet{T}, col::Colon, ind) where T <: AbstractIdData
    @inbounds for i in Base.OneTo(length(d))
        if Base.ht_keyindex(d.ds[i].id, ind) < 1 return getid_safe.(d.ds, ind) end
    end
    getid_unsafe.(d.ds, ind)
end
function getid(d::DataSet{T}, col::Colon, ind) where T <: AbstractIDResult
    @inbounds for i in Base.OneTo(length(d))
        if Base.ht_keyindex(d.ds[i].data.id, ind) < 1 return getid_safe.(d.ds, ind) end
    end
    getid_unsafe.(d.ds, ind)
end

function uniqueidlist(d::DataSet{T}, list::AbstractVector{Symbol}) where T <: AbstractIdData
    dl = Vector{Dict}(undef, 0)
    for i in d
        if list ⊆ keys(getid(i))
            subd = Dict(k => getid(i)[k] for k in list)
            if subd ∉ dl push!(dl, subd) end
        end
    end
    dl
end

function uniqueidlist(d::DataSet{T}, list::Symbol) where T <: AbstractIdData
    dl = Vector{Dict}(undef, 0)
    for i in d
        if list in keys(getid(i))
            subd = Dict(list => getid(i)[list])
            if subd ∉ dl push!(dl, subd) end
        end
    end
    dl
end

function subset(d::DataSet{T}, sort::Dict) where T <: AbstractIdData
    inds = findall(x-> sort ⊆ getid(x), getdata(d))
    if length(inds) > 0 return DataSet(getdata(d)[inds]) end
    DataSet(Vector{T}(undef, 0))
end
function subset(d::DataSet{T}, sort::Dict) where T <: AbstractIDResult
    inds = findall(x-> sort ⊆ getid(x), getdata(d))
    if length(inds) > 0 return DataSet(getdata(d)[inds]) end
    DataSet(Vector{T}(undef, 0))
end
function subset(d::DataSet{T}, inds) where T
    DataSet(getdata(d)[inds])
end
################################################################################
# metida_table from DataSet{AbstractIDResult}
################################################################################
function metida_table_(obj::DataSet{RD}; order = nothing, results = nothing, ids = nothing) where RD <: AbstractIDResult
    idset  = Set(keys(first(obj).data.id))
    resset = Set(keys(first(obj).result))
    if length(obj) > 1
        for i = 2:length(obj)
            union!(idset,  Set(keys(obj[i].data.id)))
            union!(resset, Set(keys(obj[i].result)))
        end
    end
    if !isnothing(results)
        if isa(results, Symbol) results = [results] end
        if isa(results, String) results = [Symbol(results)] end
        ressetl = isnothing(order) ? collect(intersect(resset, results)) : sortbyvec!(collect(intersect(resset, results)), order)
    else
        ressetl = isnothing(order) ? collect(resset) : sortbyvec!(collect(resset), order)
    end
    if !isnothing(ids)
        if isa(ids, Symbol) ids = [ids] end
        if isa(ids, String) ids = [Symbol(ids)] end
        ids ⊆ idset || error("Some id not in dataset!")
        idset = intersect(idset, ids)
    end
    mt1 = metida_table_((getid(obj, :, c) for c in idset)...; names = idset)
    mt2 = metida_table_((obj[:, c] for c in ressetl)...; names = ressetl)
    merge(mt1, mt2)
end
################################################################################
# TypedTables.jl interface

function TypedTables.Table(obj::AbstractDataSet; kwargs...)
    TypedTables.Table(metida_table_(obj; kwargs...))
end
function TypedTables.Table(obj::MetidaTable)
    TypedTables.Table(obj.table)
end

# MetidaFreq.jl
struct Proportion <: AbstractData
    x::Int
    n::Int
end


text:
using Documenter, FeatherLib

makedocs(
	modules = [FeatherLib],
	sitename = "FeatherLib.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/FeatherLib.jl.git"
)

ext: jl


text:
using SimpleWeightedGraphs: AbstractSimpleWeightedGraph, nv

function GCNConv(g::AbstractSimpleWeightedGraph, ch::Pair{<:Integer,<:Integer}, σ = identity;
                 init = glorot_uniform, T::DataType=Float32, bias::Bool=true)
    N = nv(g)
    b = bias ? init(ch[2], N) : zeros(T, ch[2], N)
    adj = adjacency_matrix(g)
    GCNConv(init(ch[2], ch[1]), b, normalized_laplacian(adj+I, T), σ)
end

function ChebConv(g::AbstractSimpleWeightedGraph, ch::Pair{<:Integer,<:Integer}, k::Integer;
                  init = glorot_uniform, T::DataType=Float32, bias::Bool=true)
    N = nv(g)
    b = bias ? init(ch[2], N) : zeros(T, ch[2], N)
    adj = adjacency_matrix(g)
    L̃ = T(2. / eigmax(Matrix(adj))) * normalized_laplacian(adj, T) - I
    ChebConv(init(ch[2], ch[1], k), b, L̃, k, ch[1], ch[2])
end

function GraphConv(g::AbstractSimpleWeightedGraph, ch::Pair{<:Integer,<:Integer}, aggr=:add;
                   init = glorot_uniform, bias::Bool=true)
    N = nv(g)
    b = bias ? init(ch[2], N) : zeros(T, ch[2], N)
    GraphConv(adjlist(g), init(ch[2], ch[1]), init(ch[2], ch[1]), b, aggr)
end

function GATConv(g::AbstractSimpleWeightedGraph, ch::Pair{<:Integer,<:Integer}; heads=1,
                 concat=true, negative_slope=0.2, init=glorot_uniform, bias::Bool=true)
    N = nv(g)
    b = bias ? init(ch[2], N) : zeros(T, ch[2], N)
    GATConv(adjlist(g), init(ch[2], ch[1]), b, init(2 * ch[2]), negative_slope)
end

function GatedGraphConv(g::AbstractSimpleWeightedGraph, out_ch::Integer, num_layers::Integer;
                        aggr=:add, init=glorot_uniform)
    N = nv(g)
    w = init(out_ch, out_ch, num_layers)
    gru = GRUCell(out_ch, out_ch)
    GatedGraphConv(adjlist(g), w, gru, out_ch, num_layers, aggr)
end

function EdgeConv(g::AbstractSimpleWeightedGraph, nn; aggr::Symbol=:max)
    EdgeConv(adjlist(g), nn, aggr)
end


text:
using SpecialFunctions

function map_log_gamma(v::Vector{N}) where {N <: Number}
  return map(loggamma, v)
end


text:
using Test
using Printf
using Dubins3D

include("instances3D.jl")

function print_results(data, name)
    @testset "Instances: $(name)" begin
        for (i, (qi, qf)) in enumerate(data)
            dubins = DubinsManeuver3D(qi, qf, rhomin, pitchmax)

            t = @timed DubinsManeuver3D(qi, qf, rhomin, pitchmax)

            lb = getLowerBound(qi, qf, rhomin, pitchmax)
            ub = getUpperBound(qi, qf, rhomin, pitchmax)
            best_lb = lb.length

            @test lb.length <= dubins.length + 1e-5
            @test dubins.length <= ub.length + 1e-5

            print(join([
        "$(@sprintf("%6s %d", name, i))",
        "$(@sprintf(" L:%9.2f", dubins.length))",
        #"$(@sprintf(" T: %5.3f", (t[2]) * 1000))",
        "$(@sprintf(" LB: %8.2f", lb.length))",
        "$(@sprintf(" UB: %8.2f", ub.length))",
        "$(@sprintf(" GAP[%%]: %6.2f", 100.0 * (dubins.length - best_lb) / dubins.length))",
        "\n",
      ]," "))
        end
    end
end

print_results(LONG, "Long")
print_results(SHORT, "Short")
print_results(BEARD, "Beard")
print_results(OTHERS, "Others")
print_results(HOTA, "Hota")


text:
"""
Test InductionMotor model:
This case study defines a four bus system with an infinite bus in 1,
a GENSAL in bus 2, and the (5th order model) induction motor in bus 3
The GENSAL machine has a SEXS and a HYGOV.
The disturbance is the outage of one line between buses 1 and 4
"""
##################################################
############### LOAD DATA ########################
##################################################

raw_file = joinpath(TEST_FILES_DIR, "data_tests/TVC_System_motor.raw")
dyr_file = joinpath(TEST_FILES_DIR, "data_tests/TVC_System_motor.dyr")

##################################################
############### SOLVE PROBLEM ####################
##################################################

@testset "Test 37 5th_order Ind. Motor" begin
    path = (joinpath(pwd(), "test-ind-mot_5th"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)
        time_span = (0.0, 20.0)
        perturbation_trip = BranchTrip(1.0, Line, "BUS1-BUS4-i_1")
        # Create Simulation with Constant Power
        sim_P = Simulation(ResidualModel, sys, path, time_span, perturbation_trip)

        # Motor parameters
        load = first(get_components(PSY.PowerLoad, sys))
        # Include the induction motor
        dynamic_injector = Ind_Motor(load)
        set_dynamic_injector!(load, dynamic_injector)
        sim = Simulation(ResidualModel, sys, path, time_span, perturbation_trip)

        # Test initial voltages between two systems are equivalent
        voltages_P = sim_P.x0_init[1:8]
        voltages_motor = sim.x0_init[1:8]
        @test LinearAlgebra.norm(voltages_P - voltages_motor) < 1e-4

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test37_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test37_eigvals) < 1e-3

        # Solve problem
        @test execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # No comparison since model is not available in other tools
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 37 5th_order Ind. Motor" begin
    path = (joinpath(pwd(), "test-ind-mot_5th"))
    !isdir(path) && mkdir(path)
    try
        sys = System(raw_file, dyr_file)
        time_span = (0.0, 20.0)
        perturbation_trip = BranchTrip(1.0, Line, "BUS1-BUS4-i_1")
        # Create Simulation with Constant Power
        sim_P = Simulation(MassMatrixModel, sys, path, time_span, perturbation_trip)

        # Motor parameters
        load = first(get_components(PSY.PowerLoad, sys))
        # Include the induction motor
        dynamic_injector = Ind_Motor(load)
        set_dynamic_injector!(load, dynamic_injector)
        sim = Simulation(MassMatrixModel, sys, path, time_span, perturbation_trip)

        # Test initial voltages between two systems are equivalent
        voltages_P = sim_P.x0_init[1:8]
        voltages_motor = sim.x0_init[1:8]
        @test LinearAlgebra.norm(voltages_P - voltages_motor) < 1e-4

        # Test Initial Condition
        diff = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test37_x0_init
            diff[1] += LinearAlgebra.norm(res[k] - v)
        end

        @test (diff[1] < 1e-3)

        # Obtain small signal results for initial conditions
        small_sig = small_signal_analysis(sim)
        eigs = small_sig.eigenvalues
        @test small_sig.stable

        #Test Eigenvalues
        @test LinearAlgebra.norm(eigs - test37_eigvals) < 1e-3

        # Solve problem
        @test execute!(sim, Rodas5(), dtmax = 0.005, saveat = 0.005) ==
              PSID.SIMULATION_FINALIZED
        results = read_results(sim)

        # No comparison since model is not available in other tools
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
## indexing ##

length(t::Tuple) = tuplelen(t)
endof(t::Tuple) = tuplelen(t)
size(t::Tuple, d) = d==1 ? tuplelen(t) : error("invalid tuple dimension $(d)")
getindex(t::Tuple, i::Int) = tupleref(t, i)
getindex(t::Tuple, i::Real) = tupleref(t, convert(Int, i))
getindex(t::Tuple, r::AbstractArray) = tuple([t[ri] for ri in r]...)
getindex(t::Tuple, b::AbstractArray{Bool}) = getindex(t,find(b))

## iterating ##

start(t::Tuple) = 1
done(t::Tuple, i::Int) = (length(t) < i)
next(t::Tuple, i::Int) = (t[i], i+1)

# this allows partial evaluation of bounded sequences of next() calls on tuples,
# while reducing to plain next() for arbitrary iterables.
indexed_next(t::Tuple, i::Int, state) = (t[i], i+1)
indexed_next(a::Array, i::Int, state) = (a[i], i+1)
indexed_next(I, i, state) = done(I,state) ? throw(BoundsError()) : next(I, state)

# eltype

eltype{T}(x::(T...)) = T

## mapping ##

ntuple(n::Integer, f::Function) = ntuple(f, n) # TODO: deprecate this?
ntuple(f::Function, n::Integer) =
    n<=0 ? () :
    n==1 ? (f(1),) :
    n==2 ? (f(1),f(2),) :
    n==3 ? (f(1),f(2),f(3),) :
    n==4 ? (f(1),f(2),f(3),f(4),) :
    n==5 ? (f(1),f(2),f(3),f(4),f(5),) :
    tuple(ntuple(n-2,f)..., f(n-1), f(n))

argtail(x, rest...) = rest
tail(x::Tuple) = argtail(x...)

# 0 argument function
map(f::Callable) = f()
# 1 argument function
map(f::Callable, t::())                   = ()
map(f::Callable, t::(Any,))               = (f(t[1]),)
map(f::Callable, t::(Any, Any))           = (f(t[1]), f(t[2]))
map(f::Callable, t::(Any, Any, Any))      = (f(t[1]), f(t[2]), f(t[3]))
map(f::Callable, t::(Any, Any, Any, Any)) = (f(t[1]), f(t[2]), f(t[3]), f(t[4]))
map(f::Callable, t::Tuple)                = tuple(f(t[1]), map(f,tail(t))...)
# 2 argument function
map(f::Callable, t::(),        s::())        = ()
map(f::Callable, t::(Any,),    s::(Any,))    = (f(t[1],s[1]),)
map(f::Callable, t::(Any,Any), s::(Any,Any)) = (f(t[1],s[1]), f(t[2],s[2]))
map(f::Callable, t::(Any,Any,Any), s::(Any,Any,Any)) =
    (f(t[1],s[1]), f(t[2],s[2]), f(t[3],s[3]))
map(f::Callable, t::(Any,Any,Any,Any), s::(Any,Any,Any,Any)) =
    (f(t[1],s[1]), f(t[2],s[2]), f(t[3],s[3]), f(t[4],s[4]))
# n argument function
map(f::Callable, ts::Tuple...) = tuple([f(map(t->t[n],ts)...) for n=1:length_checked_equal(ts...)]...)

function length_checked_equal(args...)
    n = length(args[1])
    for i=2:length(args)
        if length(args[i]) != n
            error("argument dimensions must match")
        end
    end
    n
end

## comparison ##

function isequal(t1::Tuple, t2::Tuple)
    if length(t1) != length(t2)
        return false
    end
    for i = 1:length(t1)
        if !isequal(t1[i], t2[i])
            return false
        end
    end
    return true
end

function ==(t1::Tuple, t2::Tuple)
    if length(t1) != length(t2)
        return false
    end
    for i = 1:length(t1)
        if !(t1[i] == t2[i])
            return false
        end
    end
    return true
end

function isless(t1::Tuple, t2::Tuple)
    n1, n2 = length(t1), length(t2)
    for i = 1:min(n1, n2)
        a, b = t1[i], t2[i]
        if !isequal(a, b)
            return isless(a, b)
        end
    end
    return n1 < n2
end

## functions ##

isempty(x::()) = true
isempty(x::Tuple) = false

revargs() = ()
revargs(x, r...) = tuple(revargs(r...)..., x)

reverse(t::Tuple) = revargs(t...)

## specialized reduction ##

# TODO: these definitions cannot yet be combined, since +(x...)
# where x might be any tuple matches too many methods.
sum(x::(Any, Any...)) = +(x...)

# NOTE: should remove, but often used on array sizes
prod(x::()) = 1
prod(x::(Any, Any...)) = *(x...)

all(x::()) = true
all(x::(Any, Any...)) = (&)(x...)

any(x::()) = false
any(x::(Any, Any...)) = |(x...)

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
import CSV

struct Account
    name::String
end

account_main = Account("Main")
account_net = Account("Net")

struct Transaction
    date::Date
    category::String
    countername::String
    description::String
    change::Dict{Account,Float64}
    balance::Dict{Account,Float64}
end

max_stars_repo_name: CrazedProgrammer/BalancePlot

max_issues_repo_name: CrazedProgrammer/BalancePlot

max_forks_repo_name: CrazedProgrammer/BalancePlot


text:
X = vec([+1 -1 -1 -1 -1 +1;
         -1 +1 -1 -1 +1 -1;
         -1 -1 +1 +1 -1 -1;
         -1 -1 +1 +1 -1 -1;
         -1 -1 +1 +1 -1 -1;
         -1 +1 -1 -1 +1 -1;
         +1 -1 -1 -1 -1 +1.])

O = vec([+1 +1 +1 +1 +1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 -1 -1 -1 -1 +1;
         +1 +1 +1 +1 +1 +1.])

F1 = vec([-1 -1 -1 -1 -1 -1;
          -1 +1 +1 +1 +1 -1;
          -1 +1 -1 -1 -1 -1;
          -1 +1 +1 +1 +1 -1;
          -1 +1 -1 -1 -1 -1;
          -1 +1 -1 -1 -1 -1;
          -1 -1 -1 -1 -1 -1.])

F2 = vec([+1 +1 +1 +1 +1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 +1 +1 +1 +1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 -1 -1 -1 -1 -1;
          +1 -1 -1 -1 -1 -1.])


text:
const MAX_CONTRACTION = 12
export MAX_CONTRACTION

#                        | s | p |    d     |    f    |
#                        ------------------------------
const axial_norm_fact = Matrix{Float64}(
                          [ 1.0 1.0    1.0        1.0   ;
                            0.0 1.0 sqrt(3.0)  sqrt(5.0);
                            0.0 1.0 sqrt(3.0)  sqrt(5.0);
                            0.0 0.0    1.0     sqrt(5.0);
                            0.0 0.0 sqrt(3.0) sqrt(15.0);
                            0.0 0.0    1.0     sqrt(5.0);
                            0.0 0.0    0.0        1.0   ;
                            0.0 0.0    0.0     sqrt(5.0);
                            0.0 0.0    0.0     sqrt(5.0);
                            0.0 0.0    0.0        1.0
                          ]
                        )
export axial_norm_fact


text:
include("../src/frequency_adaptation.jl")

# Combine Loschilov's ACD with Glasmachers ACF-CD. Initial testing indicated it is 3-20 times
# worse than just ACD. I guess the ACD needs new fitness evals in each direction of the transformation
# matrix in order to reliably update it. So maybe only use ACF when increasing the population size, i.e.
# ensure that all directions have been tried at least once and then add additional once from the frequency
# selected directions?
function acf_adaptive_coordinate_descent(fitnessfct, P, Xmin, Xmax;
  MaxEval = 1e4 * P,
  stopfitness = 1e-8,
  howOftenUpdateRotation = 1,
  numfevals = 0,
  kSuccess = 2.0,
  kUnsuccess = 0.5,
  alwaysSampleTwoPoints = false, # sample 2 points per coordinate regardless if 1st is improvement
  nonImprovementBudget = 1       # Percent of fevals we allow before restarting if not enough improvement
  )

  # Initialize mean and sigma
  xmins = Xmin * ones(P)
  xmaxs = Xmax * ones(P)
  xmean = xmins .+ rand(P) * (Xmax - Xmin)
  sigma = (xmaxs .- xmins) / 4 # Step size per principal component (NOT per dimension!)

  bestFit = 1e+100 # Start from a very high fitness to ensure we will find lower... Fix smarter...

  mu = P # Use the P best candidates among the 2*P in each update of the adaptive encoding...
  ae = AdaptiveEncoding(mu, P)
  # To compare to CD without adaptive encoding instead use:
  #ae = NoEncoding(mu, P)
  #alwaysSampleTwoPoints = true

  iter = 0
  somebetter = false

  # There can be 2*P or 2*P+1 candidate points tested before we update the encoding,
  # since we only add the 2nd point if the first was not better.
  allx = zeros(P, 2*P+1) # Each column is a candidate and we will save 2 candidates per component/iteration in the population
  allf = zeros(2*P+1)    # One fitness value per candidate
  currentpos = 0

  # Restart the search from a new random position unless the fitness has improved more than 2*stopfitness
  # in last nonImprovementBudget % of fevals.
  lastconvergencecheck = numfevals
  lastcheckbestfit = bestFit

  # We will adaptively select among the P directions
  freqadapter = FrequencyAdapter(P)

  while ((numfevals < MaxEval) && (bestFit > stopfitness))
    ix = next(freqadapter) # Get the next direction according to the adapted frequencies

    # Sample two candidate solutions and calc their fitness
    dx = sigma[ix] * ae.B[:,ix] # shift along ix'th principal component, the computational complexity is linear

    # Randomly take the minus or plus direction first.
    direction = (rand() < 0.5) ? 1 : (-1)

    x1 = xmean + direction * dx             # first point to test along ix'th principal component
    Fit1 = fitnessfct(x1)
    update!(freqadapter, ix, (bestFit - Fit1) / bestFit)
    numfevals += 1

    # Update archive with new point
    # posx1 = ix*2 - 1
    currentpos += 1
    allx[:, currentpos] = x1
    allf[currentpos] = Fit1

    # Found better point?
    foundbetter = false
    if (Fit1 < bestFit)
      bestFit = Fit1
      fevalsatbest = numfevals
      xmean = x1
      foundbetter = true
    end

    # Only do second point if first point was not better (unless parameter states we should always smaple two).
    if !foundbetter || alwaysSampleTwoPoints
      x2 = xmean - direction * dx             # second point to test is symmetric to the first one on the same ix'principal component
      Fit2 = fitnessfct(x2)
      update!(freqadapter, ix, (bestFit - Fit2) / bestFit)
      numfevals += 1
      # Update archive with new point
      currentpos += 1
      # posx2 = ix*2
      allx[:, currentpos] = x2
      allf[currentpos] = Fit2

      if (Fit2 < bestFit)
        bestFit = Fit2
        fevalsatbest = numfevals
        xmean = x2
        foundbetter = true
      end
    end

    # Adapt step-size sigma depending on the success/unsuccess of the last search step.
    # Increase step size if found a better candidate, decrease if not.
    if foundbetter
      somebetter = true
      sigma[ix] *= kSuccess
    else
      sigma[ix] *= kUnsuccess
    end

    # Update the encoding every time we have tested at least 2*P candidate points and at least one
    # better candidate found.
    if currentpos >= 2*P
      if somebetter
        somebetter = false
        fitnessperm = sortperm(allf[1:currentpos])
        population = allx[:,fitnessperm]
        update!(ae, population, howOftenUpdateRotation)
      end

      # Check for convergence/nonimprovement every now and then and restart if so.
      if (numfevals - lastconvergencecheck) >= (MaxEval * nonImprovementBudget / 100.0)
        # We make a simple prediction of the final fitness given the slope between last
        # check and this. Restart if the predicted is too far from stopfitness.
        pfitness = predicted_final_fitness(numfevals, lastconvergencecheck, MaxEval, bestFit, lastcheckbestfit)
        if pfitness > 1000*stopfitness
          println("Restart!")
          # Recursive call to restart but express how many function evals we have already done...
          xm2, bf2, nf2 = acf_adaptive_coordinate_descent(fitnessfct, P, Xmin, Xmax;
            MaxEval = MaxEval, stopfitness = stopfitness,
            howOftenUpdateRotation = howOftenUpdateRotation,
            numfevals = numfevals)
          if bf2 < bestFit
            return (xm2, bf2, nf2)
          else
            return (xmean, bestFit, nf2)
          end
        end
        lastcheckbestfit = bestFit
        lastconvergencecheck = numfevals
      end

      currentpos = 0
    end

    iter += 1
    if rand() < 0.01
      pf = predicted_final_fitness(numfevals, lastconvergencecheck, MaxEval, bestFit, lastcheckbestfit)
      println("$(numfevals): bestfit = $(bestFit) (predicted: $(pf))")
    end
  end

  return (xmean, bestFit, numfevals)
end

function predicted_final_fitness(currentfevals, lastfevals, maxfevals, currentfitness, lastfitness)
  slope = (currentfitness - lastfitness) / (currentfevals - lastfevals)
  currentfitness + slope * (maxfevals - currentfevals)
end

function generate_random_orthogonal_transformation(P)
  B = eye(P,P)
  for i = 1:P
    v = randn(P,1)
    while norm(v) < 1e-2
      v = randn(P,1)
    end
    for j = 1:i-1
      v = v - (v'*B[:,j]) .* B[:,j]
    end
    B[:,i] = v / norm(v)
  end
  return B
end

# Use this when no adaptive encoding is used, i.e. the ACD algorithm is a normal CD.
type NoEncoding
  B
  NoEncoding(mu, P) = begin
    new(eye(P, P))
  end
end

function update!(e::NoEncoding, population::Matrix{Float64}, howOftenUpdateRotation = 1)
  # Do nothing since we never update anything...
end

type AdaptiveEncoding
  mu
  P
  B
  Bo
  invB
  weights
  mucov
  alpha_p
  c1
  cmu
  cp
  pc
  pcmu
  xmean
  C
  Cold
  diagD
  ps
  iter
  AdaptiveEncoding(mu::Int64, P::Int64, B::Matrix{Float64}, population = randn(P, mu)) = begin
    # We assume B is orthogonal!
    weights = ones(mu) / mu # Can also try non-uniform weights here...
    new(mu, P, B, B, B',
      weights,
      mu, # mucov
      1, # alpha_p
      0.5 / P, # c1
      0.5 / P, # cmu
      1 / sqrt(P), # cp
      zeros(P), # pc
      zeros(P), # pcmu
      population * weights, # xmean
      eye(P), # C
      eye(P), # Cold
      ones(P), # diagD
      0, # ps
      0  # iter
    )
  end
  AdaptiveEncoding(mu::Int64, P::Int64) = begin
    AdaptiveEncoding(mu, P, generate_random_orthogonal_transformation(P))
  end
end

function update!(ae::AdaptiveEncoding, population::Matrix{Float64}, howOftenUpdateRotation = 1)
  ae.iter += 1

  #
  # Update the adaptive ancoding, see lines 8-14 on Algorithm 2 of the Loschilov ACD paper.
  #

  xold = ae.xmean # Line 8

  # Line 9 but ensure only first mu candidates of population is used => they should be the best ones
  # in fitness-sorted order.
  if size(population, 2) != ae.mu
    pop = population[:,1:ae.mu]
  else
    pop = population
  end
  ae.xmean = pop * ae.weights  # Line 9

  # Line 10: Calculate z0 but check for the denominator being zero...
  xdiff = ae.xmean .- xold
  denom = sum((ae.invB * xdiff).^2)
  if denom == 0
    z0 = zeros(ae.P)
  else
    z0 = (sqrt(ae.P) / sqrt(denom)) * xdiff
  end

  # Line 11&13: Calc zi's and then Cmu.
  #ae.cmu = 0.0
  if ae.cmu != 0.0
    xis_minus_xold = broadcast(-, pop, xold)
    denoms = sum((ae.invB * xis_minus_xold).^2, 1)
    zis = broadcast(*, (sqrt(ae.P) ./ sqrt(denoms)), xis_minus_xold)
    Cmu = broadcast(*, ae.weights', zis * zis')
  else
    Cmu = 1 # dummy just so it is defined
  end

  # Line 12: Update p
  ae.pc = (1-ae.cp) * ae.pc .+ sqrt(ae.cp*(2-ae.cp)) * z0

  # Line 14: Update C.
  ae.C = (1-ae.c1 - ae.cmu) * ae.C .+ ae.c1 * (ae.pc * ae.pc') .+ ae.cmu * Cmu

  if ae.iter % howOftenUpdateRotation == 0 || ae.iter <= 2
    # Go back to old covar matrix if latest has NaN or Inf element
    if any(isnan(ae.C)) || any(isinf(ae.C))
      ae.C = ae.Cold
    end

    # Ensure covar matrix is symmetric by copying the upper triangular part to the lower tri part.
    ae.C = (triu(ae.C) + triu(ae.C, 1)')

    # Now do the eigen decomposition. This is the costly operation...
    eigenvalues, ae.Bo = eig(ae.C)

    # Limit the condition value of C and control for extreme eigen values.
    cond = 1e12
    if minimum(eigenvalues) <= 0
      eigenvalues[eigenvalues .< 0] = 0
      tmp = maximum(eigenvalues) / cond
      ae.C = ae.C .+ tmp * eye(ae.P)
      eigenvalues += tmp * ones(ae.P,1)
    end
    if maximum(eigenvalues) > cond*minimum(eigenvalues)
      tmp = maximum(eigenvalues)/cond - minimum(eigenvalues)
      ae.C = ae.C .+ tmp * eye(ae.P)
      eigenvalues += tmp*ones(ae.P,1)
    end
    if (minimum(eigenvalues) <= 0)
      # Should never happen since we take care of this in the limiting code above...
      throw("Negative eigenvalues => skipping update")
    end

    try
      ae.diagD = sqrt(eigenvalues)
      # Use Cholesky instead??
      ae.B = ae.Bo * diagm(ae.diagD) # This is sometimes a matrix rather than a vector. Strange!
      ae.invB = diagm(1 ./ ae.diagD) * ae.Bo'
      ae.Cold = ae.C
    catch
      # Reset to old covar matrix with a random perturbation to try to get out of the problems
      ae.C = ae.Cold .+ 0.01 * rand(ae.P, ae.P)
    end
  end
end


text:
# vary transition prob, mcts_depth, and number of nodes in network
cond1 = make_condition_dict(num_nets=500,mcts_d=[10],n_bnds=[8,35],seed=5,role=:trusted)
cond4 = make_condition_dict(num_nets=500,mcts_d=[3],n_bnds=[8,35],seed=2)

exp_dict = Dict(:name=>"mturk_supplement", :conditions=>Dict(1=>cond1, 4=> cond4),
                :xQ=>Dict(:inpts=>[:tprob,:N],:epocs=>1000,:cmp=>[4],:nn_loc=>"nn_logs"))

fldr = "logs"
for x in keys(exp_dict[:conditions])
    fname = make_useful_fname(exp_dict,x)
    exp_dict[:conditions][x][:fname] = fname
    exp_dict[:conditions][x][:fldr] = fldr
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
# Vector of NamedTuples
const RowTable{T} = Vector{T} where {T <: NamedTuple}

# RowTables can't support WeakRefStrings
unweakref(x) = x
unweakreftype(T) = T
Base.@pure unweakreftypes(::Type{T}) where {T <: Tuple} = Tuple{Any[unweakreftype(fieldtype(T, i)) for i = 1:fieldcount(T)]...}

# interface implementation
istable(::Type{<:RowTable}) = true
rowaccess(::Type{<:RowTable}) = true
# a Vector of NamedTuple iterates `Row`s itself
rows(x::RowTable) = x
schema(x::Vector{NamedTuple{names, types}}) where {names, types} = Schema(names, types)
materializer(x::RowTable) = rowtable

# struct to transform `Row`s into NamedTuples
struct NamedTupleIterator{S, T}
    x::T
end
Base.IteratorEltype(::Type{<:NamedTupleIterator{S}}) where {S} = S === nothing ? Base.EltypeUnknown() : Base.HasEltype()
Base.eltype(rows::NamedTupleIterator{Schema{names, T}}) where {names, T} = NamedTuple{names, unweakreftypes(T)}
Base.IteratorSize(::Type{NamedTupleIterator{S, T}}) where {S, T} = Base.IteratorSize(T)
Base.length(nt::NamedTupleIterator) = length(nt.x)
Base.size(nt::NamedTupleIterator) = (length(nt.x),)

function Base.iterate(rows::NamedTupleIterator{Schema{names, T}}, st=()) where {names, T}
    if @generated
        vals = Tuple(:(unweakref(getproperty(row, $(fieldtype(T, i)), $i, $(Meta.QuoteNode(names[i]))))) for i = 1:fieldcount(T))
        return quote
            x = iterate(rows.x, st...)
            x === nothing && return nothing
            row, st = x
            return $(NamedTuple{names, unweakreftypes(T)})(($(vals...),)), (st,)
        end
    else
        x = iterate(rows.x, st...)
        x === nothing && return nothing
        row, st = x
        return NamedTuple{names, unweakreftypes(T)}(Tuple(unweakref(getproperty(row, fieldtype(T, i), i, names[i])) for i = 1:fieldcount(T))), (st,)
    end
end

# unknown schema case
function Base.iterate(rows::NamedTupleIterator{Nothing, T}, st=()) where {T}
    x = iterate(rows.x, st...)
    x === nothing && return nothing
    row, st = x
    names = Tuple(propertynames(row))
    return NamedTuple{names}(Tuple(unweakref(getproperty(row, nm)) for nm in names)), (st,)
end

namedtupleiterator(::Type{T}, rows::S) where {T <: NamedTuple, S} = rows
namedtupleiterator(::Type{T}, rows::S) where {T, S} = NamedTupleIterator{typeof(schema(rows)), S}(rows)

# sink function
function rowtable(itr::T) where {T}
    istable(T) || throw(ArgumentError("Vector of NamedTuples requires a table input"))
    r = rows(itr)
    return collect(namedtupleiterator(eltype(r), r))
end

function rowtable(rt::RowTable, itr::T) where {T}
    istable(T) || throw(ArgumentError("Vector of NamedTuples requires a table input"))
    r = rows(itr)
    return append!(rt, namedtupleiterator(eltype(r), r))
end

# NamedTuple of arrays of matching dimensionality
const ColumnTable = NamedTuple{names, T} where {names, T <: NTuple{N, AbstractArray{S, D} where S}} where {N, D}
rowcount(c::ColumnTable) = length(c) == 0 ? 0 : length(c[1])

# interface implementation
istable(::Type{<:ColumnTable}) = true
columnaccess(::Type{<:ColumnTable}) = true
# a NamedTuple of AbstractVectors is itself a `Columns` object
columns(x::ColumnTable) = x
schema(x::T) where {T <: ColumnTable} = Schema(names(T), _types(T))
materializer(x::ColumnTable) = columntable

_eltype(::Type{A}) where {A <: AbstractVector{T}} where {T} = T
Base.@pure function _types(::Type{NT}) where {NT <: NamedTuple{names, T}} where {names, T <: NTuple{N, AbstractVector{S} where S}} where {N}
    return Tuple{Any[ _eltype(fieldtype(NT, i)) for i = 1:fieldcount(NT) ]...}
end

getarray(x::AbstractArray) = x
getarray(x) = collect(x)

function columntable(sch::Schema{names, types}, cols) where {names, types}
    if @generated
        vals = Tuple(:(getarray(getproperty(cols, $(fieldtype(types, i)), $i, $(Meta.QuoteNode(names[i]))))) for i = 1:fieldcount(types))
        return :(NamedTuple{names}(($(vals...),)))
    else
        return NamedTuple{names}(Tuple(getarray(getproperty(cols, fieldtype(types, i), i, names[i])) for i = 1:fieldcount(types)))
    end
end
# unknown schema case
columntable(::Nothing, cols) = NamedTuple{Tuple(propertynames(cols))}(Tuple(getarray(col) for col in eachcolumn(cols)))

function columntable(itr::T) where {T}
    istable(T) || throw(ArgumentError("NamedTuple of AbstractVectors requires a table input"))
    cols = columns(itr)
    cols isa ColumnTable && return cols
    return columntable(schema(cols), cols)
end
columntable(x::ColumnTable) = x

function ctappend(ct1::NamedTuple{N1, T1}, ct2::NamedTuple{N2, T2}) where {N1, T1, N2, T2}
    if @generated
        appends = Expr(:block, Any[:(append!(ct1[$(Meta.QuoteNode(nm))], ct2[$(Meta.QuoteNode(nm))])) for nm in N1]...)
        return quote
            $appends
            return ct1
        end
    else
        foreach(nm->append!(ct1[nm], ct2[nm]), nms)
        return ct1
    end
end

columntable(ct::ColumnTable, itr) = ctappend(ct, columntable(itr))


text:
#---------------------------------------------------------
# # [Examples](@id 1-examples)
#---------------------------------------------------------

# These examples illustrate how to use `MIRTjim`.

# First we tell Julia we are using this package:

using MIRTjim: jim, prompt

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each image is displayed.

isinteractive() && jim(:prompt, true);

# ## Simple 2D image

# The simplest example is a 2D array.
# Note that `jim` is designed to show a function `f(x,y)`
# sampled as an array `z[x,y]` so the 1st index is horizontal direction.

x, y = 1:9, 1:7
f(x,y) = x * (y-4)^2
z = f.(x, y') # 9 × 7 array
jim(z ; xlabel="x", ylabel="y", title="f(x,y) = x * (y-4)^2")

# Compare with `Plots.heatmap` to see the differences
# (transpose, color, wrong aspect ratio, distractingly many ticks):

import Plots
Plots.heatmap(z, title="heatmap")

#-
isinteractive() && prompt();

# Images often should include a title, so `title =` is optional.
jim(z, "hello")

# ## OffsetArrays

# `jim` displays the axes naturally.

using OffsetArrays
zo = OffsetArray(z, (-3,-1))
jim(zo, "OffsetArray example")

# ## 3D arrays

# `jim` automatically makes 3D arrays into a mosaic.

z3 = reshape(1:(9*7*6), (9, 7, 6))
jim(z3, "3D")

# ## Arrays of images

# `jim` automatically makes arrays of images into a mosaic.

z3 = reshape(1:(9*7*6), (9, 7, 6))
z4 = [z3[:,:,(j-1)*3+i] for i=1:3, j=1:2]
jim(z4, "Arrays of images")

# ## Units

# `jim` supports units, with axis and colorbar units appended naturally,
# thanks to UnitfulRecipes.jl.

using Unitful

x = (1:9)u"m/s"
y = (1:7)u"s"
zu = x * y'
jim(x, y, zu, "units" ;
    clim=(0,40).*u"m", xlabel="rate", ylabel="time", colorbar_title="distance")

# See `UnitfulRecipes.jl` to customize the units

# ## AxisArrays

# `jim` displays the axes (names and units) naturally by default:

using AxisArrays
using Unitful: μm, s

x = (1:9)μm
y = (1:7)μm/s
za = AxisArray(x * y'; x, y)
jim(za, "AxisArray")

# ## Options

# See the docstring for `jim` for its many options.
# Here are some defaults.

jim(:defs)

# One can set "global" defaults using appropriate keywords from above list.
# Use `:push!` and `:pop!` for such changes to be temporary.

jim(:push!) # save current defaults
jim(:colorbar, :none) # disable colorbar for subsequent figures
jim(:yflip, false) # have "y" axis increase upward
jim(rand(9,7), "rand", color=:viridis) # kwargs... passed to heatmap()

#-
jim(:pop!); # restore

# ## Layout

# One can use `jim` just like `plot` with a layout of subplots.
# The `gui=true` option is useful when you want a figure to appear
# even when other code follows.
# Often it is used with the `prompt=true` option (not shown here).

p1 = jim(rand(5,7))
p2 = jim(rand(6,8), color=:viridis)
p3 = jim(rand(9,7), title="plot 3")
jim(p1, p2, p3; layout=(1,3), gui=true)

# Just for debugging:
using Pkg
Pkg.status("MIRTjim")


text:
"""
    outlier_fraction(y)

Determine the fraction of outliers in a given vector.

Parameters
----------
    y::Labels
An array containing "normal" and "outlier" classes.

Returns
----------
    outlier_fraction::Float64
The fraction of outliers.
"""
outlier_fraction(y::Labels) = sum(y .== OutlierDetection.CLASS_OUTLIER) / length(y)

"""
    n_normal(y)
Determine the amount of normal data points in a given vector.

Parameters
----------
    y::Labels
An array containing "normal" and "outlier" classes.

Returns
----------
    n_normal::Int64
The amount of normal data points.
"""
n_normal(y::Labels) = sum(y .== OutlierDetection.CLASS_NORMAL)

"""
    n_outlier(y)
Determine the amount of outlier data points in a given vector.

Parameters
----------
    y::Labels
An array containing "normal" and "outlier" classes.

Returns
----------
    outliers::Int64
The amount of outlier data points.
"""
n_outlier(y::Labels) = sum(y .== OutlierDetection.CLASS_OUTLIER)


text:
"""
```
pseudo_measurement(m::Model990{T}, TTT::Matrix{T}, RRR::Matrix{T},
                   CCC::Vector{T}; reg::Int = 1,
                   TTTs::Vector{<: AbstractMatrix{T}} = Matrix{T}[],
                   CCCs::Vector{<: AbstractVector{T}} = Vector{T}[],
                   information_set::UnitRange = reg:reg,
                   memo::Union{ForwardMultipleExpectationsMemo, Nothing} = nothing) where {T <: AbstractFloat}
```
Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement(m::Model990{T}, TTT::Matrix{T}, RRR::Matrix{T},
                            CCC::Vector{T}; reg::Int = 1,
                            TTTs::Vector{<: AbstractMatrix{T}} = Matrix{T}[],
                            CCCs::Vector{<: AbstractVector{T}} = Vector{T}[],
                            information_set::UnitRange = reg:reg,
                            memo::Union{ForwardMultipleExpectationsMemo, Nothing} = nothing) where {T <: AbstractFloat}

    endo   = m.endogenous_states
    pseudo = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    # Initialize pseudo ZZ and DD matrices
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## Flexible Output
    ZZ_pseudo[pseudo[:y_f_t],endo[:y_f_t]] = 1.

    ## Natural Rate
    ZZ_pseudo[pseudo[:NaturalRate],endo[:r_f_t]] = 1.
    DD_pseudo[pseudo[:NaturalRate]]              = 100. * (m[:rstar]-1.)

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## Output Gap
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_t]] = 1;
    ZZ_pseudo[pseudo[:OutputGap],endo[:y_f_t]] = -1;

    ## Ex Ante Real Rate
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:R_t]]  = 1;
    ZZ_pseudo[pseudo[:ExAnteRealRate],endo[:Eπ_t]] = -1;
    DD_pseudo[pseudo[:ExAnteRealRate]]             = m[:Rstarn] - 100*(m[:π_star]-1);

    ## Long Run Inflation
    ZZ_pseudo[pseudo[:LongRunInflation],endo[:π_star_t]] = 1.
    DD_pseudo[pseudo[:LongRunInflation]]                 = 100. *(m[:π_star]-1.)

    ## Marginal Cost
    ZZ_pseudo[pseudo[:MarginalCost],endo[:mc_t]] = 1.

    ## Wages
    ZZ_pseudo[pseudo[:Wages],endo[:w_t]] = 1.

    ## Flexible Wages
    ZZ_pseudo[pseudo[:FlexibleWages],endo[:w_f_t]] = 1.

    ## Hours
    ZZ_pseudo[pseudo[:Hours],endo[:L_t]] = 1.

    ## Flexible Hours
    ZZ_pseudo[pseudo[:FlexibleHours],endo[:L_f_t]] = 1.

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Z_t
    ZZ_pseudo[pseudo[:Z_t], endo[:z_t]] = 1.
    DD_pseudo[pseudo[:Z_t]]             = 100. * exp(m[:z_star]-1)

    ## ztil_t
    ZZ_pseudo[pseudo[:ztil_t], endo[:ztil_t]] = 1.

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:Rstarn]

    ## Real FFR
    ZZ_pseudo[pseudo[:RealFFR], endo[:R_t]] = 1.
    ZZ_pseudo[pseudo[:RealFFR], endo[:π_t]] = -1.
    DD_pseudo[pseudo[:RealFFR]] =  m[:Rstarn] - (100. * (m[:π_star] - 1.))

    ## Expected nominal natural rate
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:r_f_t]] = 1.
    ZZ_pseudo[pseudo[:ExpectedNominalNaturalRate], endo[:Eπ_t]]  = 1.
    DD_pseudo[pseudo[:ExpectedNominalNaturalRate]]               =  m[:Rstarn]

    # Collect indices and transforms
    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(
        inradius,
        outradius,
        center,
        ls,
        ms,
        lc,
        mc,
        theta0,
    )
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function displacement_error(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor;
    eta = +1,
)
    L = W = width

    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCellDG.lame_coefficients(stiffness, -1)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    analyticalsolution = AnalyticalSolution(
        inradius,
        outradius,
        center,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
    levelsetbasis = HermiteTensorProductBasis(2)

    dim, nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    mesh = CutCellDG.DGMesh(
        [0.0, 0.0],
        [width, width],
        [nelmts, nelmts],
        refpoints,
    )
    cgmesh = CutCellDG.CGMesh([0.0, 0.0], [width, width], [nelmts, nelmts], nf)

    levelset = CutCellDG.LevelSet(
        x -> -circle_distance_function(x, center, inradius)[1],
        cgmesh,
        levelsetbasis,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    mergedmesh =
        CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        elasticitybasis,
        cellquads,
        stiffness,
        mergedmesh,
    )
    CutCellDG.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        cellquads,
        mergedmesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_interelement_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        facequads,
        mergedmesh,
    )

    CutCellDG.assemble_coherent_interface_condition!(
        sysmatrix,
        elasticitybasis,
        interfacequads,
        stiffness,
        mergedmesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_coherent_interface_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        interfacequads,
        mergedmesh,
    )

    CutCellDG.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        analyticalsolution,
        elasticitybasis,
        facequads,
        stiffness,
        mergedmesh,
        x -> onboundary(x, L, W),
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_transformation_linear_form!(
        sysrhs,
        transfstress,
        elasticitybasis,
        facequads,
        mergedmesh,
        x -> onboundary(x, L, W),
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mergedmesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mergedmesh)

    solution = matrix \ rhs
    nodaldisplacement = reshape(solution, 2, :)

    err = mesh_L2_error(
        nodaldisplacement,
        analyticalsolution,
        elasticitybasis,
        cellquads,
        mergedmesh,
    )

    return err
end

function test_circular_interface()
    lambda1, mu1 = 100.0, 80.0
    lambda2, mu2 = 80.0, 60.0
    theta0 = -0.067
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

    width = 1.0
    penaltyfactor = 1e3

    polyorder = 2
    numqp = required_quadrature_order(polyorder) + 2

    center = [width / 2, width / 2]
    inradius = width / 4
    outradius = width

    powers = [3, 4, 5]
    nelmts = [2^p + 1 for p in powers]

    err = [
        displacement_error(
            width,
            center,
            inradius,
            outradius,
            stiffness,
            theta0,
            ne,
            polyorder,
            numqp,
            penaltyfactor,
        ) for ne in nelmts
    ]

    dx = 1.0 ./ nelmts
    u1err = [er[1] for er in err]
    u2err = [er[2] for er in err]

    u1rate = convergence_rate(dx, u1err)
    u2rate = convergence_rate(dx, u2err)

    @test all(u1rate .> 2.8)
    @test all(u2rate .> 2.8)
end

function test_corner_curved_interface()
    lambda1, mu1 = 100.0, 80.0
    lambda2, mu2 = 80.0, 60.0
    theta0 = -0.067
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

    width = 1.0
    penaltyfactor = 1e3

    polyorder = 2
    numqp = required_quadrature_order(polyorder) + 2

    center = [width, width]
    inradius = width / 2
    outradius = 2width
    powers = [3, 4, 5]
    nelmts = [2^p + 1 for p in powers]

    err = [
        displacement_error(
            width,
            center,
            inradius,
            outradius,
            stiffness,
            theta0,
            ne,
            polyorder,
            numqp,
            penaltyfactor,
        ) for ne in nelmts
    ]

    dx = 1.0 ./ nelmts
    u1err = [er[1] for er in err]
    u2err = [er[2] for er in err]

    u1rate = convergence_rate(dx, u1err)
    u2rate = convergence_rate(dx, u2err)

    @test all(u1rate .> 2.9)
    @test all(u2rate .> 2.9)
end

test_circular_interface()
test_corner_curved_interface()


text:
using Ogg
using Base.Test

testdir = dirname(@__FILE__)

# Let's start with building our own Ogg structure, writing it out to an IOBuf,
# then loading it back in again and checking everything about it we can think of

# We are going to build three streams, each with 10 packets
num_packets = 10
stream_ids = rand(Cint, 3)
packets = Dict{Clong,Vector{Vector{UInt8}}}()
granulepos = Dict{Int64,Vector{Int64}}()
for serial in stream_ids
    # The packets are all of different size
    packets[serial] = Vector{UInt8}[rand(UInt8, 100*x) for x in 1:num_packets]

    # Each packet will have a monotonically increasing granulepos, except for
    # the first two packets which are our "header" packets with granulepos == 0
    granulepos[serial] = Int64[0, 0, [20*x for x in 1:(num_packets - 2)]...]
end

# Now we write these packets out to an IOBuffer
ogg_file = IOBuffer()
save(ogg_file, packets, granulepos)

# Rewind to the beginning of this IOBuffer and load the packets back in
seekstart(ogg_file)

# This is the verbose way to load an Ogg file that leaves the decoder object
# around for us to inspect; the load() function returns only `dec.packets`
dec = OggDecoder()
Ogg.decode_all_pages(dec, ogg_file)
Ogg.decode_all_packets(dec, ogg_file)

# # Are the serial numbers the same?
@test sort(collect(keys(dec.packets))) == sort(stream_ids)

# Are the number of packets the same?
for serial in stream_ids
    @test length(dec.packets[serial]) == length(packets[serial])
end

# Are the contents of the packets the same?
for serial in stream_ids
    for packet_idx in 1:length(dec.packets[serial])
        @test dec.packets[serial][packet_idx] == packets[serial][packet_idx]
    end
end

# Let's dig deeper; let's ensure that the first two pages had length equal to
# our first two packets, proving that our header packets had their own pages:
for serial in stream_ids
    @test dec.pages[serial][1].body_len == length(packets[serial][1])
    @test dec.pages[serial][2].body_len == length(packets[serial][2])
end

# Next, let's load a known ogg stream and ensure that it's exactly as we expect
ogg_packets = load(joinpath(testdir, "zero.ogg"))

# There is only one stream, and we know its serial number
@test collect(keys(ogg_packets)) == [1238561138]
serial = first(keys(ogg_packets))

# There are four packets, the first starts with \x7fFLAC
@test length(ogg_packets[serial]) == 4
@test bytestring(ogg_packets[serial][1][2:5]) == "FLAC"

# The lengths of the packets are:
@test [length(x) for x in ogg_packets[serial]] == [51, 55, 13, 0]


text:
module PoincareInvariants

    using HDF5
    using FFTW
    using OffsetArrays
    using ProgressMeter
    using Requires

    using GeometricIntegrators
    using GeometricIntegrators.Utils

    import DoubleFloats: Double64
    import LinearAlgebra: ⋅

    include("matrix_utils.jl")

    export PoincareInvariant1st, PoincareInvariant1stCanonical,
           PoincareInvariant2nd, PoincareInvariant2ndCanonical,
           PoincareInvariant2ndApproxFun,
           PoincareInvariant2ndApproxFunCanonical,
           PoincareInvariant2ndTrapezoidal,
           evaluate_poincare_invariant,
           evaluate_poincare_invariant_correction,
           write_to_hdf5

    abstract type AbstractPoincareInvariant2nd{DT} end
    abstract type AbstractPoincareInvariant2ndApproxFun{DT,ND,NC,NV} <: AbstractPoincareInvariant2nd{DT} end

    include("poincare_invariant_1st_common.jl")
    include("poincare_invariant_1st.jl")
    include("poincare_invariant_1st_canonical.jl")
    include("poincare_invariant_2nd_approxfun.jl")
    include("poincare_invariant_2nd_trapezoidal.jl")

    # function __init__()
    #     @require ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" begin
    #         @require FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c" begin
                export PoincareInvariant2ndCOP
                include("poincare_invariant_2nd_cop.jl")
    #         end
    #     end
    # end

    const PoincareInvariant2nd = PoincareInvariant2ndApproxFun
    const PoincareInvariant2ndCanonical = PoincareInvariant2ndApproxFunCanonical

end


text:
using Test
using PICDataStructures
using LookingGlass: module_recursive_globals_names

@test isempty(detect_ambiguities(PICDataStructures))

@test isempty(detect_unbound_args(PICDataStructures))

@test all(isempty.(values(module_recursive_globals_names(
    PICDataStructures,
    constness=:nonconst,
    mutability=:all)
)))


text:
"""
    load_config(yamlfile)

Load configurations which directs creation of a polymer system and its components from a YAML file given at `yamlfile`.

The configurations are stored in a nested Dicts.
"""
load_config(yamlfile) = YAML.load_file(yamlfile)

"""
    ObjType{name}

A type representing the type of an object defined in Polymer.jl.

A `ObjType` object can be passed to the first argument of [`make`](@ref) to dispatch to the correction version of this method to create an object designated by `name`:

* :System: PolymerSystem
* :Component: Component
* :χN_map: Dict{Set{Symbol},AbstractFloat}
* :Specie: AbstractSpecie (KuhnSegment) and SmallMolecule
* :SMOL: SmallMolecule
* :BCP: BlockCopolymer
* :Block: PolymerBlock
* :BlockEnd: BlockEnd (FreeEnd and BranchPoint)
"""
struct ObjType{name} end

function make(::ObjType{:BlockEnd}, config)
    if isempty(config)
        return FreeEnd(), FreeEnd()
    end
    if length(config) == 1
        return FreeEnd(), BranchPoint(Symbol(config[1]))
    end
    e1, e2 = Symbol.(config)
    return BranchPoint(e1), BranchPoint(e2)
end

function make(::ObjType{:Block}, config, sps)
    label = Symbol(config["label"])
    if haskey(config, "segment") && !isnothing(config["segment"])
        slabel = Symbol(config["segment"])
    else
        slabel= label
    end
    segment = nothing
    for sp in sps
        if sp.label == slabel
            segment = sp
            continue
        end
    end
    f = config["length"]
    E1, E2 = make(ObjType{:BlockEnd}(), config["ends"])
    return PolymerBlock(label, segment, f, E1, E2)
end

function make(::ObjType{:BCP}, config, sps)
    label = Symbol(config["label"])
    blocks = [make(ObjType{:Block}(), c, sps) for c in config["chain"]]
    return BlockCopolymer(label, blocks)
end

function make(::ObjType{:SMOL}, config, sps)
    label = Symbol(config["label"])
    for sp in sps
        if sp.label == label
            return sp
        end
    end
end

function make(::ObjType{:Component}, config, sps)
    molecule = make(ObjType{Symbol(config["type"])}(), config, sps)
    kwargs = Dict()
    if haskey(config, "length")
        kwargs[:α] = config["length"]
    end
    if haskey(config, "volume_fraction")
        kwargs[:ϕ] = config["volume_fraction"]
    end
    return Component(molecule; kwargs...)
end

function make(::ObjType{:Specie}, config)
    label = Symbol(config["label"])
    kwargs = Dict()
    if haskey(config, "b") && !isnothing(config["b"])
        kwargs[:b] = config["b"]
    end
    if haskey(config, "M") && !isnothing(config["M"])
        kwargs[:M] = config["M"]
    end
    if haskey(config, "type")
        t = config["type"]
    else
        t = "Segment"
    end
    if t == "Segment"
        return KuhnSegment(label; kwargs...)
    else
        return SmallMolecule(label; kwargs...)
    end
end

function make(::ObjType{:χN_map}, config)
    sp1, sp2, χN = first(config)
    χN_map = Dict{Set{Symbol}, eltype(χN)}()
    for a in config
        s1, s2, χN = a
        χN_map[Set([Symbol(s1), Symbol(s2)])] = χN
    end
    return χN_map
end

"""
    make(::ObjType{:System}, config)

Make an `PolymerSystem` object.

Note: `confinement` is not implmented.
"""
function make(::ObjType{:System}, config)
    sps = [make(ObjType{:Specie}(), c) for c in config["Species"]]
    components = [make(ObjType{:Component}(), c, sps) for c in config["Components"]]
    kwargs = Dict()
    if haskey(config, "χN_map")
        χN_map = make(ObjType{:χN_map}(), config["χN_map"])
    else
        error("Must provide an interaction map.")
    end
    if haskey(config, "chain_density")
        kwargs[:C] = config["chain_density"]
    end
    return PolymerSystem(components, χN_map; kwargs...)
end

make(config) = make(ObjType{:System}(), config["System"])


text:
## input

# V: provide path for VCs files (.csv file)
# y: provide file name and column index (.csv file, with column names)
# X: provide file name and column index (.csv file)
# λgrid: provide value or use the default (vector)
# print: default is false (Bool)
# λfactor: provide vector which the unpenalized VCs set to zero other wise the
#          default is vector of 1 with length equal to number of VCs (vector)

# criterion: None, AIC, BIC, default is AIC (String)
#            if None, will output estimations under each tuning parameter

## output

function microVClasso(args...; resFile::String= "", covFile::String = "",
  Vpath::String = "", outpath::String = "",
  covIdx ::Array{Int, 1} = [3,4], resIdx ::Int = 3, idIdx::Int = 1,
  yhead::Bool = true, xhead::Bool = true, longitudinal::Bool = false,
  λgrid::Array{Float64, 1} = Float64[], print::Bool = false,
  λfactor::Array{Float64, 1} = Float64[], criterion::String = "AIC")

  # input response data: y
    if isempty(resFile)
      error("microVClasso:ymissing\n",
      "# need to input responses file name");
    else
      phenotye = readcsv(resFile)
      if yhead
        y = phenotye[2:end,resIdx];
        nper = phenotye[2:end,idIdx]
      else
        y = phenotye[:,resIdx];
        nper = phenotye[:,idIdx]
      end
      y = convert(Array{Float64, 1}, y);
    end

    n = length(y)

  # input covariates data: X

    if isempty(covFile)
      # no covariates provided
      X = ones(nPer, 1);
    else
      X = readcsv(covFile)
      if xhead
        X = X[2:end,covIdx];
      else
        X = X[:,covIdx];
      end
      X = convert(Array{Float64, 2}, X);
    end

    # input microbiome kernel matrices (variance components)

    if isempty(Vpath)
      error("microVClasso:vmisiing\n",
      "# need to input a path for kernel matrices");
    else
      vclist = readdir(Vpath)
      vcname = [split(vclist[i],".csv")[1] for i =1:length(vclist)]
      vc = []

      if longitudinal & isempty(λfactor)

        dic = countmap(nper)
        id = unique(nper)

        ZZ = cat([1,2],ones(dic[id[1]],dic[id[1]]))
        for i = 2:length(id)
          ZZ = cat([1,2],ZZ,ones(dic[id[i]],dic[id[i]]))
        end
        λfactor = vcat(0,ones(Int,length(vclist)))
        push!(vc,Hermitian(ZZ))
        vcname = vcat("ZZ",vcname)
      end

      for i=1:length(vclist)
        if endswith(Vpath,"/")
          psd = readcsv(string(Vpath,vclist[i]))
        else
          psd = readcsv(string(Vpath,"/",vclist[i]))
        end
        push!(vc, Hermitian(psd[:,:]))
      end
      push!(vc, eye(n,n)) # add the error term
    end

    vc = collect(vc)
    vc = convert(Array{Array{Float64,2},1}, vc)

    if isempty(λgrid)
      λgrid = vcat(0.1:0.1:0.6,1.35,5,10)
    end

    if longitudinal
      Lasso = VClasso(y, X, vc, λfactor,false, λgrid)
    else
      λfactor = ones(Int,length(vc)-1)
      Lasso = VClasso(y, X, vc, λfactor,false, λgrid)
    end

    σ2 = Lasso[1][1:(end-1),:]
    nvc = [countnz(σ2[:,col]) for col=1:size(σ2,2)]
    k = nvc .+ size(X, 2)
    nloglik = Lasso[3]

    AIC = 2 * k - nloglik * 2
    BIC = log(n) * k - nloglik * 2

    σ2AIC = σ2[:,indmin(AIC)]
    vcAIC = vcname[σ2AIC.!=0] # output selected name and the estimates
    outAIC = hcat(vcAIC,σ2AIC[σ2AIC.!=0])

    σ2BIC = σ2[:,indmin(BIC)]
    vcBIC = vcname[σ2BIC.!=0] # output selected name and the estimates
    outBIC = hcat(vcBIC,σ2BIC[σ2BIC.!=0])

    if endswith(outpath,"/")
      outpath = outpath
    else
      outpath = string(outpath,"/")
    end

    writedlm(string(outpath,"vcest_all.txt"),σ2)
    writedlm(string(outpath,"loglik_all.txt"),hcat(λgrid,nloglik))

    if criterion == "None"
      return σ2, nloglik
    elseif criterion == "AIC"
      writedlm(string(outpath,"AIC_all.txt"),hcat(λgrid,AIC))
      writedlm(string(outpath,"vc_AIC_",λgrid[indmin(AIC)],".txt"),outAIC)
      return σ2, nloglik, AIC, outAIC
    elseif criterion == "BIC"
      writedlm(string(outpath,"BIC_all.txt"),hcat(λgrid,BIC))
      writedlm(string(outpath,"vc_BIC_",λgrid[indmin(BIC)],".txt"),outBIC)
      return σ2, nloglik, BIC, outBIC
    elseif criterion == "Both"
      writedlm(string(outpath,"AIC_BIC_all.txt"),hcat(λgrid,AIC,BIC))
      writedlm(string(outpath,"vc_AIC_",λgrid[indmin(AIC)],".txt"),outAIC)
      writedlm(string(outpath,"vc_BIC_",λgrid[indmin(BIC)],".txt"),outBIC)
      return σ2, nloglik, AIC, outAIC, BIC, outBIC
    end

  end


text:
using StemCellModels
using DataFrames
using Pkg
Pkg.build("RCall")
using RCall
using ProgressMeter
using StatsBase
using Statistics
using Random
using StatsBase
using ArgParse
R"""
library(ggplot2)
library(cowplot)
library(tidyverse)
library(viridis)
library(jcolors)
""";

# This script contains some functions to do a Maximum Likelihood fit of our model
include("optim.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--examplefitsout"
        help = "File to save example fits"
    "--exampledifferentbins"
        help = "File to save example fits"
    "--nsamples"
        help = "Number of samples to take"
        arg_type = Int
        default = 10
end

parsed_args = parse_args(ARGS, s)

Random.seed!(12345)

#function to calculate interval dN/dS from simulations
function intervaldnds_data(Msel, Mneut, mmin, Nmax; mgap = 5.0, mup = 1.0, mud = 1.0)
    #calculate cumulative dn/ds across cohort, normalizing for mutation rates
    dndscumsum = Float64[]
    dnvec = Float64[]
    dsvec = Float64[]
    driverVAFt = filter(x -> x >= mmin, Msel)
    passengerVAFt = filter(x -> x >= mmin, Mneut)
    freqs = collect(mmin:mgap:Nmax)
    for f in freqs
        dn = length(filter(x -> x< f, driverVAFt))/mud
        ds = length(filter(x -> x< f, passengerVAFt))/mup
        push!(dndscumsum, dn/ds)
        push!(dnvec, dn)
        push!(dsvec, ds)
    end
    return dndscumsum[2:end], freqs[2:end], dnvec[2:end], dsvec[2:end]
end

#function to simulate population of cells with different paramters
function simulatepopulation(;Δ = 0.0, ρ = 100.0, Amin = 1 ./ ρ, gap = 0.02, tend = 20.0,
        Ncells = 1000, r = 0.5, λ = 1.0,
        sample = false, nsample = 10, Amax = 0.0)

    rangeN = 1:10^4
    Nits = 100
    Cntd = zeros(length(rangeN), Nits)
    Cntp = zeros(length(rangeN), Nits)
    Msel = Int64[]
    Mneut = Int64[]
    SMtog = SkinStemCellModel(Ncells, Δmut = Δ, μp = 0.001, μd = 0.001, tend = tend, r = r, λ = λ)
    for i in 1:Nits
        scst = runsimulation(SMtog, progress = false, finish = "time", onedriver = true, restart = true)
        append!(Msel, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[2]), rangeN)))
        append!(Mneut, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[1]), rangeN)))
    end

    Asel = Msel ./ ρ
    Aneut = Mneut ./ ρ
    Asel = filter(x -> x >= Amin, Asel)
    Aneut = filter(x -> x >= Amin, Aneut)

    if sample
        prevlength = length(Asel)
        Asel = Asel[StatsBase.sample(1:length(Asel), nsample, replace = false)]
        Aneut = Aneut[StatsBase.sample(1:length(Aneut),
                Int64(round(nsample * (length(Aneut) / prevlength))), replace = false)]
    end

    if Amax == 0.0
        Amax = maximum(Asel)
    end

    x1 = intervaldnds_data(Asel, Aneut, Amin, Amax; mgap = gap, mup = SMtog.μp, mud = SMtog.μd)

    DFinterval = DataFrame(dnds = x1[1], A = x1[2], dn = x1[3] * SMtog.μd / Nits,
    ds = x1[4] * SMtog.μd / Nits)

    return DFinterval, SMtog, [DFinterval[:dn][end] * Nits, DFinterval[:ds][end] * Nits]
end

#####################################################
# Power to recover paramters
#####################################################
println("Running simulation to test power to recover paramters")

#create some empty vectors and dataframes to store data
deltavec = Float64[]
rlambdavec = Float64[]
nmutsdn = Float64[]
nmutsds = Float64[]
itvec = Int64[]
sampvec = Float64[]
rsq = Float64[]
myDF = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Int64, Int64],
[:dnds, :A, :dndsfit, :deltafit, :lambdarfit, :rsq, :deltatrue, :lambdartrue, :its, :nsamples], 0)

myDF = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Int64, Int64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda, :rsq, :deltatrue, :lambdartrue,
    :its, :nsamples], 0)

#number od mutations to be samples
for nsamples in [5, 8, 10, 25, 50, 100, 250]
    #repeat 100 times
    for i in 1:parsed_args["nsamples"]
        Random.seed!(i)
        Δ = 0.25
        r = 0.5
        lambda = 1.0
        rlambda = r * lambda
        DF1, SM, nmuts = simulatepopulation(;Δ = Δ, tend = 30.0, Amin = 0.05, ρ = 100.0,
                            λ = lambda, r = r, sample = true, nsample = nsamples + 1)
        println(nmuts)
        if isnan(DF1[:dnds][1]) | (DF1[:dnds][1] == Inf)
            continue
        end
        x = LLoptimizationresults(DF1[:dnds], DF1[:A]; t = 30.0, Amin = 0.05, ρ = 100.0)
        x.DF[:deltatrue] = Δ
        x.DF[:lambdartrue] = rlambda
        x.DF[:its] = i
        x.DF[:nsamples] = nsamples
        append!(myDF, x.DF)
        push!(deltavec, x.Δ)
        push!(rlambdavec, x.rλ)
        push!(itvec, i)
        push!(nmutsdn, nmuts[1])
        push!(nmutsds, nmuts[2])
        push!(sampvec, nsamples)
        push!(rsq, x.DF[:rsq][1])
    end
end

DFfit = DataFrame(iterations = itvec,
nmutsdn = nmutsdn,
nmutsds = nmutsds,
delta = deltavec,
rlambda = rlambdavec,
nsamples = sampvec,
rsq = rsq)

DFfit[:product] = DFfit[:delta] .* DFfit[:rlambda];

println("Write data to file")
@rput DFfit
R"""
write_csv(DFfit, $(parsed_args["powerout"]))
""";


text:
int main() {
     int i = foo(true);
     return 0 ;
}

int foo(boolean b) { b = true; }

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
export topological_sort

# Returns a topolically sorted ordering of the states in the supplied wfst.
function topological_sort(wfst::Wfst)
    graph = deepcopy(wfst)
    l = Any[]
    s = states_with_no_in_arcs(wfst)
    while length(s) > 0
        n = pop!(s)
        push!(l, n)
        debug(string("s: ", s))
        debug(string("l: ", l))
        # Remove arcs that lead from n to other nodes
        for arc in graph.arcs
            if arc.from == n
                debug(string("arc to delete: ", arc))
                graph.arcs = setdiff(graph.arcs, Set([arc]))
                debug(string("graph.arcs ", graph.arcs))
                # Add to the queue those other states have other incoming arcs.
                if arc.to in states_with_no_in_arcs(graph)
                    s = union(s, Set([arc.to]))
                end
            end
        end
    end
    debug(string("graph.arcs: ", graph.arcs))
    if length(graph.arcs) > 0
        error("Graph is cyclic.")
    else
        return l
    end
end

# Returns a set of states that have no arcs coming in.
function states_with_no_in_arcs(wfst::Wfst)
    states_with_in_arcs = Set([arc.to for arc in wfst.arcs])
    return setdiff(wfst.states, states_with_in_arcs)
end


text:
"""
    CoulombReactionField(; cutoff_dist, matrix_14, solvent_dielectric, nl_only,
                            weight_14, coulomb_const, force_unit, energy_unit)

The Coulomb electrostatic interaction modified using the reaction field approximation.
"""
struct CoulombReactionField{D, S, W, T, F, E, D2, K, R} <: GeneralInteraction
    cutoff_dist::D
    matrix_14::BitArray{2}
    solvent_dielectric::S
    nl_only::Bool
    weight_14::W
    coulomb_const::T
    force_unit::F
    energy_unit::E
    sqdist_cutoff::D2
    krf::K
    crf::R
    krf_14::K
    crf_14::R
end

function CoulombReactionField(;
                    cutoff_dist,
                    matrix_14,
                    solvent_dielectric=78.3,
                    nl_only=false,
                    weight_14=1.0,
                    coulomb_const=coulombconst,
                    force_unit=u"kJ * mol^-1 * nm^-1",
                    energy_unit=u"kJ * mol^-1")
    sqdist_cutoff = cutoff_dist ^ 2
    krf = (1 / (cutoff_dist ^ 3)) * ((solvent_dielectric - 1) / (2 * solvent_dielectric + 1))
    crf = (1 /  cutoff_dist     ) * ((3 * solvent_dielectric) / (2 * solvent_dielectric + 1))
    krf_14 = (1 / (cutoff_dist ^ 3)) * 0
    crf_14 = (1 /  cutoff_dist     ) * 0
    return CoulombReactionField{typeof(cutoff_dist), typeof(solvent_dielectric), typeof(weight_14),
                                typeof(coulomb_const), typeof(force_unit), typeof(energy_unit),
                                typeof(sqdist_cutoff), typeof(krf), typeof(crf)}(
        cutoff_dist, matrix_14, solvent_dielectric, nl_only, weight_14, coulomb_const,
        force_unit, energy_unit, sqdist_cutoff, krf, crf, krf_14, crf_14)
end

@inline @inbounds function force(inter::CoulombReactionField,
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size)
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    if r2 > inter.sqdist_cutoff
        return ustrip.(zero(coord_i)) * inter.force_unit
    end

    coulomb_const = inter.coulomb_const
    qi, qj = atom_i.charge, atom_j.charge
    r = √r2
    i, j = atom_i.index, atom_j.index
    if inter.matrix_14[i, j]
        # 1-4 interactions do not use the reaction field approximation
        krf = inter.krf_14
    else
        krf = inter.krf
    end

    f = (coulomb_const * qi * qj) * (inv(r) - 2 * krf * r2) * inv(r2)

    return f * dr
end

@inline @inbounds function potential_energy(inter::CoulombReactionField,
                                    s::Simulation,
                                    i::Integer,
                                    j::Integer)
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)

    if r2 > inter.sqdist_cutoff
        return ustrip(zero(s.timestep)) * inter.energy_unit
    end

    coulomb_const = inter.coulomb_const
    qi, qj = s.atoms[i].charge, s.atoms[j].charge
    r = √r2
    if inter.matrix_14[i, j]
        # 1-4 interactions do not use the reaction field approximation
        krf = inter.krf_14
        crf = inter.crf_14
    else
        krf = inter.krf
        crf = inter.crf
    end

    return (coulomb_const * qi * qj) * (inv(r) + krf * r2 - crf)
end


text:
using Colors
using Statistics
using LinearAlgebra
using ImageSegmentation
using ColorTypes, ColorVectorSpace
export ECW,
    FRCRGBD,
    Zeboudj,
    ValuesEntropy,
    LiuYangF,
    FPrime,
    ErdemMethod,
    Q

import Base:
    -,
    iterate

(-)(a::Lab{T1}, b::Lab{T2}) where T1 where T2 =
    Lab{T1}(comp1(a)-comp1(b),
            comp2(a)-comp2(b),
            comp3(a)-comp3(b))

iterate(a::Lab) = (comp1(a), 1)
function iterate(a::Lab, i::Integer)
    if i == 1
        return (comp2(a), i+1)
    elseif i == 2
        return (comp3(a), i+1)
    else
        return nothing
    end
end
function evaluate(algo, img, segmented_image::SegmentedImage)
    return evaluate(algo, img, labels_map(segmented_image))
end

# On Selecting the Best Unsupervised Evaluation Techniques for Image Segmentation
"""
The use of visible color difference in the quantitative evaluation of color image segmentation,
Hsin-Chia Chen and Sheng-Jyh Wang
"""
struct ECW
    threshold::Float64
    ECW(;threshold::Float64) = new(threshold)
end
function evaluate(c::ECW, image::AbstractArray{Color, N}, segments::SegmentedImage) where T<:Integer where {Color<:Lab, N}
    segments_mean = Dict(i=>convert(Lab, segment_mean(segments,i))
                     for i in segments.segment_labels)
    R = length(segments_mean)
    mean_segments = collect(map(
            i->segments_mean[i],
            segments.image_indexmap))
    E_intra = sum(norm.(image-mean_segments).<  c.threshold)
    C = 1/6
    total_sum = 0
    bool_map = [segments.image_indexmap .== i for i=1:R]
    for i=1:R-1
        a = bool_map[i]
        for j=i+1:R
            b = bool_map[j]
            if (norm(segments_mean[i] - segments_mean[j])  > c.threshold)
                 Kij = sum((a[1:end-1,1:end] + b[2:end,1:end]).==2)  + sum((a[1:end,1:end-1] + b[1:end,2:end]).==2)
                total_sum+=Kij
            end
        end
    end
    E_inter= (2*total_sum)/(C*prod(size(segments.image_indexmap)))
    return 0.5*(E_intra +  E_inter)
end

function evaluate(c::ECW, image::AbstractArray{Color, N}, segments::SegmentedImage) where {Color<:RGB, N}
    return evaluate(c, convert.(Lab, image), segments)
end

"""
Zéboudj, Rachid. Filtrage, seuillage automatique, contraste et contours: du pré-traitement à l'analyse d'image. Diss. Saint-Etienne, 1988.
Unsupervised Evaluation of Image Segmentation Application to Multi-spectral Images

"This contrast takes into account the internal and external contrast of the regions measured in the neighborhood of each pixel"
"""
struct Zeboudj
    radius::Integer
    Zeboudj(;radius::Integer=3) = new(radius)
end

function evaluate(c::Zeboudj, image::AbstractArray{Color, M}, segments::SegmentedImage) where {Color<:Colorant, M}
    image = convert.(Gray, image)
    N = maximum(segments.segment_labels)
    inside = zeros(N)
    outside = zeros(N)
    segment_sizes = zeros(N)
    border_lengths = zeros(N)
    R = CartesianIndices(segments.image_indexmap)
    I1, Iend = first(R), last(R)
    W = CartesianIndex(c.radius,c.radius)
    for I in CartesianIndices(size(image))
        current_label = segments.image_indexmap[I]
        segment_sizes[current_label] += 1
        max_inside = -1
        max_outside = -1
        is_a_border_point = false
        for J in max(I1, I-W):min(Iend, I+W)
            if current_label != segments.image_indexmap[J]
                is_a_border_point = true
                max_outside = max(max_outside, abs(image[I] - image[J]))
            else
                max_inside = max(max_inside, abs(image[I] - image[J]))
            end
        end
        if (is_a_border_point)
            border_lengths[current_label]+=1
        end
        if max_inside == -1
            max_inside = 0
        end
        if max_outside == -1
            max_outside = 0
        end
        inside[current_label] += max_inside
        outside[current_label] += max_outside;
    end
    inside=inside ./ segment_sizes
    outside=outside ./ border_lengths

    C = zeros(N)
    for j=1:N
        if (inside[j] > 0 ) &&  (inside[j] < outside[j])
            C[j] = 1 - (inside[j]/outside[j])
        elseif (inside[j]==0)
            C[j] = outside[j]
        else
            C[j] = 0
        end
    end
    return (1/(prod(size(image)))) *sum(segment_sizes.*C)
end
@doc raw"""
An Entropy-based Objective Evaluation Method for Image Segmentation
Hui Zhang*, Jason E. Fritts and Sally A. Goldman

Given an image I of ($n \times m$) , $S_I = nm$

$H_v(R_j) = - \sum\limits_{m \in V^{(v)}_j} \dfrac{L_j(m)}{S_j) log(\dfrac{L_j(m)}{S_j)
)$

$ H_l(I) = - \sum\limits_{j=1}^N \dfrac{S_j}{S_I} log(\dfrac{S_j}{S_I}) $
"""
struct ValuesEntropy
end
function components(a::Vector{T}) where T<:Color
    return (comp1.(a), comp2.(a), comp3.(a))
end
function components(a::Vector{T}) where T<:TransparentColor
    return (comp1.(a), comp2.(a), comp3.(a), alpha.(a))
end
function components(a::Vector{T}) where T<:Number
    return a
end
function evaluate(c::ValuesEntropy, image::AbstractArray{C, N}, segments::SegmentedImage) where {C<:Colorant, N}
    Hr = 0
    Hl = 0
    Si = prod(size(image))
    for j in segments.segment_labels
        segment = segments.image_indexmap .== j
        Sj = segments.segment_pixel_count[j]
        comps = components(image[segment])
        Hl += (Sj/Si)*(log(Sj/Si))
        H = 0
        for k=1:length(comps)
            for val in unique(comps[k])
                Lj = count(comps[k].==val)
                H += (Lj / Sj)*(log(Lj / Sj))
            end
        end
        Hr += (Sj/Si)*(-H)
    end
    Hl = -Hl
    return Hr + Hl
end
@doc raw"""
Multiresolution Color Image Segmentation
Jianqing Liu and Yee-Hong Yang, Senior Member, IEEE

$F(I) = \sqrt{R} \times \sum\limits_{i=1}^R \dfrac{e_i^2}{\sqrt{A_i}}$
"""
struct LiuYangF
end

function scale_factor(c::LiuYangF, image::AbstractArray{C,N}, segments::SegmentedImage) where {C<:Colorant, N}
    return sqrt(maximum(segments.segment_labels))
end
@doc raw"""
´´´julia
function color_error_sum(image::Matrix, segments::Matrix{Integer)
´´´
Computes
$\sum\limits_{i=1}^R \dfrac{e_i^2}{\sqrt{A_i}}$
where $R$ is the number of regions in the segmented image, $A_i$ is the area, or the number of pixels of the ith region $i$, and $e_i$  the color error of region $i$. e is defined as the sum of the Euclidean distance of the color vectors between the original image and the segmented image of each pixel in the region.

"""
function color_error_sum(c, image::AbstractArray{C, N}, segments::SegmentedImage) where {C<:Colorant, N}
    labels = segments.segment_labels
    out = 0
    for i in labels
        pixels_in_segment = segments.image_indexmap .== i
        A_i = segments.segment_pixel_count[i]
        e_i = sum((norm.(mean(image[pixels_in_segment]) - image[pixels_in_segment])).^2)
        out += e_i/sqrt(A_i)
    end
    return (1/(prod(size(image))))*sqrt(scale_factor(c, image, segments))*out
end
function evaluate(c::LiuYangF, image::AbstractArray{C, N}, segments::SegmentedImage) where {C<:Colorant, N}
    return scale_factor(c, image, segments)*color_error_sum(c, image, segments)
end

"""
Quantitative evaluation of color image segmentation results
M. Borsotti a, P. Campadelli a,2, R. Schettini b,
"""
struct FPrime
end

function scale_factor(c::FPrime, image::AbstractArray{C, N}, segments::SegmentedImage) where {C<:Colorant, N}
    segment_sizes = values(segments.segment_pixel_count)
    scale = 0
    for A in unique(segment_sizes)
        R_A = count(segment_sizes .== A)
        scale += R_A^(1 + (1/A))
    end
    return (1/prod(size(image))) *scale
end
function evaluate(c::FPrime, image::AbstractArray{C, N}, segments::SegmentedImage)  where {C<:Colorant, N}
    return scale_factor(c, image, segments)*color_error_sum(c, image, segments)
end

@doc raw"""
Quantitative evaluation of color image segmentation results
M. Borsotti a, P. Campadelli a,2, R. Schettini b,

$Q(I) = \dfrac{1}{10000(N \times M)} \sqrt{R} \times \sum\limits_{i=1}^R \left[  \dfrac{e_i^2}{1+\log A_i} + \left(  \dfrac{R(A_i)}{A_i} \right)^2 \right]$

where $R$ is the number of regions in the segmented image, $A_i$ is the area, or the number of pixels of the ith region $i$, and $e_i$  the color error of region $i$. e is defined as the sum of the Euclidean distance of the color vectors between the original image and the segmented image of each pixel in the region, while $R(A_i)$ represents the number of regions having an area equal to $A_i$.
"""
struct Q
end
function evaluate(c::Q, image::AbstractArray{C, N}, segments::SegmentedImage) where {C<:Colorant, N}
    segment_sizes = collect(values(segments.segment_pixel_count))
    out = 0
    for i in segments.segment_labels
        pixels_in_segment = segments.image_indexmap .== i
        A_i = segments.segment_pixel_count[i]
        e_i = sum((norm.(mean(image[pixels_in_segment]) - image[pixels_in_segment])).^2)
        R_A_i = count(segment_sizes.==A_i)
        out += e_i^2/(1 + log(A_i)) + (R_A_i/A_i)^2
    end
    return (1/(10000*prod(size(image))))*sqrt(maximum(segments.segment_labels))*out
end

"""
Performance Measures for Video Object Segmentation and Tracking
Erdem, Sankur, Tekalp
"""
struct ErdemMethod
    M::Integer # Neighborhood size to evaluate
    L::Integer # Length of the normal line
end
function mean_value(img::Matrix{T}, pos, w) where T
    R = CartesianIndices(size(img))
    I, L = first(R), last(R)
    sum = zeros(T)
    range = CartesianIndices(max(I, pos-w), min(L, pos+w))
    for N in  range
        sum += img[N]
    end
    return sum / prod(size(range))
end

function evaluate(c::ErdemMethod, image::AbstractArray{C, N}, segments::SegmentedImage) where {C<:Colorant, N}
    (h, w) = size(segments.image_indexmap)
    sum = 0
    (inside, outside) = normal_lines_extremes(boundary, c.L)
    for i=1:length(inside)
        val_inside = mean_value(image, inside[i,:], c.M)
        val_outside = mean_value(iamge, outside[i,:], c.M)
        delta_color = norm(val_outside - val_inside)/(sqrt(3*255^2))
        sum += delta_color
    end
    return 1 - (sum/length(inside))
end

# RGB-D Images
struct FRCRGBD
end
function color_std(colors)
    return std(channelview(colors))
end

function evaluate(c::FRCRGBD,
                  color_image::AbstractArray{T},
                  range_image::Matrix,
                  seg_image::SegmentedImage) where T<:Colorant
    return evaluate(c,
                    convert.(Lab, color_image),
                    range_image,
                    labels_map(seg_image))
end
function evaluate(c::FRCRGBD,
                  color_image::Matrix{L},
                  range_image::Matrix,
                  segments::Matrix{<:Integer}) where L <: Lab
    color_image = channelview(color_image)
    N = size(segments, 1) * size(segments, 2)
    K = 0
    unique_segments = unique(segments)
    params = Dict()
    sigma_w = mapwindow(color_std,
                        colorview(RGB, color_image),
                        (3,3))

    for i in unique_segments
        mask = segments .== i
        indices = findall(mask)
        valsIi = color_image[:, indices]
        valsRi = range_image[indices]
        S_star = erode(mask)
        params[i] = Dict(
            :stdI=> std(valsIi),
            :stdR=> std(valsRi),
            :meanI=>dropdims(
                mean(valsIi, dims=[2, 3]),
                dims=2),
            :meanR=>mean(valsRi),
            :n=>sum(mask),
            :sigma_t=>sum(sigma_w[findall(S_star)])/sum(S_star),
            :n_s_star=>sum(S_star)
        )
    end
    DIntraI = 0
    DInterI = 0
    DIntraD = 0
    DInterD = 0
    for i in unique_segments
        if (params[i][:n_s_star] > 0)
            K = K + 1
            DIntraI += max(params[i][:stdI] - params[i][:sigma_t], 0)*(params[i][:n]/N)
            DIntraD += params[i][:stdR]*(params[i][:n]/N)
            for j in unique_segments
                if i != j
                    DInterI += norm(params[i][:meanI] - params[j][:meanI])
                    DInterD += norm(params[i][:meanR] - params[j][:meanR])
                end
            end
        else
            @warn("Segment to small")
        end
    end
    DInterI = DInterI / (K*(K-1) + eps())
    DInterD = DInterD / (K*(K-1) + eps())
    QColor = (DInterI - DIntraI) / 2
    QDepth = (DInterD - DIntraD) / 2
    return QColor + 3*QDepth
end


text:
using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq

using Distributions: Exponential, Geometric, cdf
using Random: randsubseq
using StatsBase: sample

# use a bit array to identify which ones will fire and which will not.
struct queued_updates
  matches::Vector{ACSetTransformation{T}} where {T}
  valid::BitVector
  L::ACSetTransformation
  R::ACSetTransformation

  function queued_updates(matches::Vector{ACSetTransformation{T}}, L::ACSetTransformation, R::ACSetTransformation) where {T}
      # no events are scheduled to fire upon creation
      return new(matches, falses(length(matches)), L, R)
  end
end

# Bernoulli sample of matches in `m` (all same probability)
function sample_matches(m::queued_updates, r::AbstractFloat, Δt)
  p = cdf(Exponential(), r * Δt)
  fire = randsubseq(1:length(m.valid), p)
  m.valid[fire] .= true # these will fire
end

# Bernoulli sample of matches in `m` (unique probabilities)
function sample_matches(m::queued_updates, r::AbstractVector{R}, Δt) where {R <: AbstractFloat}
    p = cdf(Exponential(), r * Δt)
    runif = rand(length(m.valid))
    fire = Int64[]
    for i in 1:length(m.valid)
        if runif[i] < p[i]
            push!(fire, i)
        end
    end

    m.valid[fire] .= true
end

# update matches; return nothing if invalid (is there a way to do this in a type-stable way?)

"""
Convert a match L->G to a match L->H using a partial morphism G->H, if possible.
       L ========= L
     m ↓           ↓ m'
       G <-- K --> H
"""
function postcompose_partial(kg::ACSetTransformation, kh::ACSetTransformation, m::ACSetTransformation)
  m′ = Dict()
  # morphism (m) is all the functions between Obs in L -> G
  for (m_ob, m_fn) in pairs(components(m)) # m_ob: the ob, m_fn: the function for ob(L)->ob(G)
    m′_fn = Int[] # the new function in m′ between L -> H
    for x in collect(m_fn) # the actual function values of m_fn
      kx = findfirst(==(x), collect(kg[m_ob])) # find elem in K getting sent to the same thing as x
      if isnothing(kx) # it could be that thing was deleted, in which case there is no match
        return nothing
      else
        push!(m′_fn, kh[m_ob](kx)) # x now gets sent here in H
      end
    end
    m′[m_ob] = m′_fn # put the newly constructed fn for that ob into m′
  end
  return ACSetTransformation(dom(m), codom(kh); m′...)
end

# fire all events which are scheduled to fire this time step
function fire_events(state::ACSet, events::Vector{queued_updates})
  newstate = state
  # for each event type
  for i in 1:length(events)
    if sum(events[i].valid) == 0
      continue
    end
    # for each match
    for j in 1:length(events[i].matches)
        if events[i].valid[j] === false
            continue
        else
            # apply rewrite
            _, kg, _, kh = rewrite_match_maps(events[i].L, events[i].R, events[i].matches[j])
            newstate = codom(kh)
            # update remaining matches post rewrite
            if j < length(events[i].matches)
                for k in j+1:length(events[i].matches)
                # if the event wasn't going to happen anyway, continue
                if events[i].valid[k] === false
                    continue
                else
                    m = postcompose_partial(kg, kh, events[i].matches[k])
                    if isnothing(m)
                        events[i].valid[k] = false
                    else
                        events[i].matches[k] = m
                    end
                end
                end
            end
            # update remaining matches for other event types post rewrite
            if i < length(events)
                for k in i+1:length(events)
                    for l in 1:length(events[k].matches)
                    # if the event wasn't going to happen anyway, continue
                    if events[k].valid[l] === false
                        continue
                    else
                        m = postcompose_partial(kg, kh, events[k].matches[l])
                        if isnothing(m)
                            events[k].valid[l] = false
                        else
                            events[k].matches[l] = m
                        end
                    end
                    end
                end
            end
        end
    end
  end
  return newstate
end