text:
# Clayton, N. S. and Dally, J. and Gilbert, J. and Dickinson, A. (2005),
# Food Caching by Western Scrub-Jays (Aphelocoma californica) Is Sensitive to
# the Conditions at Recovery. Journal of Experimental Psychology: Animal
# Behavior Processes, 2:115-124

function _summarize(::Union{Experiment{:Clayton05_exp1}, Experiment{:Clayton05_exp2}},
                       data)
    combine(groupby(data, [:group, :trial, :action, :foodtype]),
	   df -> DataFrame(μ = mean(df.items), sem = sem(df.items)))
end
function statistical_tests(exp::Experiment{:Clayton05_exp1}, data)
    tests = Test[]
    id = "group × foodtype first trial"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :trial .== 1),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :foodtype,
                                     id = :id, data = :items),
                      exp.tests))
    id = "group × trial × foodtype"
    push!(tests, Test(id,
                      splitplotanova(data,
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = [:trial, :foodtype],
                                     id = :id, data = :items),
                      exp.tests))
    for g in groupby(data, :group)
        groupname = g.group[1]
        id = "Group $groupname"
        push!(tests, Test(id,
                          splitplotanova(g,
                                         exp.tests[id].locals,
                                         betweenfactors = [],
                                         withinfactors = [:trial, :foodtype],
                                         id = :id, data = :items),
                          exp.tests))
        if groupname == "consistent degrade"
            id = "Group $groupname: more worms than nuts in trial 1"
            push!(tests, Test(id,
                              anova(@where(g, :trial .== 1),
                                    @formula(items ~ foodtype + id), exp.tests[id].locals),
                              exp.tests))
            id = "Group $groupname: more nuts than worms in trials 3-6"
            push!(tests, Test(id,
                              splitplotanova(@where(g, (|).([:trial .== i
                                                             for i in 3:6]...)),
                                             exp.tests[id].locals,
                                             betweenfactors = [],
                                             withinfactors = [:trial, :foodtype],
                                             id = :id, data = :items),
                              exp.tests))
            for subg in groupby(g, :foodtype)
                foodtype = subg.foodtype[1]
                id = "Group consistent degrade: $foodtype"
                push!(tests, Test(id,
                                  anova(subg, @formula(items ~ trial + id), exp.tests[id].locals),
                                  exp.tests))
            end
        end
    end
    for g in groupby(data, :foodtype)
        id = "Foodtype $(g.foodtype[1])"
        push!(tests, Test(id,
                          splitplotanova(g,
                                         exp.tests[id].locals,
                                         betweenfactors = :group,
                                         withinfactors = :trial,
                                         id = :id, data = :items),
                          exp.tests))
    end
    TestCollection(tests)
end
function foodtypelabel(f)
    if f === Peanut
        "nuts"
    elseif f === Waxworm || f === Mealworm
        "worms"
    elseif f === Pineapple
        "pineapple"
    elseif f === Salami
        "salami"
    end
end
@inline function Clayton05_trial!(m, results, group, id, trial;
                                  longRIduration = 100u"hr",
                                  longRI = id > 2 && trial in (1, 4, 5, 7, 9) ||
                                           id <= 2 && trial in (2, 3, 6, 8, 10),
                                  foodtype1 = Peanut, foodtype2 = Waxworm,
                                  treatment = group in ("consistent degrade", "partial degrade") ?  "degrade" : group == "replenish" ? "none" : "pilfer",
                                  pilfer = group == "consistent pilfer" ||
                                           group == "partial pilfer" && longRI,
                                  degrade = group == "consistent degrade" ||
                                            group == "partial degrade" && longRI)
    add!(m, MaintenanceDiet) # Sunday 20:00
    wait!(m, 24u"hr")
    remove!(m, Any)          # Monday 20:00
    wait!(m, 12u"hr")
    add!(m, foodtype1, 50)      # Tuesday 8:00
    add!(m, foodtype2, 50)
    tray = Tray(1)
    add!(m, tray)
    wait!(m, 15u"minute")
    remove!(m, Any)
    if longRI
        add!(m, MaintenanceDiet)
        wait!(m, longRIduration - 4u"hr")
        remove!(m, Any)
    end
    wait!(m, 4u"hr")
    for foodtype in [foodtype1, foodtype2]
        push!(results, [group, trial, "cache",
                        foodtypelabel(foodtype),
                        id,
                        group in ("consistent degrade", "consistent pilfer"),
                        treatment,
                        countcache(tray, foodtype)])
    end
    degrade && degrade!(tray, foodtype2)
    pilfer && pilfer!(tray, foodtype2)
    remove!(m, Any)
    add!(m, tray)
    wait!(m, 10u"minute")
    remove!(m, Any)         # Tuesday 12:25 or Saturday 12:25
    add!(m, MaintenanceDiet)
    wait!(m, 7u"hr" + 35u"minute" + (longRI ? 1u"d" : longRIduration - 4u"hr"))
end
function run!(::Experiment{:Clayton05_exp1}, models)
    results = DataFrame(group = String[], trial = Int[], action = String[],
                        foodtype = String[],
                        id = Int[], consistent = Bool[], treatment = String[],
                        items = Int[])
    id = 0
    for group in ["replenish", "consistent degrade", "partial degrade"]
        for i in 1:4
            id += 1
            m = models[id]
            for trial in 1:6
                Clayton05_trial!(m, results, group, id, trial)
            end
        end
    end
    results
end

function statistical_tests(exp::Experiment{:Clayton05_exp2}, data)
    tests = Test[]
    id = "group × foodtype first trial"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :trial .== 1),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :foodtype,
                                     id = :id, data = :items),
                      exp.tests))
    id = "Group × Food Type × Trial for replenish, consistent degrade and partial degrade"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "consistent pilfer",
                                                  :group .!= "partial pilfer"),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = [:trial, :foodtype],
                                     id = :id, data = :items),
                      exp.tests))
    for foodtype in ("nuts", "worms")
        id = "Group × Trial for $foodtype"
        push!(tests, Test(id,
                          splitplotanova(@where(data, :group .!= "consistent pilfer",
                                                      :group .!= "partial pilfer",
                                                      :foodtype .== foodtype),
                                         exp.tests[id].locals,
                                         betweenfactors = :group,
                                         withinfactors = :trial,
                                         id = :id, data = :items),
                          exp.tests))
    end
    for g in groupby(data, [:group, :foodtype])
        groupname = g.group[1]
        foodtype = g.foodtype[1]
        (groupname == "consistent pilfer" || groupname == "partial pilfer") && continue
        id = "Change in caching $foodtype $groupname"
        push!(tests, Test(id,
                          anova(g, @formula(items ~ trial + id), exp.tests[id].locals),
                          exp.tests))
    end
    id = "no difference between worms cache by replenish and partial degrade in trials 3-8"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "consistent pilfer",
                                                  :group .!= "partial pilfer",
                                                  :group .!= "consistent degrade",
                                                  :foodtype .== "worms",
                                                  3 .<= :trial .<= 8),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :trial,
                                     id = :id, data = :items),
                      exp.tests))
    id = "more worms cache by consistent degrade than partial degrade in trials 3-8"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "consistent pilfer",
                                                  :group .!= "partial pilfer",
                                                  :group .!= "replenish",
                                                  :foodtype .== "worms",
                                                  3 .<= :trial .<= 8),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :trial,
                                     id = :id, data = :items),
                      exp.tests))
    id = "more nuts cache by replenish than partial degrade in trials 4-8"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "consistent pilfer",
                                                  :group .!= "partial pilfer",
                                                  :group .!= "consistent degrade",
                                                  :foodtype .== "nuts",
                                                  4 .<= :trial .<= 8),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :trial,
                                     id = :id, data = :items),
                      exp.tests))
    id = "more nuts cache by consistent degrade than replenish in trials 4-8"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "consistent pilfer",
                                                  :group .!= "partial pilfer",
                                                  :group .!= "partial degrade",
                                                  :foodtype .== "nuts",
                                                  4 .<= :trial .<= 8),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = :trial,
                                     id = :id, data = :items),
                      exp.tests))
    id = "First trial pilfer"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "replenish",
                                                  :trial .== 1),
                                     exp.tests[id].locals,
                                     betweenfactors = :group,
                                     withinfactors = [:foodtype],
                                     id = :id, data = :items),
                      exp.tests))
    id = "compare pilfer, degrade"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :group .!= "replenish"),
                                     exp.tests[id].locals,
                                     betweenfactors = [:treatment, :consistent],
                                     withinfactors = [:foodtype, :trial],
                                     id = :id, data = :items),
                      exp.tests))
    TestCollection(tests)
end
function run!(::Experiment{:Clayton05_exp2}, models)
    results = DataFrame(group = String[], trial = Int[], action = String[],
                        foodtype = String[],
                        id = Int[], consistent = Bool[], treatment = String[],
                        items = Int[])
    id = 0
    for group in ["replenish", "consistent degrade", "partial degrade",
                  "consistent pilfer", "partial pilfer"]
        for i in 1:(group in ["consistent pilfer", "partial pilfer"] ? 5 : 4)
            id += 1
            m = models[id]
            for trial in 1:8
                Clayton05_trial!(m, results, group, id, trial,
                                 longRIduration = 28u"hr")
            end
        end
    end
    results
end

function _summarize(::Experiment{:Clayton05_exp3}, results)
    res = combine(groupby(results, [:group, :trial, :action, :foodtype]),
       d -> DataFrame(μ = mean(d.items), n = length(d.items)))
    for foodtype in ("pineapple", "salami")
        data = @where(results, :group .== "partial degrade", :treatment .== foodtype)
        push!(res, ["partial degrade", missing, foodtype, "cache",
                        mean(data.items), length(data.items)])
    end
    res
end
# key:      here - paper
#          group - group
#      treatment - foodtype
#          trial - trial
#       foodtype - degrade
function statistical_tests(exp::Experiment{:Clayton05_exp3}, data)
    tests = Test[]
    id = "overall"
    push!(tests, Test(id,
                      splitplotanova(data,
                                     exp.tests[id].locals,
                                     betweenfactors = [:group, :treatment],
                                     withinfactors = [:foodtype, :trial],
                                     id = :id, data = :items),
                      exp.tests))
    id = "trial 1"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :trial .== 1),
                                     exp.tests[id].locals,
                                     betweenfactors = [:group, :treatment],
                                     withinfactors = :foodtype, id = :id,
                                     data = :items),
                      exp.tests))
    for g in groupby(data, :group)
        group = g.group[1]
        push!(tests, Test(group,
                          splitplotanova(g,
                                         exp.tests[group].locals,
                                         betweenfactors = :treatment,
                                         withinfactors = [:foodtype, :trial],
                                         id = :id, data = :items),
                          exp.tests))
        if group == "consistent degrade"
            for subg in groupby(g, :trial)
                id = "$group trial $(subg.trial[1])"
                push!(tests, Test(id,
                                  splitplotanova(subg,
                                                 exp.tests[id].locals,
                                                 betweenfactors = :treatment,
                                                 withinfactors = :foodtype,
                                                 id = :id, data = :items),
                                  exp.tests))
            end
        end
    end
    TestCollection(tests)
end
function run!(::Experiment{:Clayton05_exp3}, models)
    results = DataFrame(group = String[], trial = Union{Int,Missing}[],
                        action = String[],
                        foodtype = String[],
                        id = Int[], consistent = Bool[], treatment = String[],
                        items = Int[])
    id = 0
    for group in ["consistent degrade", "partial degrade"]
        n = group == "consistent degrade" ? 6 : 8
        for i in 1:n
            id += 1
            m = models[id]
            for trial in 1:10
                Clayton05_trial!(m, results, group, id, trial,
                                 longRIduration = 28u"hr",
                                 foodtype1 = i <= n/2 ? Pineapple : Salami,
                                 foodtype2 = i <= n/2 ? Salami : Pineapple,
                                 treatment = i <= n/2 ? "salami" : "pineapple")
            end
        end
    end
    for i in 1:nrow(results)
        results.foodtype[i] = results.foodtype[i] == results.treatment[i] ? "degrading" : "nondegrading"
    end
    results
end

function _summarize(::Experiment{:Clayton05_exp4}, results)
    combine(groupby(results, [:trial, :action, :foodtype]),
        d -> DataFrame(μ = mean(d.items), n = length(d.items)))
end
function statistical_tests(exp::Experiment{:Clayton05_exp4}, data)
    tests = Test[]
    id = "overall"
    push!(tests, Test(id,
                      splitplotanova(data,
                                     exp.tests[id].locals,
                                     betweenfactors = [],
                                     withinfactors = [:foodtype, :trial],
                                     id = :id, data = :items),
                      exp.tests))
    id = "trial 4"
    push!(tests, Test(id,
                      anova(@where(data, :trial .== 4),
                            @formula(items ~ foodtype + id), exp.tests[id].locals),
                      exp.tests))
    for g in groupby(data, :foodtype)
        push!(tests, Test(g.foodtype[1],
                          anova(g, @formula(items ~ trial + id),
                                exp.tests[g.foodtype[1]].locals),
                          exp.tests))
    end
    push!(tests, Test("worms 1-2",
                      is_significantly_bigger(@where(data, :foodtype .== "worms"),
                                              condbigger = 1, condsmaller = 2,
                                              conditionkey = :trial,
                                              valuekey = :items, paired = true),
                      exp.tests))
    TestCollection(tests)
end
function run!(::Experiment{:Clayton05_exp4}, models)
    results = DataFrame(group = String[], trial = Union{Int,Missing}[],
                        action = String[],
                        foodtype = String[],
                        id = Int[], consistent = Bool[], treatment = String[],
                        items = Int[])
    for i in 1:6
        m = models[i]
        for trial in 1:4
            # TODO: how to deal with test at the end of trial 4?
            Clayton05_trial!(m, results, "", i, trial,
                             longRI = true, degrade = true,
                             longRIduration = trial in (1, 3) ? 76u"hr" : 52u"hr",
                             foodtype2 = Mealworm)
        end
    end
    results
end


text:
import JSON
export BrainFlowModelParams

abstract type BrainFlowMetric end
struct Relaxation <: BrainFlowMetric end
struct Concentration <: BrainFlowMetric end
Base.Integer(::Relaxation) = 0
Base.Integer(::Concentration) = 1
Base.convert(::Type{BrainFlowMetric}, s::AbstractString) = BrainFlowMetric(Val(Symbol(lowercase(s))))
BrainFlowMetric(::Val{:relaxation}) = Relaxation()
BrainFlowMetric(::Val{:concentration}) = Concentration()

abstract type BrainFlowClassifier end
struct Regression <: BrainFlowClassifier end
struct Knn <: BrainFlowClassifier end
struct Svm <: BrainFlowClassifier end
struct Lda <: BrainFlowClassifier end
Base.Integer(::Regression) = 0
Base.Integer(::Knn) = 1
Base.Integer(::Svm) = 2
Base.Integer(::Lda) = 3
Base.convert(::Type{BrainFlowClassifier}, s::AbstractString) = BrainFlowClassifier(Val(Symbol(lowercase(s))))
BrainFlowClassifier(::Val{:regression}) = Regression()
BrainFlowClassifier(::Val{:knn}) = Knn()
BrainFlowClassifier(::Val{:svm}) = Svm()
BrainFlowClassifier(::Val{:lda}) = Lda()

@Base.kwdef mutable struct BrainFlowModelParams
    metric::BrainFlowMetric = Relaxation()
    classifier::BrainFlowClassifier = Regression()
    file::String = ""
    other_info::String = ""
end

function JSON.json(params::BrainFlowModelParams)
    d = Dict(
        "metric" => Integer(params.metric),
        "classifier" => Integer(params.classifier),
        "file" => params.file,
        "other_info" => params.other_info,
        )
    return JSON.json(d)
end

@brainflow_rethrow function prepare(params::BrainFlowModelParams)
    input_json = JSON.json(params)
    ccall((:prepare, ML_MODULE_INTERFACE), Cint, (Ptr{UInt8},), input_json)
    return
end

@brainflow_rethrow function release(params::BrainFlowModelParams)
    input_json = JSON.json(params)
    ccall((:release, ML_MODULE_INTERFACE), Cint, (Ptr{UInt8},), input_json)
    return
end

@brainflow_rethrow function predict(data, params::BrainFlowModelParams)
    input_json = JSON.json(params)
    val = Vector{Float64}(undef, 1)
    ccall((:predict, ML_MODULE_INTERFACE), Cint, (Ptr{Float64}, Cint, Ptr{Float64}, Ptr{UInt8}),
        data, length(data), val, input_json)
    value = val[1]
    return value
end

ext: jl


text:
# This test reproduces a case, where all starting strategies failed to
# produce a reasonable estimate. Estimated τ diverges.
module DebugTest

using DataFrames, LinearAlgebra, Random, Test, WiSER

function createvlmm(t, k, j)
    p  = 5    # number of fixed effects, including intercept
    q  = 2    # number of random effects, including intercept
    l  = 5    # number of WS variance covariates, including intercept
    q◺ = ◺(q)

    # true parameter values
    global βtrue = [0.1; 6.5; -3.5; 1.0; 5]
    global τtrue = [0.0; 0.5; -0.2; 0.5; 0.0]
    Σγ    = Matrix(Diagonal([2.0; 1.2]))
    δγω   = [0.2; 0.1]
    σω    = [1.0]

    Σγω   = [Σγ δγω; δγω' σω]
    Lγω   = cholesky(Symmetric(Σγω), check = false).L
    global Lγ    = Lγω[1:q, 1:q]
    lγω   = Lγω[q + 1, 1:q]
    lω    = Lγω[q + 1, q + 1]

    vechLγ = vech(Lγ)
    # generate data
    γω = Vector{Float64}(undef, q + 1)
    z  = similar(γω) # hold vector of iid std normal

    intervals = zeros(p + l, 2) #hold intervals
    curcoverage = zeros(p + l) #hold current coverage resutls
    trueparams = [βtrue; τtrue] #hold true parameters

    #simulation parameters
    samplesizes = collect(1000:1000:6000)
    ns = [10; 25; 50; 100; 1000]

    m = samplesizes[t]
    ni = ns[k] # number of observations per individual
    obsvec = Vector{WSVarLmmObs{Float64}}(undef, m)
    println("rep $j obs per person $ni samplesize $m")
    Random.seed!(j + 100000k + 1000t)
    for i in 1:m
        # first column intercept, remaining entries iid std normal
        X = Matrix{Float64}(undef, ni, p)
        X[:, 1] .= 1
        @views Distributions.rand!(Normal(), X[:, 2:p])
        # first column intercept, remaining entries iid std normal
        Z = Matrix{Float64}(undef, ni, q)
        Z[:, 1] .= 1
        @views Distributions.rand!(Normal(), Z[:, 2:q])
        # first column intercept, remaining entries iid std normal
        W = Matrix{Float64}(undef, ni, l)
        W[:, 1] .= 1
        @views Distributions.rand!(Normal(), W[:, 2:l])
        # generate random effects: γω = Lγω * z
        mul!(γω, Lγω, Distributions.rand!(Normal(), z))
        # generate y
        μy = X * βtrue + Z * γω[1:q]
        @views vy = exp.(W * τtrue .+ dot(γω[1:q], lγω) .+ γω[end])
        y = rand(MvNormal(μy, Diagonal(vy)))
        # form a VarLmmObs instance
        obsvec[i] = WSVarLmmObs(y, X, Z, W)
    end
    # form VarLmmModel
    vlmm = WSVarLmmModel(obsvec);
    return vlmm
end

vlmm1 = createvlmm(1, 1, 203) # trouble case
# vlmm1 = createvlmm(2, 1, 233)

### If you revert back to initialize intercept only these cases will also fail:
# ts = [2,   2,   2,   4,   1,   3,   2, 1,   1,   2]
# ks = [1,   1,   1,   1,   2,   1,   2, 1,   1,   1]
# js = [35, 53, 109, 148, 168, 174, 100, 7, 203, 233]

# solver = Ipopt.IpoptSolver(print_level=0)
solver = Ipopt.IpoptSolver(print_level=0, mehrotra_algorithm="yes", max_iter=100)
# solver = NLopt.NLoptSolver(algorithm = :LN_BOBYQA,
#     ftol_rel = 1e-12, ftol_abs = 1e-8, maxeval = 10000)
# solver = NLopt.NLoptSolver(algorithm = :LD_SLSQP, maxeval = 4000)
# solver = NLopt.NLoptSolver(algorithm = :LD_MMA, maxeval = 4000)
# solver = NLopt.NLoptSolver(algorithm = :LD_LBFGS, maxeval = 4000)

@info "starting point by LS init_ls!"
WiSER.init_ls!(vlmm1)
println("β"); display(vlmm1.β); println()
println("τ"); display(vlmm1.τ); println()
println("Lγ"); display(vlmm1.Lγ); println()
@show norm(vlmm1.β - βtrue)
@show norm(vlmm1.τ[2:end] - τtrue[2:end])
@show norm(vlmm1.Lγ - Lγ)

@info "unweighted NLS fit (start from LS)"
@time init_mom!(vlmm1, solver)
println("β"); display(vlmm1.β); println()
println("τ"); display(vlmm1.τ); println()
println("Lγ"); display(vlmm1.Lγ); println()
@show norm(vlmm1.β - βtrue)
@show norm(vlmm1.τ[2:end] - τtrue[2:end])
@show norm(vlmm1.Lγ - Lγ)

# @info "weighted NLS fit (start from true τ and LS β, Σγ)"
# init_ls!(vlmm1)
# # vlmm1.τ .= [-0.1; 0.5; -0.2; 0.5; 0.0]
# # vlmm1.τ .= [0.27, 0.14, -0.04, 0.13, -0.02] # LS (x5 works)
# # vlmm1.τ .= [2.09, 0.13, -0.019, 0.15, -0.0007] # log-LS (x5 works)
# fit!(vlmm1, solver; init = vlmm1, runs = 2)
# println("β"); display(vlmm1.β); println()
# println("τ"); display(vlmm1.τ); println()
# println("Lγ"); display(vlmm1.Lγ); println()
# show(vlmm1)

@info "weighted fit (start from LS)"
@time fit!(vlmm1, solver; init = init_ls!(vlmm1), runs = 2)
println("β"); display(vlmm1.β); println()
println("τ"); display(vlmm1.τ); println()
println("Lγ"); display(vlmm1.Lγ); println()
@show norm(vlmm1.β - βtrue)
@show norm(vlmm1.τ[2:end] - τtrue[2:end])
@show norm(vlmm1.Lγ - Lγ)

@info "weighted fit (start from MoM)"
@time fit!(vlmm1, solver; init = init_mom!(vlmm1, solver), runs = 2)
println("β"); display(vlmm1.β); println()
println("τ"); display(vlmm1.τ); println()
println("Lγ"); display(vlmm1.Lγ); println()
@show norm(vlmm1.β - βtrue)
@show norm(vlmm1.τ[2:end] - τtrue[2:end])
@show norm(vlmm1.Lγ - Lγ)

end


text:
"""
TODO: Description of the current population model.
"""
module population
importall common
using   core, _deme_, utility, species, breeder
export  Population, PopulationDefinition

"""
The base type used by all population definitions.
"""
abstract PopulationDefinition

"""
Used to hold the contents of a population, where a population is modelled
as a collection of demes (which may, or may not, be of different species).
"""
type Population
  demes::Vector{Deme}

  Population() = new([])
  Population(dl::Vector{Deme}) = new(dl)
end

"""
Produces the offspring for a given population at each generation.
"""
breed!(p::Population) = for deme in p.demes; breed!(deme); end

"""
Performs any necessary fitness scaling and processing on each of the
individuals within a given population.
"""
scale!(p::Population) = for deme in p.demes; scale!(deme); end
#scale!(d::Deme) = scale!(d.species.fitness, contents(d))

"""
TODO: IMPLEMENT SCALING.
"""
scale!(d::Deme) = true

# Load each of the population types.
include("population/simple.jl")
include("population/complex.jl")
end


text: println("I'm excited to learn Julia!")

max_stars_repo_name: PushpneetSingh/Hello-world

max_issues_repo_name: PushpneetSingh/Hello-world

max_forks_repo_name: PushpneetSingh/Hello-world


text:
#= This file is a part of URIValidator. License is MIT:
https://github.com/rafpyprog/URIValidator.jl/blob/master/LICENSE.md =#

using Base.Test
using URIValidator

splited = Dict("scheme"=>"http",
               "authority"=>"julialang.org",
               "path"=>"",
               "query"=>"",
               "fragment"=>"")

function test_invalid_argument_exception()
    n = deepcopy(splited)
    pop!(n, "fragment", nothing) #remove one component
    result = URIComponents(n)
end

function test_uri_component_type()
    result = URIComponents(splited)
    return isa(result, URIComponents)
end

function test_components_fields()
    result = URIComponents(splited)
    check_components = (
        isequal(result.scheme, splited["scheme"]),
        isequal(result.authority, splited["authority"]),
        isequal(result.path, splited["path"]),
        isequal(result.query, splited["query"]),
        isequal(result.fragment, splited["fragment"]),
        )
    return all(check_components)
end

function test_components_tuple()
    test_case = URIComponents(splited)
    components = test_case.components
    return isequal(length(components), 5)
end

@testset "URIComponents" begin
    @test test_uri_component_type()
    @test test_components_fields()
    @test test_components_tuple()
    @test_throws ArgumentError test_invalid_argument_exception()
end


text:
using Random

@testset "Testing random_mc.jl" begin
    @testset "Test random_markov_chain" begin
        n, k = 5, 3
        mc_dicts = (Dict("P" => random_markov_chain(n).p, "k" => n),
                    Dict("P" => random_markov_chain(n, k).p, "k" => k))
        for d in mc_dicts
            P = d["P"]
            @test size(P) == (n, n)
            @test all(x->(count(!iszero, x)==d["k"]),
                      [P[i, :] for i in 1:size(P)[1]]) == true
        end

        seed = 1234
        rngs = [MersenneTwister(seed) for i in 1:2]
        mcs = random_markov_chain.(rngs, n, k)
        @test mcs[2].p == mcs[1].p
    end

    @testset "Test random_stochastic_matrix" begin
        n, k = 5, 3
        Ps = (random_stochastic_matrix(n), random_stochastic_matrix(n, k))
        for P in Ps
            @test all(P .>= 0) == true
            @test all(x->isapprox(sum(x), 1),
                      [P[i, :] for i in 1:size(P)[1]]) == true
        end

        seed = 1234
        rngs = [MersenneTwister(seed) for i in 1:2]
        Ps = random_stochastic_matrix.(rngs, n, k)
        @test Ps[2] == Ps[1]
    end

    @testset "Test random_stochastic_matrix with k=1" begin
        n, k = 3, 1
        P = random_stochastic_matrix(n, k)
        @test all((P .== 0) .| (P .== 1))
        @test all(x->isequal(sum(x), 1),
                  [P[i, :] for i in 1:size(P)[1]]) == true
    end

    @testset "Test errors properly thrown" begin
        # n <= 0
        @test_throws ArgumentError random_markov_chain(0)

        # k <= 0
        @test_throws ArgumentError random_markov_chain(2, 0)

        # k > n
        @test_throws ArgumentError random_markov_chain(2, 3)
    end

    @testset "Test random_discrete_dp" begin
        num_states, num_actions = 5, 4
        num_sa = num_states * num_actions
        k = 3
        ddp = random_discrete_dp(num_states, num_actions; k=k)

        # Check shapes
        @test size(ddp.R) == (num_states, num_actions)
        @test size(ddp.Q) == (num_states, num_actions, num_states)

        # Check ddp.Q[:, a, :] is a stochastic matrix for all actions `a`
        @test all(ddp.Q .>= 0) == true
        for a in 1:num_actions
            P = reshape(ddp.Q[:, a, :], (num_states, num_states))
            @test all(x->isapprox(sum(x), 1),
                      [P[i, :] for i in 1:size(P)[1]]) == true
        end

        # Check number of nonzero entries for each state-action pair
        @test sum(ddp.Q .> 0, dims = 3) ==
            ones(Int, (num_states, num_actions, 1)) * k

        seed = 1234
        rngs = [MersenneTwister(seed) for i in 1:2]
        ddps = random_discrete_dp.(rngs, num_states, num_actions)
        @test ddps[2].R == ddps[1].R
        @test ddps[2].Q == ddps[1].Q
        @test ddps[2].beta == ddps[1].beta
    end

    @testset "Test random_probvec" begin
        k = 5
        x = random_probvec(k)
        @test size(x) == (k,)
    end
end  # @testset


text:
module Jieba
using PyCall

export jieba, analyse, pseg
export cut, cut_for_search

const jieba, analyse, pseg = PyNULL(), PyNULL(), PyNULL()

# https://github.com/JuliaPy/PyCall.jl#using-pycall-from-julia-modules
function __init__()
    copy!(jieba, pyimport("jieba"))
    copy!(analyse, pyimport("jieba.analyse"))
    copy!(pseg, pyimport("jieba.posseg"))
end

cut(msg; use_paddle = false, cut_all = false, HMM = false) = jieba.lcut(
    msg, use_paddle = use_paddle, cut_all = cut_all, HMM = HMM
)
cut_for_search(msg; HMM = false) = jieba.lcut_for_search(msg, HMM = HMM)

end

max_stars_repo_name: lotcher/Jieba

max_issues_repo_name: lotcher/Jieba

max_forks_repo_name: lotcher/Jieba


text:
using Documenter, IterTools

makedocs(
    modules = [IterTools],
    format = :html,
    sitename = "IterTools",
    pages = Any[
        "Introduction" => "index.md",
        "Function index" => "functionindex.md"
        ])

deploydocs(
    repo = "github.com/JuliaCollections/IterTools.jl.git",
    target = "build",
    julia  = "1.0",
    osname = "linux",
    deps = nothing,
    make = nothing
)


text:
# Plotting functionalities.

using PyPlot
using Dierckx: Spline2D

export plotdata, plotlogdata, plot, scatter, contour, contourf, plot_surface,
   tricontourf, plot_trisurf, streamplot, streamslice, quiver, cutplot

include("utility.jl")

@static if matplotlib.__version__ >= "3.3"
   matplotlib.rc("image", cmap="turbo") # set default colormap
end

"""
    plotlogdata(data, head, func; plotmode="line")

Plot information from log file.
# Input arguments
- `data::Array`: output data.
- `head::NamedTuple`: header info.
- `func::String`: variables for plotting.
- `plotmode::String`: type of plotting ["line","scatter"].
"""
function plotlogdata(data, head::NamedTuple, func::AbstractString; plotmode="line")

   vars     = split(func)
   plotmode = split(plotmode)

   for (ivar, var) in enumerate(vars)
      varIndex_ = findfirst(x->x==lowercase(var), lowercase.(head.variables))
      isnothing(varIndex_) && error("$(var) not found in file header variables!")

      figure()
      if plotmode[ivar] == "line"
         plot(data[1,:],data[varIndex_,:])
      elseif plotmode[ivar] == "scatter"
         scatter(data[1,:],data[varIndex_,:])
      else
         throw(ArgumentError("unknown plot mode $(plotmode[ivar])!"))
      end
      xlabel(head.variables[1])
      ylabel(head.variables[varIndex_])
      title("log file data")
   end

end

"""
    plotdata(data, func, kwargs)

Plot the variable from SWMF output.

`plotdata(data, "p", plotmode="contbar")`

`plotdata(data, "p", plotmode="grid")`

`plotdata(data, func, plotmode="trimesh",plotrange=plotrange, plotinterval=0.2)`

# Input arguments
- `data::Data`: output data.
- `func::String`: variables for plotting.
- `plotmode::String`: (optional) type of plotting ["cont","contbar"]...
- `plotrange::Vector`: (optional) range of plotting.
- `plotinterval`: (optional) interval for interpolation.
- `level`: (optional) level of contour.
- `innermask`: (optional) Bool for masking a circle at the inner boundary.
- `dir`: (optional) 2D cut plane orientation from 3D outputs ["x","y","z"].
- `sequence`: (optional) sequence of plane from - to + in that direction.
- `multifigure`: (optional) 1 for multifigure display, 0 for subplots.
- `verbose`: (optional) display additional information.
- `density`: (optional) density for streamlines.
- `stride`: (optional) quiver strides in number of cells.
Right now this can only deal with 2D plots or 3D cuts. Full 3D plots may be supported in the
future.
"""
function plotdata(data::Data, func::AbstractString; dir="x", plotmode="contbar",
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, sequence=1, multifigure=true,
   getrangeonly=false, level=0, innermask=false, verbose=false, density=1.0,
   stride=10)

   x, w = data.x, data.w
   plotmode = split(plotmode)
   vars     = split(func)
   ndim     = data.head.ndim
   nvar     = length(vars)

   if verbose || getrangeonly
      @info "============ PLOTTING PARAMETERS ==============="
      @info "wnames = $(data.head.wnames)"
      wmin = Vector{Float64}(undef,nvar)
      wmax = Vector{Float64}(undef,nvar)
      # Display min & max for each variable
      for (ivar,var) in enumerate(vars)
         if occursin(";",var) continue end # skip the vars for streamline
         varIndex_ = findindex(data, var)
         if ndim == 1
            wmin[ivar] = minimum(w[:,varIndex_])
            wmax[ivar] = maximum(w[:,varIndex_])
         elseif ndim == 2
            wmin[ivar] = minimum(w[:,:,varIndex_])
            wmax[ivar] = maximum(w[:,:,varIndex_])
         end
         println("Min & Max value for $(var) :$(wmin[ivar])",", $(wmax[ivar])")
      end
      if getrangeonly return wmin, wmax end
   end

   ## plot multiple variables with same plotmode
   if length(plotmode) < nvar
      [push!(plotmode, plotmode[i]) for i = 1:nvar-length(plotmode)]
   end

   ## Plot
   if ndim == 1
      for (ivar,var) in enumerate(vars)
         varIndex_ = findindex(data, var)
         if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
         if !occursin("scatter",plotmode[ivar])
            plot(x, w[:,varIndex_])
         else
            scatter(x, w[:,varIndex_])
         end
         if occursin("grid", plotmode[ivar])
            grid(true)
         end
         xlabel("x"); ylabel("$(var)")
         dim = [0.125, 0.013, 0.2, 0.045]
         str = @sprintf "it=%d, time=%4.2f" data.head.it data.head.time
         at = matplotlib.offsetbox.AnchoredText(str,
            loc="lower left", prop=Dict("size"=>8), frameon=true,
            bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
         at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
         ax.add_artist(at)
      end
   elseif ndim == 2
      for (ivar,var) in enumerate(vars)
         occursin("over", plotmode[ivar]) && (multifigure = false)
         if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end
         if !occursin(";",var)
            varIndex_ = findindex(data, var)
         end

         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar",
            "contlog","contbarlog")

            Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; griddim=2, innermask)

            # More robust method needed!
            if plotmode[ivar] ∈ ["contbar", "contbarlog"]
               if level == 0
                  c = contourf(Xi, Yi, Wi)
               else
                  c = contourf(Xi, Yi, Wi, level)
               end
            elseif plotmode[ivar] ∈ ["cont", "contlog"]
               c = contour(Xi, Yi, Wi)
            elseif plotmode[ivar] ∈ ["surfbar", "surfbarlog"]
               c = plot_surface(Xi, Yi, Wi)
            end

            occursin("bar", plotmode[ivar]) && colorbar()
            occursin("log", plotmode[ivar]) &&
            ( c.locator = matplotlib.ticker.LogLocator() )
            title(data.head.wnames[varIndex_])

         elseif plotmode[ivar] ∈ ("trimesh","trisurf","tricont","tristream")
            X = vec(x[:,:,1])
            Y = vec(x[:,:,2])
            W = vec(w[:,:,varIndex_])

            # This needs to be modified!!!
            if !all(isinf.(plotrange))
               xyIndex = X .> plotrange[1] .& X .< plotrange[2] .&
                  Y .> plotrange[3] .& Y .< plotrange[4]
               X = X[xyIndex]
               Y = Y[xyIndex]
               W = W[xyIndex]
            end

            if plotmode[ivar] == "trimesh"
               triang = matplotlib.tri.Triangulation(X, Y)
               c = ax.triplot(triang)
            elseif plotmode[ivar] == "trisurf"
               c = ax.plot_trisurf(X, Y, W)
            elseif plotmode[ivar] == "tricont"
               c = ax.tricontourf(X, Y, W)
               fig.colorbar(c,ax=ax)
            elseif plotmode[ivar] == "tristream"
               throw(ArgumentError("tristream not yet implemented!"))
            end

            title(data.head.wnames[varIndex_])

         elseif plotmode[ivar] ∈ ("stream","streamover")
            VarStream  = split(var,";")
            VarIndex1_ = findindex(data, VarStream[1])
            VarIndex2_ = findindex(data, VarStream[2])

            if data.head.gencoord # Generalized coordinates
               xrange = @view x[:,:,1]
               yrange = @view x[:,:,2]
               if any(isinf.(plotrange))
                  if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end
                  if plotrange[2] ==  Inf plotrange[2] = maximum(xrange) end
                  if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end
                  if plotrange[4] ==  Inf plotrange[4] = maximum(yrange) end
               end

               # Create grid values first.
               xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
               yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)

               # The PyCall here can be potentially replaced with Spline2D.
               # Perform linear interpolation of the data (x,y) on grid(xi,yi)
               triang = matplotlib.tri.Triangulation(xrange[:], yrange[:])

               Xi, Yi = meshgrid(xi, yi)

               W = w[:,1,VarIndex1_]
               interpolator = matplotlib.tri.LinearTriInterpolator(triang, W)
               v1 = interpolator(Xi, Yi)

               W = w[:,1,VarIndex2_]
               interpolator = matplotlib.tri.LinearTriInterpolator(triang, W)
               v2 = interpolator(Xi, Yi)

            else # Cartesian coordinates
               xrange = @view x[:,1,1]
               yrange = @view x[1,:,2]
               if all(isinf.(plotrange))
                  Xi, Yi = meshgrid(xrange, yrange)
                  v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'
               else
                  if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end
                  if plotrange[2] ==  Inf plotrange[2] = maximum(xrange) end
                  if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end
                  if plotrange[4] ==  Inf plotrange[4] = maximum(yrange) end

                  w1, w2 = w[:,:,VarIndex1_], w[:,:,VarIndex2_]

                  xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
                  yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)

                  Xi, Yi = meshgrid(xi, yi)

                  spline = Spline2D(xrange, yrange, w1)
                  v1 = spline(Xi[:], Yi[:])
                  v1 = reshape(v1, size(Xi))

                  spline = Spline2D(xrange, yrange, w2)
                  v2 = spline(Xi[:], Yi[:])
                  v2 = reshape(v2, size(Xi))
               end
            end

            s = streamplot(Xi, Yi, v1, v2; color="w", linewidth=1.0, density)

         elseif occursin("quiver", plotmode[ivar])
            VarQuiver  = split(var, ";")
            VarIndex1_ = findindex(data, VarQuiver[1])
            VarIndex2_ = findindex(data, VarQuiver[2])

            X, Y = x[:,1,1], x[1,:,2]
            v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'

            @views Xq, Yq = X[1:stride:end], Y[1:stride:end]
            v1q = @view v1[1:stride:end, 1:stride:end]
            v2q = @view v2[1:stride:end, 1:stride:end]

            q = quiver(Xq, Yq, v1q, v2q, color="w")

         elseif occursin("grid", plotmode[ivar])
            # This does not take subdomain plot into account!
            X, Y = x[:,:,1], x[:,:,2]
            scatter(X, Y, marker=".", alpha=0.6)
            title("Grid illustration")
         else
            throw(ArgumentError("unknown plot mode: $(plotmode[ivar])"))
         end

         xlabel(data.head.variables[1]); ylabel(data.head.variables[2])
         dim = [0.125, 0.013, 0.2, 0.045]
         str = @sprintf "it=%d, time=%4.2f" data.head.it data.head.time
         at = matplotlib.offsetbox.AnchoredText(str,
            loc="lower left", prop=Dict("size"=>8), frameon=true,
            bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
         at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
         ax.add_artist(at)
         # recover status
         occursin("over", plotmode[ivar]) && (multifigure = true)
      end

   else # 2D cut from 3D output; now only for Cartesian output
      X = @view x[:,:,:,1]
      Y = @view x[:,:,:,2]
      Z = @view x[:,:,:,3]
      for (ivar,var) in enumerate(vars)
         if plotmode[ivar] ∈ ("surf","surfbar","surfbarlog","cont","contbar", "contlog",
               "contbarlog")

            varIndex_ = findindex(data, var)

            if ivar == 1 || multifigure fig, ax = subplots() else ax = gca() end

            W = w[:,:,:,varIndex_]

            if dir == "x"
               cut1 = @view X[sequence,:,:]
               cut2 = @view Y[sequence,:,:]
               W    = @view W[sequence,:,:]
            elseif dir == "y"
               cut1 = @view X[:,sequence,:]
               cut2 = @view Z[:,sequence,:]
               W    = @view W[:,sequence,:]
            elseif dir == "z"
               cut1 = @view X[:,:,sequence]
               cut2 = @view Y[:,:,sequence]
               W    = @view W[:,:,sequence]
            end
         elseif plotmode[ivar] ∈ ("stream","streamover")
            varStream  = split(var,";")
            VarIndex1_ = findindex(data, varStream[1])
            VarIndex2_ = findindex(data, varStream[2])

            v1 = @view w[:,:,:,VarIndex1_]
            v2 = @view w[:,:,:,VarIndex2_]

            if dir == "x"
               cut1 = @view Y[sequence,:,:]
               cut2 = @view Z[sequence,:,:]
               v1   = v1[sequence,:,:]'
               v2   = v2[sequence,:,:]'
            elseif dir == "y"
               cut1 = @view X[:,sequence,:]
               cut2 = @view Z[:,sequence,:]
               v1   = v1[:,sequence,:]'
               v2   = v2[:,sequence,:]'
            elseif dir == "z"
               cut1 = @view X[:,:,sequence]
               cut2 = @view Y[:,:,sequence]
               v1   = v1[:,:,sequence]'
               v2   = v2[:,:,sequence]'
            end
            cut1, cut2 = cut1', cut2'
         end

         if !all(isinf.(plotrange))
            cut1, cut2, v1, v2 = subsurface(cut1, cut2, v1, v2, plotrange)
         end

         if plotmode[ivar] ∈ ("surf", "surfbar", "surfbarlog", "cont", "contbar", "contlog",
               "contbarlog")
            if level == 0
               c = ax.contourf(cut1, cut2, W)
            else
               c = ax.contourf(cut1, cut2, W, level)
            end
            fig.colorbar(c, ax=ax)
            title(data.head.wnames[varIndex_])

         elseif plotmode[ivar] ∈ ("stream", "streamover")
            xi = range(cut1[1,1], stop=cut1[1,end],
               step=(cut1[1,end] - cut1[1,1]) / (size(cut1,2) - 1))
            yi = range(cut2[1,1], stop=cut2[end,1],
               step=(cut2[end,1] - cut2[1,1]) / (size(cut2,1) - 1))

            Xi, Yi = meshgrid(xi, yi)

            s = streamplot(Xi, Yi, v1, v2; color="w", linewidth=1.0, density)
         end

         if dir == "x"
            xlabel("y"); ylabel("z")
         elseif dir == "y"
            xlabel("x"); ylabel("z")
         elseif dir == "z"
            xlabel("x"); ylabel("y")
         end

         ax = gca()
         dim = [0.125, 0.013, 0.2, 0.045]
         str = @sprintf "it=%d, time=%4.2f" data.head.it data.head.time
         at = matplotlib.offsetbox.AnchoredText(str,
            loc="lower left", prop=Dict("size"=>8), frameon=true,
            bbox_to_anchor=(0., 1.), bbox_transform=ax.transAxes)
         at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
         ax.add_artist(at)
      end
   end

end

"""
    cutplot(data, var, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1,
       level=20)

2D plane cut contourf of 3D box data.
"""
function cutplot(data::Data, var::AbstractString, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf],
   dir="x", sequence=1, level=20)

   x, w = data.x, data.w
   varIndex_ = findindex(data, var)

   X = @view x[:,:,:,1]
   Y = @view x[:,:,:,2]
   Z = @view x[:,:,:,3]

   W = w[:,:,:,varIndex_]

   if dir == "x"
      cut1 = @view X[sequence,:,:]
      cut2 = @view Y[sequence,:,:]
      W    = @view W[sequence,:,:]
   elseif dir == "y"
      cut1 = @view X[:,sequence,:]
      cut2 = @view Z[:,sequence,:]
      W    = @view W[:,sequence,:]
   elseif dir == "z"
      cut1 = @view X[:,:,sequence]
      cut2 = @view Y[:,:,sequence]
      W    = @view W[:,:,sequence]
   end

   if !all(isinf.(plotrange))
      cut1, cut2, W = subsurface(cut1, cut2, W, plotrange)
   end
   if isnothing(ax) ax = plt.gca() end
   c = ax.contourf(cut1, cut2, W, level)

   title(data.head.wnames[varIndex_])

   if dir == "x"
      xlabel("y"); ylabel("z")
   elseif dir == "y"
      xlabel("x"); ylabel("z")
   elseif dir == "z"
      xlabel("x"); ylabel("y")
   end
   c
end

"""
    streamslice(data::Data, var, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf], dir="x",
       sequence=1; kwargs...)

Plot streamlines on 2D slices of 3D box data. Variable names in `var` string must be
separated with `;`.
"""
function streamslice(data::Data, var::AbstractString, ax=nothing;
   plotrange=[-Inf,Inf,-Inf,Inf], dir="x", sequence=1, kwargs...)

   x,w = data.x, data.w
   varStream  = split(var, ";")
   VarIndex1_ = findindex(data, varStream[1])
   VarIndex2_ = findindex(data, varStream[2])

   X = @view x[:,:,:,1]
   Y = @view x[:,:,:,2]
   Z = @view x[:,:,:,3]

   v1 = @view w[:,:,:,VarIndex1_]
   v2 = @view w[:,:,:,VarIndex2_]

   if dir == "x"
      cut1 = @view X[sequence,:,:]
      cut2 = @view Y[sequence,:,:]
      v1   = v1[sequence,:,:]
      v2   = v2[sequence,:,:]
   elseif dir == "y"
      cut1 = @view X[:,sequence,:]
      cut2 = @view Z[:,sequence,:]
      v1   = v1[:,sequence,:]
      v2   = v2[:,sequence,:]
   elseif dir == "z"
      cut1 = @view X[:,:,sequence]
      cut2 = @view Y[:,:,sequence]
      v1   = v1[:,:,sequence]
      v2   = v2[:,:,sequence]
   end

   if !all(isinf.(plotrange))
      cut1, cut2, v1, v2 = subsurface(cut1, cut2, v1, v2, plotrange)
   end

   xi = range(cut1[1,1], stop=cut1[end,1],
      step = (cut1[end,1] - cut1[1,1]) / (size(cut1,1) - 1))
   yi = range(cut2[1,1], stop=cut2[1,end],
      step = (cut2[1,end] - cut2[1,1]) / (size(cut2,2) - 1))

   Xi, Yi = meshgrid(xi, yi)
   if isnothing(ax) ax = plt.gca() end
   s = ax.streamplot(Xi, Yi, v1', v2'; kwargs...)

   if dir == "x"
      xlabel("y"); ylabel("z")
   elseif dir == "y"
      xlabel("x"); ylabel("z")
   elseif dir == "z"
      xlabel("x"); ylabel("y")
   end
   s
end

"""
    plot(data, var, ax=nothing; kwargs...)

Wrapper over `plot` in matplotlib.
"""
function PyPlot.plot(data::Data, var::AbstractString, ax=nothing; kwargs...)
   x, w = data.x, data.w
   varIndex_ = findindex(data, var)
   if isnothing(ax) ax = plt.gca() end
   c = ax.plot(x, w[:,varIndex_]; kwargs...)
end

"""
    scatter(data, var, ax=nothing; kwargs...)

Wrapper over `scatter` in matplotlib.
"""
function PyPlot.scatter(data::Data, var::AbstractString, ax=nothing; kwargs...)
   x, w = data.x, data.w
   varIndex_ = findindex(data, var)
   if isnothing(ax) ax = plt.gca() end
   c = ax.scatter(x, w[:,varIndex_]; kwargs...)
end

"""
    contour(data, var, levels=0; ax=nothing, plotrange=[-Inf,Inf,-Inf,Inf],
       plotinterval=0.1, innermask=false, kwargs...)

Wrapper over `contour` in matplotlib.
"""
function PyPlot.contour(data::Data, var::AbstractString, levels=0; ax=nothing,
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)

   Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; innermask)

   if isnothing(ax) ax = plt.gca() end

   if levels != 0
      c = ax.contour(Xi, Yi, Wi, levels; kwargs...)
   else
      c = ax.contour(Xi, Yi, Wi; kwargs...)
   end
end

"""
    contourf(data, var, levels=0; ax=nothing, plotrange=[-Inf,Inf,-Inf,Inf],
       plotinterval=0.1, innermask=false, kwargs...)

Wrapper over `contourf` in matplotlib.
"""
function PyPlot.contourf(data::Data, var::AbstractString, levels=0; ax=nothing,
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)

   Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; innermask)

   if isnothing(ax) ax = plt.gca() end

   if levels != 0
      c = ax.contourf(Xi, Yi, Wi, levels; kwargs...)
   else
      c = ax.contourf(Xi, Yi, Wi; kwargs...)
   end
end

"""
    tricontourf(data, var, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1,
       kwargs...)

Wrapper over `tricontourf` in matplotlib.
"""
function PyPlot.tricontourf(data::Data, var::AbstractString, ax=nothing;
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...)

   x, w = data.x, data.w
   varIndex_ = findindex(data, var)

   X = vec(x[:,:,1])
   Y = vec(x[:,:,2])
   W = vec(w[:,:,varIndex_])

   # This needs to be modified!!!
   if !all(isinf.(plotrange))
      xyIndex = X .> plotrange[1] .& X .< plotrange[2] .&
         Y .> plotrange[3] .& Y .< plotrange[4]
      X = X[xyIndex]
      Y = Y[xyIndex]
      W = W[xyIndex]
   end
   if isnothing(ax) ax = plt.gca() end
   ax.tricontourf(X, Y, W)
end

"""
    plot_trisurf(data::Data, var::String, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf],
       kwargs...)

Wrapper over `plot_trisurf` in matplotlib.
"""
function PyPlot.plot_trisurf(data::Data, var::AbstractString, ax=nothing;
   plotrange=[-Inf,Inf,-Inf,Inf], kwargs...)

   x, w = data.x, data.w
   varIndex_ = findindex(data, var)

   X = vec(x[:,:,1])
   Y = vec(x[:,:,2])
   W = vec(w[:,:,varIndex_])

   # This needs to be modified!!!
   if !all(isinf.(plotrange))
      xyIndex = X .> plotrange[1] .& X .< plotrange[2] .&
         Y .> plotrange[3] .& Y .< plotrange[4]
      X = X[xyIndex]
      Y = Y[xyIndex]
      W = W[xyIndex]
   end
   if isnothing(ax) ax = plt.gca() end
   ax.plot_trisurf(X, Y, W)
end

"""
    plot_surface(data, var; plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1,
       innermask=false, kwargs...)

Wrapper over `plot_surface` in matplotlib.
"""
function PyPlot.plot_surface(data::Data, var::AbstractString;
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, innermask=false, kwargs...)

   Xi, Yi, Wi = getdata(data, var, plotrange, plotinterval; griddim=2, innermask)
   plot_surface(Xi, Yi, Wi; kwargs...)
end

"""
    streamplot(data, var, ax=nothing; plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1,
       kwargs...)

Wrapper over `streamplot` in matplotlib. `streamplot` does not have **kwargs in the API, but
it supports `density`, `color`, and some other keywords.
"""
function PyPlot.streamplot(data::Data, var::AbstractString, ax=nothing;
   plotrange=[-Inf,Inf,-Inf,Inf], plotinterval=0.1, kwargs...)

   x, w = data.x, data.w
   VarStream  = split(var,";")
   wnames = lowercase.(data.head.wnames)
   VarIndex1_ = findfirst(x->x==lowercase(VarStream[1]), wnames)
   VarIndex2_ = findfirst(x->x==lowercase(VarStream[2]), wnames)

   if data.head.gencoord # generalized coordinates
      X, Y = vec(x[:,:,1]), vec(x[:,:,2])
      if any(isinf.(plotrange))
         if plotrange[1] == -Inf plotrange[1] = minimum(X) end
	      if plotrange[2] ==  Inf plotrange[2] = maximum(X) end
         if plotrange[3] == -Inf plotrange[3] = minimum(Y) end
	      if plotrange[4] ==  Inf plotrange[4] = maximum(Y) end
      end

      # Create grid values first.
      xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
      yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)

      # Is there a triangulation method in Julia?
      tr = matplotlib.tri.Triangulation(X, Y)
      Xi, Yi = meshgrid(xi, yi)

      interpolator = matplotlib.tri.LinearTriInterpolator(tr, w[:,1,VarIndex1_])
      v1 = interpolator(Xi, Yi)

      interpolator = matplotlib.tri.LinearTriInterpolator(tr, w[:,1,VarIndex2_])
      v2 = interpolator(Xi, Yi)

   else # Cartesian coordinates
      if all(isinf.(plotrange))
         Xi, Yi = x[:,:,1]', x[:,:,2]'
         v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'
      else
         xrange = @view x[:,1,1]
         yrange = @view x[1,:,2]
	      if plotrange[1] == -Inf plotrange[1] = minimum(xrange) end
	      if plotrange[2] ==  Inf plotrange[2] = maximum(xrange) end
         if plotrange[3] == -Inf plotrange[3] = minimum(yrange) end
         if plotrange[4] ==  Inf plotrange[4] = maximum(yrange) end

	      w1, w2 = w[:,:,VarIndex1_], w[:,:,VarIndex2_]

         xi = range(plotrange[1], stop=plotrange[2], step=plotinterval)
         yi = range(plotrange[3], stop=plotrange[4], step=plotinterval)

         Xi, Yi = meshgrid(xi, yi)

         spline = Spline2D(xrange, yrange, w1)
	      v1 = spline(Xi[:], Yi[:])
	      v1 = reshape(v1, size(Xi))

	      spline = Spline2D(xrange, yrange, w2)
	      v2 = spline(Xi[:], Yi[:])
	      v2 = reshape(v2, size(Xi))
      end
   end
   if isnothing(ax) ax = plt.gca() end
   ax.streamplot(Xi, Yi, v1, v2; kwargs...)
end

"""
    quiver(data, var, ax=nothing; stride=10, kwargs...)

Wrapper over `quiver` in matplotlib. Only supports Cartesian grid for now.
"""
function PyPlot.quiver(data::Data, var::AbstractString, ax=nothing;
   stride::Integer=10, kwargs...)
   x, w = data.x, data.w
   VarQuiver  = split(var, ";")
   VarIndex1_ = findindex(data, VarQuiver[1])
   VarIndex2_ = findindex(data, VarQuiver[2])

   @views X, Y = x[:,1,1], x[1,:,2]
   @views v1, v2 = w[:,:,VarIndex1_]', w[:,:,VarIndex2_]'

   @views Xq, Yq = X[1:stride:end], Y[1:stride:end]
   v1q = @view v1[1:stride:end, 1:stride:end]
   v2q = @view v2[1:stride:end, 1:stride:end]

   if isnothing(ax) ax = plt.gca() end
   ax.quiver(Xq, Yq, v1q, v2q; kwargs...)
end

"Set colorbar norm and ticks."
function set_colorbar(colorscale, vmin, vmax, data=[1.0])
   if colorscale == :linear || any(<(0), data)
      colorscale == :log && @warn "Nonpositive data detected: use linear scale instead!"
      v1 = isinf(vmin) ? minimum(x->isnan(x) ? +Inf : x, data) : vmin
      v2 = isinf(vmax) ? maximum(x->isnan(x) ? -Inf : x, data) : vmax
      nticks = 7
      levels = matplotlib.ticker.MaxNLocator(nbins=255).tick_values(v1, v2)
      cnorm = matplotlib.colors.BoundaryNorm(levels, ncolors=256, clip=true)
      ticks = range(v1, v2, length=nticks)
   else # logarithmic
      datapositive = data[data .> 0.0]
      v1 = isinf(vmin) ? minimum(datapositive) : vmin
      v2 = isinf(vmax) ? maximum(x->isnan(x) ? -Inf : x, data) : vmax

      cnorm = matplotlib.colors.LogNorm(vmin=v1, vmax=v2)
      ticks = matplotlib.ticker.LogLocator(base=10,subs=collect(0:9))
   end

   cnorm, ticks
end


text:
__precompile__()

module Led
	include("./Types.jl")

	using PyCall
	using .Types

	struct LedConfig
		uc::PyObject
		sectors::Vector{Int}
		ledCount::Int
		ledOffsets::Vector{Int}
		colorSimilarityThreshold::Int
	end

	function init_serial(port, sectors, colorSimilarityThreshold = 5)
		uc = pyimport("serial").Serial(port, 115200, timeout=0.2)
		ledOffsets = Array{Int}(undef, 4)

		for i in 1:4
			ledOffsets[i] = sum(sectors[2:i])
		end

		return LedConfig(
			uc,
			sectors,
			sum(sectors),
			ledOffsets,
			colorSimilarityThreshold
		)
	end

	function clear_frame(config::LedConfig)::Vector{RGB}
		for i in 1:config.ledCount
			config.uc.write([UInt8(i), 0xFF, 0xFF, 0xFF])
		end

		config.uc.write([0xFF])
		return fill((0xFF, 0xFF, 0xFF), config.ledCount)
	end

	function send_frame(config::LedConfig, data::Vector{RGB}, previousFrame=Vector{RGB}())::Vector{RGB}
		for i in 1:size(data, 1)
			p = data[i]

			if (length(previousFrame) >= i)
				if (isColorSimilar(p, previousFrame[i], config.colorSimilarityThreshold))
					continue
				end
			end

			config.uc.write(UInt8.([i - 1, p[1], p[2], p[3]]))
		end

		config.uc.write([0xFF])
		return data
	end

	function isColorSimilar(a::RGB, b::RGB, threshold = 5)
		return (
		   abs(convert(Int32, a[1]) - b[1]) <= threshold
		&& abs(convert(Int32, a[2]) - b[2]) <= threshold
		&& abs(convert(Int32, a[3]) - b[3]) <= threshold
		)
	end

	export init_serial, clear_frame, send_frame, LedConfig
end


text:
module HomogeneousCoordinates

	using Interpolations
	using StaticArrays

	export transform

	"""
		Not sure if this is true for all homogeneous coordinates, but this requires
		the matrix [4,4] element to be equal to one.
	"""
	function transform(a::AbstractArray{<: Real, 3}, b_size::NTuple{N, Int64}, Ma::AbstractMatrix{<: Real}, Mb::AbstractMatrix{<: Real}) where N
		# TODO generated for any dimension

		(rows, columns, slices) = size(a)

		b = zeros(b_size)
		a_itp = extrapolate(
			interpolate(
				a,
				BSpline(Linear())
			),
			0
		)

		Mba = (inv(Ma) * Mb)[1:3, :]
		for i in CartesianIndices(b_size)
			b_vec = SVector{4, Float64}(i[1], i[2], i[3], 1)
			a_vec = Mba * b_vec
			#a_vec = a_vec ./ a_vec[4]
			b[i] = a_itp(a_vec...)
		end
		return b
	end

	# This is not yet useful, it can be useful if you only want to change the coordinate system,
	# but not into the system of another scan. In that case the target array size is not defined,
	# and it needs to be determined by transforming the corners. Like that, no empty space is stored.
	# Such a function must then also return the actual indices of the zero index point of the resulting array.
	#b_size = let
	#	corners = @SMatrix [
	#		0 rows	0		rows	0		rows	0		rows;
	#		0 0		columns columns	0		0		columns	columns;
	#		0 0		0		0		slices	slices	slices	slices;
	#		1 1		1		1		1		1		1		1
	#	]
	#	transformed_corners = (inv(Mb) * Ma)[1:3, :] * corners
	#	b_size = map(
	#		limits -> ceil(Int, limits[2] - limits[1]),
	#		dropdims(extrema(transformed_corners; dims=2); dims=2)
	#	)
	#end

end


text:
"""
assuming that you load PowerDynamics beforehand
"""

## patching PD.jl with ControlledPowerGrid

using LightGraphs: AbstractGraph

# stack dynamics with higher level controller
struct ControlledPowerGrid
    graph:: G where G <: AbstractGraph
    nodes
    lines
    controller # oop function
end

# Constructor
function ControlledPowerGrid(control, pg::PowerGrid, args...)
    ControlledPowerGrid(pg.graph, pg.nodes, pg.lines, (u, p, t) -> control(u, p, t, args...))
end

import PowerDynamics: PowerGrid

PowerGrid(pg::PowerGrid) = pg
PowerGrid(pg::ControlledPowerGrid) = PowerGrid(pg.graph, pg.nodes, pg.lines)

import PowerDynamics: State, PowerGridSolution

State(pg::ControlledPowerGrid, vec) = State(PowerGrid(pg), vec)
PowerGridSolution(vec, pg::ControlledPowerGrid) = PowerGridSolution(vec, PowerGrid(pg))

import PowerDynamics: systemsize

systemsize(cpg::ControlledPowerGrid) = systemsize(PowerGrid(cpg.graph, cpg.nodes, cpg.lines))

## define ControlledPowerGrid rhs

using OrdinaryDiffEq: ODEFunction
import PowerDynamics: rhs, get_current

function rhs(cpg::ControlledPowerGrid)
    open_loop_dyn = PowerGrid(cpg.graph, cpg.nodes, cpg.lines) |> rhs
    function cpg_rhs!(du, u, p, t)
        # Get controller state
        control_output = cpg.controller(u, p, t)
        # Calculate the derivatives, passing the controller state through
        open_loop_dyn(du, u, control_output, t)
    end
    ODEFunction{true}(cpg_rhs!, mass_matrix=open_loop_dyn.mass_matrix, syms=open_loop_dyn.syms)
end

using NetworkDynamics: network_dynamics, StaticEdgeFunction

function rhs(pg::PowerGrid)
    network_dynamics(map(construct_vertex, pg.nodes), map(construct_edge, pg.lines), pg.graph; parallel=true)
end

get_current(state::State, n) = begin
    vertices = map(construct_vertex, state.grid.nodes)
    edges = map(construct_edge, state.grid.lines)
    sef = StaticEdgeFunction(vertices,edges,state.grid.graph, parallel=true)
    (e_s, e_d) = sef(state.vec, Nothing, 0)
    total_current(e_s[n], e_d[n])
end

## operation point

import PowerDynamics: find_operationpoint
using NLsolve: nlsolve, converged
using ForwardDiff: jacobian
using LinearAlgebra: eigvals
using PowerDynamics: AbstractNode
using OrdinaryDiffEq: ODEFunction
using DifferentialEquations: DynamicSS
using NLsolve: nlsolve, converged
using LinearAlgebra: pinv

struct RootRhs_ic
    rhs
    mpm
end
function (rr::RootRhs_ic)(x)
    dx = similar(x)
    rr.rhs(dx, x, nothing, 0.)
    rr.mpm * dx .- dx
end

function RootRhs_ic(of::ODEFunction)
    mm = of.mass_matrix
    @assert mm != nothing
    mpm = pinv(mm) * mm
    RootRhs_ic(of.f, mpm)
end

"""
Makes an type specific initial guess to help the operation point search.
The voltage is of all nodes is fixed to the voltage of the first SlackAlgebraic
in the system. The other symbols are set to zero.
Inputs:
    pg: Power grid, a graph containg nodes and lines
Outputs:
    guess: Type specific initial guess
"""
function initial_guess(pg)
    if SlackAlgebraic ∉ pg.nodes .|> typeof
        @warn "There is no slack bus in the system to balance powers."
    end

    sl = findfirst(SlackAlgebraic  ∈  pg.nodes .|> typeof)
    slack = pg.nodes[sl]

    type_guesses = guess.(pg.nodes, slack.U)
    return vcat(type_guesses...)
end

function initial_guess(pg, v_vector)
    type_guesses = guess.(pg.nodes, v_vector)
    return vcat(type_guesses...)
end

"""
guess(::Type{SlackAlgebraic}) = [slack.U, 0.]         #[:u_r, :u_i]
guess(::Type{PQAlgebraic}) = [slack.U, 0.]            #[:u_r, :u_i]
guess(::Type{ThirdOrderEq}) = [slack.U, 0., 0., 0.]   #[:u_r, :u_i, :θ, :ω]
"""
function guess(n::AbstractNode, slack_voltage)
    voltage = zeros(dimension(n))
    voltage[1] = real(slack_voltage)  #[:u_r, :u_i]
    voltage[2] = imag(slack_voltage)
    return voltage
end

function guess(n::SlackAlgebraic, slack_voltage)
    @assert n.U ≈ slack_voltage
    return [real(n.U), imag(n.U)]
end

function find_steady_state(pg, op_guess=nothing)
    if op_guess isa Nothing
        op_guess = initial_guess(pg)
    end
    ode = rhs(pg)
    op_prob = ODEProblem(ode, op_guess, (0.0, 1.0))
    sol = solve(
        SteadyStateProblem(op_prob),
        DynamicSS(Rodas5(); abstol = 1e-8, reltol = 1e-6, tspan = Inf),
    )
    return State(pg, sol.u)
end

function find_operationpoint(pg::ControlledPowerGrid, ic_guess = nothing)
    if SlackAlgebraic ∉ pg.nodes .|> typeof
        @warn "There is no slack bus in the system to balance powers. Currently not making any checks concerning assumptions of whether its possible to find the fixed point"
    end
    if SwingEq ∈ pg.nodes .|> typeof
        throw(OperationPointError("found SwingEq node but these should be SwingEqLVS (just SwingEq is not yet supported for operation point search)"))
    end

    if ic_guess === nothing
        system_size = systemsize(pg)
        ic_guess = ones(system_size)
    end

    rr = RootRhs_ic(rhs(pg))
    nl_res = nlsolve(rr, ic_guess)
    if converged(nl_res) == true
        return State(pg, nl_res.zero)
    else
        throw(OperationPointError("Failed to find initial conditions on the constraint manifold!"))
    end
end

function check_operationpoint(pg, fixed_point::State; p=nothing, t0=0.)
    ode = rhs(pg)
    M = ode.mass_matrix
    f!(dx, x) = ode(dx, x, p, t0)
    j(x) = (dx = similar(x); jacobian(f!, dx, x))
    λs = eigvals(Array(M), j(fixed_point.vec)) .|> real
    λ = λs[isfinite.(λs)] |> extrema
    println("Jacobian spectrum \nmin : ", first(λ), "\nmax : ",last(λ), "\nstable : ", last(λ)==0)
    return λ
end
check_operationpoint(op::State) = check_operationpoint(op.grid, op.vec)

## ADN global control

function ADN(pg::PowerGrid, DGUnit_Type, P_ref, Q_ref)
    slack_idx = findfirst([node isa SlackAlgebraic for node in pg.nodes])

    Y = PiModel(first(pg.lines))
    sl_func = u -> conj(- Y[1,1] - Y[1, 2]*complex(u[3:4]...))

    args = (pg, slack_idx, P_ref, Q_ref, sl_func)

    ControlledPowerGrid(
        ADN_control,
        pg,
        args...
    )
end

function ADN_control(u, p, t, pg, slack_idx, P_ref, Q_ref, sl_func)

    # this works but raises annoying ND.jl warning every time
    #obs = State(pg, u) # TODO probably better to use graphdata from ND?
    #S_slack = obs[slack_idx, :s]

    # take the high-voltage side as the measurement reference
    # S will be positive since node currents are positive by convention

    S_slack = sl_func(u)

    ΔP = P_ref(t) - real(S_slack)
    ΔQ = Q_ref(t) - imag(S_slack)

    # return ND parameter tuple
    # if p isa Nothing
    #     p_cont = ((ΔP, ΔQ), nothing)
    # else
    #     node_pars, edge_pars = p
    #
    #     p_cont = (
    #         [(ΔP, ΔQ, par...) for par in node_pars],
    #         edge_pars,
    #     )
    # end
    return ((ΔP, ΔQ), nothing)
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
module Precompile1

using KeywordCalls

struct Foo{N,T}
    nt::NamedTuple{N,T}
end

@kwcall Foo(a,b)

Foo(nt::NamedTuple{(:a,:b), Tuple{A,B}}) where {A,B} = Foo{(:a,:b), Tuple{A,B}}(nt)

export Foo

end # module


text:
"""
    ConstantGainEKSPF()

Constant gain from [1], Algorithm 1.

### Reference

[1] Venugopal, M., Vasu, R. M., & Roy, D. (2016). An Ensemble Kushner-Stratonovich-Poisson Filter for Recursive Estimation in Nonlinear Dynamical Systems. IEEE Transactions on Automatic Control, 61(3), 823–828. https://doi.org/10.1109/TAC.2015.2450113
"""
struct ConstantGainEKSPF <: GainEstimationMethod end

function solve!(eq::PoissonEquation, method::ConstantGainEKSPF)
    n, N, m = size(eq.gain)
    htilde = Htilde(eq)
    @inbounds for j in 1:m
        for i in 1:n
            g = 0.
            for l in 1:N
                g += eq.positions[i, l] * htilde[j, l]
            end
            for k in 1:N
                eq.gain[i, k, j] = g / (N * eq.mean_H[j, 1])
            end
        end
    end
    return eq.gain
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using Documenter, Twiddle

makedocs(
    format = Documenter.HTML(),
    sitename = "Twiddle.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Counting bits" => "bitcounting.md",
            "Reference Manual" => "reference.md",
            "Working with Nibbles" => "nibbles.md",
            "Glossary" => "glossary.md"
        ]
    ],
    authors = "Ben J. Ward & Luis Yanes."
)
deploydocs(
    repo = "github.com/BenJWard/Twiddle.jl.git",
    deps = nothing,
    make = nothing,
    push_preview = true
)


text:
module MyExample

greet() = print("Hello World!")

include("extra_file.jl")
export my_f

end # module


text:
using Convex

function readmylines(rd)
    readline(rd) #Initialize time
    readline(rd) #Dashes
    l1 = readline(rd)
    readline(rd) #Dashes
    l2 = readline(rd)[1:7]
    l3 = readline(rd)[1:7]
    l4 = readline(rd)
    return l1, l2, l3, l4
end

#Expected outputs
o11 = " Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | cg  | time"
o12 = " Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time"
o2  = "   100|"
o3  = "   200|"
o4  = "Found solution i=200"

Random.seed!(10)
n = 500
A = sprandn(n, 2n, 0.1)
x̄ = randn(2n)
b = A*x̄

x = Variable(2n)
p = minimize(norm(A*x-b), sum(x) == sum(x̄))

origstdout = stdout

#With cg output
s = GAPA(0.8, 0.9, direct=false, verbose=2, debug=0, eps=1e-8, checki=100)

rd,wr = redirect_stdout()
solve!(p,s)
l1, l2, l3, l4 = readmylines(rd)
redirect_stdout(origstdout)

@test l1 == o11
@test l2 == o2
@test l3 == o3
@test l4 == o4

@test evaluate(sum(x)-sum(x̄))       ≈ 0.0 atol=1e-8
@test evaluate(maximum(abs(A*x-b))) ≈ 0.0 atol=1e-8

#Without cg output
s = GAPA(0.8, 0.9, direct=true, verbose=2, debug=0, eps=1e-8, checki=100)

rd,wr = redirect_stdout()
solve!(p,s)
l1, l2, l3, l4 = readmylines(rd)
redirect_stdout(origstdout)

@test l1 == o12
@test l2 == o2
@test l3 == o3
@test l4 == o4

@test evaluate(sum(x)-sum(x̄))       ≈ 0.0 atol=1e-8
@test evaluate(maximum(abs(A*x-b))) ≈ 0.0 atol=1e-8


text:
export MachOSection, MachOSections, MachOSectionRef

"""
    MachOSection

Mach-O `Section` type, containing information about the section name, segment
name, size, address, etc... of the `Section.`
"""
abstract type MachOSection{H <: MachOHandle} <: Section{H} end
@io struct MachOSection32{H <: MachOHandle} <: MachOSection{H}
    sectname::fixed_string{UInt128}
    segname::fixed_string{UInt128}
    addr::UInt32
    size::UInt32
    offset::UInt32
    align::UInt32
    reloff::UInt32
    nreloc::UInt32
    flags::UInt32
    reserved1::UInt32
    reserved2::UInt32
end

@io struct MachOSection64{H <: MachOHandle} <: MachOSection{H}
    sectname::fixed_string{UInt128}
    segname::fixed_string{UInt128}
    addr::UInt64
    size::UInt64
    offset::UInt32
    align::UInt32
    reloff::UInt32
    nreloc::UInt32
    flags::UInt32
    reserved1::UInt32
    reserved2::UInt32
end

function read(oh::MachOHandle, ::Type{ST}) where {ST <: MachOSection}
    type_func = section_header_type
    size_func = section_header_size
    return read_struct(oh, type_func, size_func, "Section Header")
end

section_name(s::MachOSection) = "$(s.segname),$(s.sectname)"
section_size(s::MachOSection) = s.size
section_offset(s::MachOSection) = s.offset
section_address(s::MachOSection) = s.addr
is64bit(s::MachOSection) = isa(s, MachOSection64)

"""
    MachOSections

Mach-O `Section` object collection type.  Mach-O sections are split out over
multiple segments, one per load command.  As an example, most executables have
at least two segments, `__DATA` and `__TEXT`, each of which have multiple
sections embedded within them.  A `MachOSections` object is created from
multiple segments (in this case, realized as multiple load commands) which
contain the necessary sections, and these sections can then be accessed as
desired.
"""
struct MachOSections{H <: MachOHandle, S <: MachOSection} <: Sections{H}
    handle::H
    sections::Vector{S}
end

# Special section_header_type() overrides for segments:
section_header_type(::MachOSegment32Cmd) = MachOSection32
section_header_type(::MachOSegment64Cmd) = MachOSection64
section_header_type(seg::MachOSegment) = section_header_type(MachOLoadCmd(seg))
@derefmethod section_header_type(seg::MachOSegmentRef)

function MachOSections(segs::MachOSegments)
    oh = handle(segs)

    # Now, load the sections off of each segment:
    sections = Any[]
    for seg in segs
        # Seek to the beginning of the section headers for this segment
        seg_size = packed_sizeof(typeof(MachOLoadCmd(deref(seg))))
        sections_pos = position(MachOLoadCmdRef(seg)) + seg_size
        seek(oh, sections_pos)

        # For each section header, read it in:
        for sect_idx in 1:segment_num_sections(seg)
            sect = read(oh, section_header_type(seg))
            push!(sections, sect)
        end
    end

    # Re-run type inference on this array, so it packs nicely
    sections = [s for s in sections]

    # Bake'em away, Toys!
    return MachOSections(oh, sections)
end

handle(sections::MachOSections) = sections.handle
endof(sections::MachOSections) = endof(sections.sections)
getindex(sections::MachOSections, idx) = getindex(sections.sections, idx)

Sections(segs::MachOSegments) = MachOSections(segs)
Sections(oh::MachOHandle) = Sections(Segments(oh))

"""
    MachOSectionRef

Mach-O `SectionRef` type
"""
struct MachOSectionRef{H <: MachOHandle} <: Sections{H}
    sections::MachOSections{H}
    section::MachOSection{H}
    idx::UInt32
end

function SectionRef(sections::MachOSections, s::MachOSection, idx)
    return MachOSectionRef(sections, s, UInt32(idx))
end
@derefmethod is64bit(s::MachOSectionRef)


text:
name = "Ryan"

language = "Julia"

message = "Hi, my name is $name and I am practicing the $language language."

println(message)

another_language = "PHP"

println("$language uses the same string interpolation syntax as $another_language")

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
using Printf

const SRC_PATH = "../../src"

include(joinpath(SRC_PATH, "NEAT/config.jl"))
include(joinpath(SRC_PATH, "NEAT/global.jl"))
include(joinpath(SRC_PATH, "NEAT/genome.jl"))
include(joinpath(SRC_PATH, "NEAT/chromosome.jl"))
include(joinpath(SRC_PATH, "MutualInformation/information-plane.jl"))
include(joinpath(SRC_PATH, "NEAT/fitness.jl"))
include(joinpath(SRC_PATH, "NEAT/species.jl"))
include(joinpath(SRC_PATH, "NEAT/population.jl"))
include(joinpath(SRC_PATH, "NEAT/simulation_state.jl"))
include(joinpath(SRC_PATH, "NEAT/simulation-plot.jl"))
include(joinpath(SRC_PATH, "NEAT/checkpoint.jl"))
include(joinpath(SRC_PATH, "NEAT/neat.jl"))
#include(joinpath(SRC_PATH, "MLPClassifier/MLPClassifier.jl"))
include(joinpath(SRC_PATH, "FileInterface/fileInterface.jl"))

const DEFAULT_DATASET_PATH = joinpath("..", SRC_PATH, "datasets")

function run_NEAT_epochs(config_file::String, generations::Int=1, report::Bool=false)
    params =  loadConfig(config_file)
    cf = Config(params)
    print("checkpoint filename: ", cf.checkpoint_filename)
    filename = string(cf.checkpoint_filename, ".jld2")
    if isfile(filename)
        g = loadFromFile(filename, "Global")
    else
        g = Global(cf)
    end
    state = SimulationState(g)
    epoch(state, generations, report)
    return state
end

function print_best_solution_perfomance(state::SimulationState)
    best = state.population.best_fitness[end]
    dt_setup = state.setup.dataset_setup
    println("Best solution's fitness: ", best.fitness)
    mlp = convert_to_FluxNet(best, dt_setup.labels)
    error = error_rate(zip(dt_setup.input_train, dt_setup.output_train), mlp)
    acc = accuracy(zip(dt_setup.input_train, dt_setup.output_train), mlp)
    println("Accuracy after NEAT: ", acc)
    println("Error rate after NEAT: ", error)
end

#=
params0 =  loadConfig("test_demo_config.txt")
cf0 = Config(params0)
g0 = Global(cf0)
=#
#ch = create_minimal_chromossome(g)
#mutate_add_layer!(ch, g)
#mutate_add_node!(ch, g)
#mutate_remove_node!(ch, g)
#mutate_weights!(ch, g)
#mutate_enable_bit!(ch)
#evaluation = Evaluation(0,2)
#evaluation(ch)

#population0 = Population(g0)
#epoch(g0, population0, 10)
#time: 0.2 secs

#=
println("--> Dataset demo 2 (overlap)")
p1 = run_NEAT_epochs("test_overlap_config.txt", 1000, false)
#plot_neat_fitness(p1, imagename="Figures/fitness-demo2-t.png")
#plot_neat_IP(p1, imagename="Figures/IP-demo2-t.png")
#plot_neat_species(p1, imagename="Figures/species-demo2-t.png")
print_best_solution_perfomance(p1)
=#
#=
Time: 86.14915895462036 seconds.
Best solution's fitness: 0.42653352353780294
Accuracy after NEAT: 0.598
Error rate after NEAT: 0.402
=#

#=
println("--> Dataset 1")
p2 = run_NEAT_epochs("test_dataset1_config.txt", 1000, false)
#plot_neat_fitness(p2, imagename="Figures/fitness-dt1.png")
#plot_neat_IP(p2, imagename="Figures/IP-dt1.png")
#plot_neat_species(p2, imagename="Figures/species-dt1.png")
print_best_solution_perfomance(p2)
#=
Time: 1550.5301270484924 seconds.
Best solution's fitness: 0.4184397163120568
Accuracy after NEAT: 0.59
Error rate after NEAT: 0.41
=#

println("--> Dataset 2")
p3 = run_NEAT_epochs("test_dataset2_config.txt", 1000, false)
plot_neat_fitness(p3, imagename="Figures/fitness-dt2.png")
plot_neat_IP(p3, imagename="Figures/IP-dt2.png")
plot_neat_species(p3, imagename="Figures/species-dt2.png")
print_best_solution_perfomance(p3)

println("--> Dataset 3")
p4 = run_NEAT_epochs("test_dataset3_config.txt", 1000, false)
plot_neat_fitness(p4, imagename="Figures/fitness-dt3.png")
plot_neat_IP(p4, imagename="Figures/IP-dt3.png")
plot_neat_species(p4, imagename="Figures/species-dt3.png")
print_best_solution_perfomance(p4)

println("--> Dataset 4")
p5 = run_NEAT_epochs("test_dataset4_config.txt", 1000, false)
plot_neat_fitness(p5, imagename="Figures/fitness-dt4.png")
plot_neat_IP(p5, imagename="Figures/IP-dt4.png")
plot_neat_species(p5, imagename="Figures/species-dt4.png")
print_best_solution_perfomance(p5)
=#

#=
println("--> Dataset 5")
p6 = run_NEAT_epochs("test_dataset5_config.txt", 1000, false)
plot_neat_fitness(p6, imagename="Figures/fitness-dt5.png")
plot_neat_IP(p6, imagename="Figures/IP-dt5.png")
plot_neat_species(p6, imagename="Figures/species-dt5.png")
print_best_solution_perfomance(p6)
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
##############################################################
## Test error propagation for the different error functions ##
##############################################################
@testset "error propagation" begin

    ε = 100*eps(Float64)
    @polyvar x y z

    ###############
    # Translation #
    ###############
    # Coefficients
    Cspher = SphericalHarmonicCoefficients(zeros(25));
    Cspher[3,-2] = 1/sqrt(15)*sqrt((4*pi)/(2*3+1));
    Csolid = SphericalHarmonicCoefficients(zeros(25),1.0,true);
    Csolid[3,-2] = 1/sqrt(15);

    # Translation
    v = [1,-1,2]

    coeffsTransSpher = SphericalHarmonicExpansions.errorTranslation(Cspher,v)
    coeffsTransSolid = SphericalHarmonicExpansions.errorTranslation(Csolid,v)

    # Correct solution
    CtransSpher = zeros(25)
    CtransSpher[1:11] = sqrt(4pi/3).*[2*sqrt(3), 2, 1, 2, 2/sqrt(5), 1/sqrt(5), 0, 1/sqrt(5), 0, 0, 1/sqrt(35)]
    CtransSpher = SphericalHarmonicCoefficients(CtransSpher)
    CtransSolid = zeros(25)
    CtransSolid[1:11] = [2, 2, 1, 2, 2/sqrt(3), 1/sqrt(3), 0, 1/sqrt(3), 0, 0, 1/sqrt(15)]
    CtransSolid = SphericalHarmonicCoefficients(CtransSolid,1.0,true)

    # Test
    @test isapprox(CtransSpher,coeffsTransSpher)
    @test isapprox(CtransSolid,coeffsTransSolid)

    ##############
    # Quadrature #
    ##############
    # 3-design:
    coordinates = [1. 0. 0.;
                   -1. 0. 0.;
                   0. 1. 0.;
                   0. -1. 0.;
                   0. 0. 1.;
                   0. 0. -1.]

    # Calculating the error of the coefficients up to l=1
    δ = ones(6)
    εquadr = SphericalHarmonicExpansions.errorSphericalQuadrature(δ,coordinates,1);

    # Correct solution
    εcorr = [sqrt(4*pi); sqrt(4*pi/3); sqrt(4*pi/3); sqrt(4*pi/3)]

    # Test
    @test isapprox(εquadr.c,εcorr,atol=ε)
end


text:
using PackageCompiler

scripts_dir = joinpath(@__DIR__, "..", "scripts")
toml_file = joinpath(scripts_dir, "Project.toml")
script = joinpath(scripts_dir, "exploration.jl")
precompiles_file = joinpath(scripts_dir, "precompiles.jl")
blacklist = ["CentroidalTrajOpt", "QPWalkingControl", "QPControl"]#, "AtlasRobot"]
sysimg_dest = joinpath(scripts_dir, "sysimg.so")

@info "Generating precompile statements"
PackageCompiler.snoop(nothing, toml_file, script, precompiles_file, false, blacklist)
@info "Done generating precompile statements"

sleep(1)
@info "Compiling system image"
new_sysimg, old_sysimg = compile_incremental(toml_file, precompiles_file)
@info "Done compiling system image"
cp(new_sysimg, sysimg_dest, force=true)


text:
module IRTools

export @code_ir

module Inner

  using MacroTools
  using MacroTools: @q, prewalk, postwalk
  import ..IRTools

  export @code_ir

  @nospecialize

  include("reflection/reflection.jl")
  include("utils.jl")

  include("ir/ir.jl")
  include("ir/utils.jl")
  include("ir/wrap.jl")
  include("ir/print.jl")
  include("ir/parse.jl")

  include("reflection/utils.jl")
  include("reflection/dynamo.jl")

  include("passes/passes.jl")
  include("passes/inline.jl")
  include("passes/cps.jl")
  include("passes/relooper.jl")
  include("passes/stackifier.jl")
  include("passes/registers.jl")

  include("interpret.jl")
  include("eval.jl")

end

let exports = :[
      # IR
      IRTools, IR, Block, BasicBlock, Variable, Statement, Branch, Pipe, CFG, Slot, branch, var, stmt, arguments, argtypes,
      branches, undef, unreachable, isreturn, isconditional, block!, deleteblock!, branch!, argument!, return!,
      canbranch, returnvalue, returntype, emptyargs!, deletearg!, block, blocks, successors, predecessors,
      xcall, exprtype, exprline, isexpr, insertafter!, explicitbranch!, prewalk, postwalk,
      prewalk!, postwalk!, finish, substitute!, substitute,
      # Passes/Analysis
      definitions, usages, dominators, domtree, domorder, domorder!, renumber,
      merge_returns!, expand!, prune!, ssa!, inlineable!, log!, pis!, func, evalir,
      Simple, Loop, Multiple, reloop, stackify, functional, cond, WorkQueue,
      Graph, liveness, interference, colouring, inline,
      # Reflection, Dynamo
      Meta, Lambda, meta, dynamo, transform, refresh, recurse!, self,
      varargs!, slots!,
      ].args
  append!(exports, Symbol.(["@code_ir", "@dynamo", "@meta"]))
  for x in exports
    @eval const $x = Inner.$x
  end
  @eval module All
    $([:(import ..IRTools: $x) for x in exports]...)
    export $(exports...)
  end
end

end # module


text:
using Accumulo
using Base.Test

const SEP = "\n" * "="^20 * "\n"
print_header(message) = println(SEP, message, SEP)

function print_table_config(session, tbl)
    cfg = table_config(session, tbl)
    println("Configuration of $tbl:")
    for (n,v) in cfg
        println("\t$n => $v")
    end
end

function print_table_iters(session, tbl)
    tbliters = iters(session, tbl)
    println("Iterators on $tbl: $tbliters")

    for (iname, iscopes) in tbliters
        for iscope in iscopes
            iprops = iter(session, tbl, iname, iscope)
            println("Iterator $iname on $tbl in scope $iscope: $iprops")
        end
    end
end

function print_info(session)
    print_header("SERVER INFO")
    tbls = tables(session)
    println("Tables: $tbls")

    diskusage = table_du(session, tbls...)
    println("Disk updage of $tbls: $diskusage")

    for tbl in tbls
        println("Table: $tbl")
        println("===========================")

        diskusage = table_du(session, tbl)
        println("Disk updage of $tbl: $diskusage")

        print_table_config(session, tbl)
        print_table_iters(session, tbl)

        tblconstr = constraints(session, tbl)
        println("Constraints on $tbl: $tblconstr")
    end
end

function test_table_admin(session)
    print_header("TABLE ADMIN")
    try
        println("creating mytable...")
        @test !table_exists(session, "mytable")
        table_create(session, "mytable")
        @test table_exists(session, "mytable")

        println("renaming mytable to mytable1...")
        table_rename(session, "mytable", "mytable1")
        @test table_exists(session, "mytable1")
        @test !table_exists(session, "mytable")

        println("cloning mytable1 to mytable2...")
        table_clone(session, "mytable1", "mytable2")
        @test table_exists(session, "mytable2")

        println("deleting mytable1...")
        table_delete(session, "mytable1")
        @test !table_exists(session, "mytable1")
        @test table_exists(session, "mytable2")

        println("deleting mytable2...")
        table_delete(session, "mytable2")
        @test !table_exists(session, "mytable2")
    finally
        for tbl in ("mytable", "mytable1", "mytable2")
            println("cleaning up $tbl...")
            table_exists(session, tbl) && table_delete(session, tbl)
        end
    end
end

function test_table_readwrite(session)
    N = 10
    TNAME = "mytable"

    function _scan_verify(valpfx)
        println("creating scanner...")
        scanner(session, TNAME) do scan
            println("\tstarting iteration...")
            nrecs = 0
            for rec in records(scan)
                row = bytestring(rec.key.row)
                colfam = bytestring(rec.key.colFamily)
                colqual = bytestring(rec.key.colQualifier)
                colvis = bytestring(rec.key.colVisibility)
                val = bytestring(rec.value)
                println("\trow($(row)), column($(colfam):$(colqual)), visibility($(colvis)), value($(val))")
                rownum = parse(Int, row[4:end])
                @test colfam == "colf"
                @test colqual == "colq"
                @test colvis == ""
                if isempty(valpfx)
                    @test val == ""
                else
                    @test val == "$valpfx$rownum"
                end
                nrecs += 1
            end
            println("\tfinished iteration ($nrecs records)...")
        end
        println("closed scanner")
    end

    function _create()
        println("creating $TNAME...")
        @test !table_exists(session, TNAME)
        table_create(session, TNAME; versioning=false)
        @test table_exists(session, TNAME)
        table_versions!(session, TNAME, 1)
        #print_table_config(session, TNAME)
        #print_table_iters(session, TNAME)
    end

    function _write()
        println("creating writer...")
        batch_writer(session, TNAME) do writer
            println("\tstarting batch...")
            upd = batch()
            println("\tcollecting mutations...")
            for rownum in 1:N
                update(upd, "row$rownum", "colf", "colq", "", "value$rownum")
            end
            println("\tupdating mutations...")
            update(writer, upd)
            println("\tflushing writer...")
            flush(writer)
        end
        println("closed writer")
    end

    function _update()
        println("creating conditional batch writer...")
        conditional_batch_writer(session, TNAME) do writer
            println("\tstarting batch...")
            upd = conditional_batch()
            println("\tcollecting mutations...")
            for rownum in 1:N
                row = "row$rownum"
                val = "value$rownum"
                update(where(upd, row, "colf", "colq"; value=val), "colf", "colq", "", "newvalue$rownum")
                #update(where(upd, row, "colf", "colq", "public"; value=val), "colf", "colq", "public", "newvalue$rownum")
            end
            println("\tupdating mutations...")
            res = update(writer, upd)
            status_counts = Dict(ConditionalStatus.ACCEPTED => 0, ConditionalStatus.REJECTED => 0, ConditionalStatus.VIOLATED => 0, ConditionalStatus.UNKNOWN => 0, ConditionalStatus.INVISIBLE_VISIBILITY => 0)
            for v in values(res)
                status_counts[v] += 1
            end
            println("\tstatus counts: $status_counts")
            println("\tflushing writer...")
            flush(writer)
        end
        println("closed writer")
    end

    function _delete()
        println("creating writer for deletes...")
        batch_writer(session, TNAME) do writer
            println("\tstarting batch...")
            upd = batch()
            println("\tcollecting mutations...")
            for rownum in 1:N
                delete(upd, "row$rownum", "colf", "colq")
            end
            println("\tupdating mutations...")
            update(writer, upd)
            println("\tflushing writer...")
            flush(writer)
        end
        println("closed writer")
    end

    print_header("TABLE READ WRITE")
    try
        _create()
        _write()
        _scan_verify("value")
        _update()
        _scan_verify("newvalue")
        _delete()
        _scan_verify("")
    finally
        println("cleaning up...")
        table_exists(session, TNAME) && table_delete(session, TNAME)
    end
end

session = AccumuloSession("127.0.1.1", 42424, AccumuloAuthSASLPlain("root", "root"))
print_info(session)
test_table_admin(session)
test_table_readwrite(session)
close(session)


text:
"""
    EmptyChain

It's a dummy chain class that does not store anything.
"""
struct EmptyChain <: AbstractChain{Nothing,Nothing} end
length(ch::EmptyChain) = 0
dim(ch::EmptyChain) = 0
keys(ch::EmptyChain) = Nothing[]
values(ch::EmptyChain) = Nothing[]
getindex(ch::EmptyChain, i) = nothing

"""
    Chain{IX<:Integer, R}

The chain of integer elements for `R`-module.

This chain implementation uses a dictionary for storing elements and coefficients.
"""
mutable struct Chain{IX<:Integer, R} <: AbstractChain{IX,R}
    dim::Int
    cells::Dict{IX,R}
end

Chain(d::Int, elems::AbstractVector, coefs::AbstractVector) = Chain(d, Dict(zip(elems, coefs)))
Chain(elems::AbstractVector, coefs::AbstractVector) = Chain(0, elems, coefs)
Chain(d::Int, ::Type{IX}, ::Type{R}) where {R, IX<:Integer} = Chain(d, IX[], R[])
Chain(::Type{IX}, ::Type{R}) where {R, IX<:Integer} = Chain(0, IX, R)
Chain(d::Int, ::Type{R}) where {R} = Chain(d, UInt, R)
Chain(::Type{R}) where {R} = Chain(0, UInt, R)

# implement interface
dim(ch::Chain) = ch.dim
length(ch::Chain) = length(ch.cells)
copy(ch::Chain{IX,R}) where {R, IX<:Integer} = Chain{IX,R}(ch.dim, copy(ch.cells))
keys(ch::Chain{IX,R}) where {R, IX<:Integer} = keys(ch.cells)
values(ch::Chain{IX,R}) where {R, IX<:Integer} = values(ch.cells)
getindex(ch::Chain{IX,R}, k::IX) where {R, IX<:Integer} = get(ch.cells, k, zero(R))
setindex!(ch::Chain{IX,R}, c::R, k::IX) where {R, IX<:Integer} = setindex!(ch.cells, c, k)

function push!(ch::Chain{IX,R}, e::Pair{IX,R}) where {R, IX<:Integer}
    (k,v) = e
    if k ∈ ch
        ch.cells[k] += v
    else
        push!(ch.cells, e)
    end
    ch
end

function iterate(ch::Chain{IX,R}, state...) where {R, IX<:Integer}
    y = iterate(ch.cells, state...)
    y === nothing && return nothing
    return (y[1], y[2])
end

function simplify(ch::Chain)
    for (k,v) in ch.cells
        v == 0 && delete!(ch.cells, k)
    end
    return ch
end


text:
"""
    changeSymbol!(expr, namedTuple)

Changes the symbols of the named tuple found in the expression into their actual
values.
"""
function changeSymbol!(expr, namedTuple)
    function changeRec!(list)
        for i=1:length(list)
            if typeof(list[i]) <: Expr
                changeRec!(list[i].args)
            elseif typeof(list[i]) <: Symbol
                if haskey(namedTuple, list[i])
                    list[i] = namedTuple[list[i]]
                end
            end
        end
    end # function

    changeRec!(expr.args)
end # function

"""
    evaluate(phenotype::Expr, varsDict::NamedTuple)

Evaluates a phenotype using a named tuple to change the symbols of the variables
into their actual values.

# Arguments
- `phenotype::Expr`: phenotype of a GP tree, obtained by [`getPhenotype`](@ref).
- `varsDict::NamedTuple`: named tuple with the variables used for a GP problem
    and their respective values.

# Returns
The evaluation of the phenotype with the variables' values.

See documentation for [`FitnessFunction`](@ref) for instructions about how to
evaluate GP individuals.
"""
function evaluate(phenotype::Expr, varsDict::NamedTuple)
    newPhenotype = copy(phenotype)
    changeSymbol!(newPhenotype, varsDict)
    return eval(newPhenotype)
end # function

"""
    evaluate(phenotype::Expr, variables::Array{String}, varValues...)

Evaluates a phenotype using a named tuple to change the symbols of the variables
into their actual values.

# Arguments
- `phenotype::Expr`: phenotype of a GP tree, obtained by [`getPhenotype`](@ref).
- `variables::Array{String}`: names of the variables of the problem.
- `varValues`: values of the variables of the problem.

# Returns
The evaluation of the phenotype with the variables' values.

See documentation for [`FitnessFunction`](@ref) for instructions about how to
evaluate GP individuals.
"""
function evaluate(phenotype::Expr, variables::Array{String}, varValues...)
    nVars = length(variables)
    vars = [Meta.parse(x) for x in variables]
    vars = (vars...,)

    varsdict = NamedTuple{vars}(varValues)

    return evaluate(phenotype, varsdict)
end # function

"""
    evaluate(phenotype::Expr)

Evaluates a phenotype of a GP problem that has no variables in it.

# Arguments
- `phenotype::Expr`: phenotype of a GP tree, obtained by [`getPhenotype`](@ref).

# Returns
The evaluation of the phenotype.

See documentation for [`FitnessFunction`](@ref) for instructions about how to
evaluate GP individuals.
"""
function evaluate(phenotype::Expr)
    return eval(phenotype)
end # function

"""
    compareFunctions(genotype::GPGenotype,
                    gpExperimentInfo::GPExperimentInfo, objectives::Array,
                    values...)

This GP-exclusive fitness function compares the objective values given with the
ones obtained by the individual's evaluation using the values for the variables.
The fitness value that this function calculates is the MSE between the results of
the genotype for the given variable values and the objective values.

`Self-provided Arguments` are provided by the library, so only `User Arguments` must be provided.

# Self-provided Arguments
- `genotype::GPGenotype`: genotype of the individual that is going to be evaluated.
- `gpExperimentInfo::GPExperimentInfo`: information about the GP experiment.

# User Arguments
- `objectives::Array`: objective values that are going to be compared.
- `values`: values of the variables. Must be in the same order as the variables
    are in the GP experiment information (see documentation of [`GPExperimentInfo`](@ref)
    and [`getVariables`](@ref)).

# Returns
The MSE between the objective values and the results of evaluating the genotype
with the given variables values.

# Examples
```jdoctests
julia> objectives = [0, 1, 4, 9, 16, 25] #Objective values
6-element Array{Int64,1}:
[...]

julia> xValues = [0, 1, 2, 3, 4 , 5] #Values for variable x
6-element Array{Int64,1}:
[...]

julia> setEvaluator([FitnessFunction(compareFunctions, objectives, xValues, weight=-1)])
Evaluator(...)
```

See also: [`setEvaluator`](@ref), [`FitnessFunction`](@ref)
"""
function compareFunctions(genotype::GPGenotype,
                          gpExperimentInfo::GPExperimentInfo, objectives::Array,
                          values...)

    vars = getVariables(gpExperimentInfo)
    nVars = length(vars)

    vars = [Meta.parse(x) for x in vars]
    vars = (vars...,)
    phenotype = getPhenotype(genotype)

    values = [i for i in values]
    nValues = 1
    y = 0
    y_est = 0
    acc = 0
    failures = 0

    if nVars == 0

        y = objectives[1]

        try
            y_est = evalPhenotype(phenotype)

            if isnan(y_est) || abs(y_est) == Inf
                failures += 1
            else
                acc += (abs(y - y_est))^2
            end
        catch
            failures += 1
        end

    else
        nValues = length(values[1])
        varValues = Array{Number}(undef, nVars)
        for i=1:nValues

            for j=1:nVars
                varValues[j] = values[j][i]
            end
            varsdict = NamedTuple{vars}(varValues)

            y = objectives[i]

            try
                y_est = evaluate(phenotype, varsdict)

                if isnan(y_est) || abs(y_est) == Inf
                    failures += 1
                else
                    acc += (abs(y - y_est))^2
                end
            catch e

                failures += 1
            end
        end
    end

    if failures > (nValues * 0.15)
        fitness = Inf
    else
        mse = acc/(nValues - failures)
        fitness = (acc + failures * mse * 3) / nValues
    end

    return fitness
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia