text:
using Test
using LinearAlgebra
using SparseArrays

using GAFramework
using GAFramework.CoordinateGA
using GAFramework.PermGA
using GAFramework.NetalGA
using GAFramework.MagnaGA

function test1()
    model = FunctionModel(x -> x[1]==0.0 ? 0.0 : x[1] * sin(1/x[1]), [-1.0], [1.0])
    state = GAState(model, ngen=50, npop=300, elite_fraction=0.01,
                    params=Dict(:mutation_rate=>0.1, :print_fitness_iter=>10))
    best = ga!(state)
    x = best.value[1]
    y = best.objvalue
    println("$best $x $y")
    abs(abs(x) - 0.223126) < 0.1 && abs(y - (-0.217219)) < 0.1
end

function test2()
    model = FunctionModel(x -> any(z -> z==0.0, x) ? 0.0 : dot(x, sin.(1 ./x)),
                            -ones(Float64,15), ones(Float64,15))
    # do simulated annealing when mutating
    state = GAState(model, ngen=500, npop=300, elite_fraction=0.01,
        params=Dict(:mutation_rate=>0.1,:sa_rate=>0.1,:sa_k=>1, :sa_lambda=>1/1000,:sa_maxiter=>1000,
            :print_fitness_iter=>50))
    best = ga!(state)
    x = best.value
    y = best.objvalue
    println("$best $x $y")
    all(abs.(abs.(x) .- 0.222549) .< 0.1) && abs(y - (-0.21723*15)) < 0.1
end

function test3()
    model = FunctionModel(x -> 10length(x) + sum(z -> z^2 - 10cos(2pi*z), x),
                            -5.12 .* ones(15), 5.12 .* ones(15))
    state = GAState(model, ngen=500, npop=300, elite_fraction=0.01,
        params=Dict(:mutation_rate=>0.1, :sa_rate=>0.1, :sa_k=>1,
            :sa_lambda=>1/1000, :sa_maxiter=>1000, :print_fitness_iter=>50))
    best = ga!(state)
    x = best.value
    y = best.objvalue
    println("$best $x $y")
    all(abs.(x) .< 0.1) && abs(y) < 0.1
end

function test4()
    rows = [3, 4, 5, 7, 9, 3, 7, 8, 9, 10, 1, 2, 6, 7, 8, 1, 6, 1, 6, 7, 9, 10, 3, 4, 5, 8, 1, 2, 3, 5, 9, 2, 3, 6, 1, 2, 5, 7, 2, 5]
    cols = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10]
    G1 = sparse(rows, cols, ones(length(rows)))
    perm = [3, 2, 4, 8, 6, 1, 5, 10, 9, 7]
    G2 = G1[invperm(perm),invperm(perm)]
    model = NetalignModel(G1,G2)
    st = GAState(model, ngen=100, npop=6_000, elite_fraction=0.1,
        params = Dict(:print_fitness_iter=>20))
    best = ga!(st)
    fitness(best) == 1
end

function test5()
    rows = [3, 4, 5, 7, 9, 3, 7, 8, 9, 10, 1, 2, 6, 7, 8, 1, 6, 1, 6, 7, 9, 10, 3, 4, 5, 8, 1, 2, 3, 5, 9, 2, 3, 6, 1, 2, 5, 7, 2, 5]
    cols = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10]
    G1 = sparse(rows, cols, ones(length(rows)))
    perm = [3, 2, 4, 8, 6, 1, 5, 10, 9, 7]
    G2 = G1[invperm(perm),invperm(perm)]
    model = NetalModel(G1,G2)
    st = GAState(model, ngen=100, npop=6_000, elite_fraction=0.1,
        params = Dict(:print_fitness_iter=>20))
    best = ga!(st)
    fitness(best) == 1
end

function test6()
    rows = [3, 4, 5, 7, 9, 3, 7, 8, 9, 10, 1, 2, 6, 7, 8, 1, 6, 1, 6, 7, 9, 10, 3, 4, 5, 8, 1, 2, 3, 5, 9, 2, 3, 6, 1, 2, 5, 7, 2, 5]
    cols = [1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 10, 10]
    G2 = sparse(rows, cols, ones(length(rows)))
    f = [3, 2, 4, 8, 6, 1, 5, 10, 9, 7]
    G1 = G2[f[1:7],f[1:7]]
    S = rand(7,10)/10
    for i in 1:7 S[i,f[i]] = 1.0 end
    model = MagnaModel(G1,G2,S,0.5)
    st = GAState(model, ngen=100, npop=6_000, elite_fraction=0.1,
        params = Dict(:print_fitness_iter=>20))
    best = ga!(st)
    best.f[1:7] == f[1:7]
end

@test test6()
@test test1()
@test test2()
@test test3()
@test test4()
@test test5()


text:
using Romberg, Test

@testset "Romberg.jl" begin
    # Test interfaces
    x = range(0, stop=π, length=2^8+1)
    y = sin.(x)

    @test romberg(x, y) == romberg(x, y, maxeval=9)

    # Test ispow2(length(x) - 1) == false
    x = range(0, stop=π, length=2^8)
    y = sin.(x)
    @test romberg(x, y)[1] ≈ 2   rtol=1e-12

    x = range(0, stop=π, length=2^8-1)
    y = sin.(x)
    @test romberg(x, y)[1] ≈ 2   rtol=1e-9

    # Test length(x) == 1
    x = 0:0
    y = sin.(x)
    @test romberg(x, y)[1] == 0

    # Test length(x) == 2
    x = range(0, stop=1, length=2)
    y = x.^2
    @test romberg(x,y)[1] == 0.5   # integration is inaccurate

    # Test length(x) == 3
    x = range(0, stop=1, length=3)
    y = x.^2
    @test romberg(x, y)[1] ≈ 1/3   rtol=3e-16  # should be exact up to roundoff error

    # Integrate different functions
    x = range(0, stop=π, length=2^8+1)
    y = sin.(x)
    @test romberg(x, y)[1] ≈ 2    rtol=1e-15

    x = range(0, stop=1, length=2^8+1)
    y = x.^3
    @test romberg(x, y)[1] ≈ 1/4   rtol=3e-16  # should be exact up to roundoff error

    x = range(0, stop=π/2, length=2^5+1)
    y = sin.(x).^2
    @test romberg(x, y)[1] ≈ π/4    rtol=1e-15

    m = 3
    n = 4
    x = range(0, stop=π, length=2^8+1)
    y = sin.(m*x).*cos.(n*x)
    v = romberg(x, y)
    @test v[1] ≈ 2*m/(m^2 - n^2)  rtol=1e-13
    @test v[2] < 1e-10

    # the following are singular integrands where the underlying
    # theory of Romberg integration starts to break down:

    # this one requires lots of samples...
    a = 15
    x = range(0, stop=a, length=2^16+1)
    y = sqrt.(a^2 .- x.^2)
    @test romberg(x, y)[1] ≈ π*a^2/4   rtol=1e-8
    # we can do much better if we put in the correct convergence
    # rate for this singular integrand, but this requires some
    # understanding of the theory that most people won't have…
    x = range(0, stop=a, length=2^7+1)
    y = sqrt.(a^2 .- x.^2)
    @test romberg(x, y, power=0.5)[1] ≈ π*a^2/4   rtol=1e-8

    # tricky to integrate
    x = range(1e-15, stop=1, length=2^16+1)
    y = log.(x)./(1 .+ x)
    @test romberg(x, y)[1] ≈ -π^2/12  atol=1e-4

    # make sure it works for abstractly-typed y and integer Δx
    @test romberg(1, Any[3,3,3,3,3,3,3])[1] == 6*3

    @test @inferred(romberg(1, [1//2, 1//4])) === (0.375, 0.375)
    @test @inferred(romberg(1, [[1//2,1//1], [1//4,1//1]])) == ([0.375,1.0], hypot(0.375,1.0))
end

ext: jl


text:
using Random
SansSerifFonts = ["Trebuchet MS", "Heiti TC", "微軟正黑體", "Arial Unicode MS", "Droid Fallback Sans", "sans-serif", "Helvetica", "Verdana", "Hei",
"Arial", "Tahoma", "Trebuchet MS", "Microsoft Yahei", "Comic Sans MS", "Impact", "Segoe Script", "STHeiti", "Apple LiGothic", "MingLiU", "Ubuntu",
"Segoe UI", "DejaVu Sans", "DejaVu Sans Mono", "Noto Sans CJK"]
SerifFonts = ["Baskerville", "Times New Roman", "華康儷金黑 Std", "華康儷宋 Std",  "DFLiKingHeiStd-W8", "DFLiSongStd-W5", "DejaVu Serif", "SimSun",
    "Hiragino Mincho Pro", "LiSong Pro", "新細明體", "serif", "Georgia", "STSong", "FangSong", "KaiTi", "STKaiti", "Courier New"]
CandiFonts = union(SansSerifFonts, SerifFonts)
CandiWeights = ["", " Regular", " Normal", " Medium", " Bold", " Light"]
function checkfonts(fonts::AbstractVector)
    fname = tempname()
    r = Bool[]
    open(fname, "w") do f
        redirect_stderr(f) do
            p = position(f)
            for font in fonts
                rendertext("a", 1 + rand(), font=font) # 相同字体相同字号仅warning一次，故首次执行最准
                # flush(f) #https://en.cppreference.com/w/cpp/io/c/fseek The standard C++ file streams guarantee both flushing and unshifting
                seekend(f)
                p2 = position(f)
                push!(r, p2 == p)
                p = p2
            end
        end
    end
    return r
end
checkfonts(f) = checkfonts([f]) |> only
function filterfonts(;fonts=CandiFonts, weights=CandiWeights)
    candi = ["$f$w" for w in weights, f in fonts] |> vec
    candi[checkfonts(candi)]
end
AvailableFonts = filterfonts()
push!(AvailableFonts, "")

Schemes_colorbrewer = filter(s -> occursin("colorbrewer", colorschemes[s].category), collect(keys(colorschemes)))
Schemes_colorbrewer =  filter(s -> (occursin("Accent", String(s))
        || occursin("Dark", String(s))
        || occursin("Paired", String(s))
        || occursin("Pastel", String(s))
        || occursin("Set", String(s))
        || occursin("Spectral", String(s))
        || occursin("BuPu", String(s))
        || occursin("PuOr", String(s))
        || occursin("RdPu", String(s))
        ), Schemes_colorbrewer)
Schemes_seaborn = filter(s -> occursin("seaborn", colorschemes[s].category), collect(keys(colorschemes)))
Schemes_tableau = filter(s -> occursin("tableau", colorschemes[s].category), collect(keys(colorschemes)))
Schemes_cvd = filter(s -> occursin("cvd", colorschemes[s].category), collect(keys(colorschemes)))
Schemes = [Schemes_colorbrewer..., Schemes_seaborn..., Schemes_tableau..., Schemes_cvd...]

function displayschemes()
    for scheme in Schemes
        display(scheme)
        colors = Render.colorschemes[scheme].colors
        display(colors)
    end
end
function randomscheme()
    if rand() < 0.95
        scheme = rand(Schemes)
        c = Render.colorschemes[scheme].colors
        colors = randsubseq(c, rand())
        colors = isempty(colors) ? c : colors
        @show (scheme, length(colors))
    else
        colors = rand((0, 1, 0, 1, 0, 1, (0, 1), rand(), (rand(), rand())))
        @show colors
    end
    (colors...,)
end
function randommask(sz::Number=800; kargs...)
    s = sz * sz * (0.5 + rand() / 2)
    ratio = (0.5 + rand() / 2)
ratio = ratio > 0.9 ? 1.0 : ratio
    h = round(Int, sqrt(s * ratio))
    w = round(Int, h / ratio)
    randommask(w, h; kargs...)
end
function randommask(sz; kargs...)
    randommask(sz...; kargs...)
end
function randommask(w, h, args...; maskshape=:rand, kargs...)
    ran = Dict(box => 0.2, squircle => 0.7, ellipse => 1, :rand => rand())[maskshape]
    if ran <= 0.2
        return randombox(w, h, args...; kargs...)
    elseif ran <= 0.7
        return randomsquircle(w, h, args...; kargs...)
    else
        return randomellipse(w, h, args...; kargs...)
    end
end
function randombox(w, h, r=:rand; kargs...)
    if r == :rand
        r = rand() * 0.5 - 0.05 # up to 0.45
        r = r < 0. ? 0. : r # 10% for 0.
        r = round(Int, h * r)
    end
    println("shape(box, $w, $h, $r", join([", $k=$(repr(v))" for (k, v) in kargs]), ")")
    return shape(box, w, h, r; kargs...)
end
function randomsquircle(w, h; rt=:rand, kargs...)
    if rt == :rand
        if rand() < 0.8
            rt = rand()
        else
            ran = rand()
            if ran < 0.5
                rt = 2
            else
    rt = 1 + 1.5rand()
            end
        end
    end
    println("shape(squircle, $w, $h, rt=$rt", join([", $k=$(repr(v))" for (k, v) in kargs]), ")")
    return shape(squircle, w, h, rt=rt; kargs...)
end
function randomellipse(w, h; kargs...)
    println("shape(ellipse, $w, $h", join([", $k=$(repr(v))" for (k, v) in kargs]), ")")
    return shape(ellipse, w, h; kargs...)
end
function randomangles()
    a = rand((-1, 1)) .* rand((0, (0, 90), (0, 90, 45), (0, 90, 45, -45), (0, 45, -45), (45, -45), -90:90, 0:90))
    println("angles = ", a)
    a
end
function randommaskcolor(colors)
    colors = parsecolor.(unique(colors))
    try
        g = Gray.(colors) |> sort
        m = g[1]
        M = g[end]
        if length(g) > 1
            d = diff(g)
            I = maximum(d)
            i = findlast(isequal(I), d)
        else
            I = 0
            i = -1
        end
        # @show I, m, M
        if I > 3(1 - M) && I > 3m
            middle = (g[i] + g[i + 1]) / 2
            th1 = clamp(max(g[i] + 0.15, middle - rand(0:0.001:0.05)), 0, middle)
            th2 = clamp(min(g[i + 1] - 0.15, middle + rand(0:0.001:0.05)), middle, 1)
            default = middle
        elseif sum(g) / length(g) < 0.7 && (m + M) / 2 < 0.7 && !(m > 2(1 - M))# 明亮
            th1 = clamp(max(M + 0.15, rand(0.85:0.001:1.0)), 0, 1)
            th2 = clamp(th1 + 0.1, 0, 1)
            default = 1.0
        else    # 黑暗
            th2 = clamp(min(m - 0.15, rand(0.0:0.001:0.3)), 0, 1) # 对深色不敏感，+0.15
            th1 = clamp(th2 - 0.15, 0, 1)
            default = 0.0
        end
        maskcolor = rand((default, (rand(th1:0.001:th2), rand(th1:0.001:th2), rand(th1:0.001:th2))))
        # @show maskcolor
        return maskcolor
    catch e
        @show e
    @show "colors sum failed", colors
        return "white"
    end
end
function randomlinecolor(colors)
    if rand() < 0.8
        linecolor = rand((colors[1], colors[1], rand(colors)))
    else
        linecolor = (rand(), rand(), rand(), min(1., 0.5 + rand() / 2))
    end
    linecolor
end
randomoutline() = rand((0, 0, 0, rand((1, 2, 3, 4, 5))))
function randomfont()
    font = rand(AvailableFonts)
    @show font
    font
end

max_stars_repo_name: guo-yong-zhi/WordCloud

max_issues_repo_name: guo-yong-zhi/WordCloud

max_forks_repo_name: guo-yong-zhi/WordCloud


text:
module LabeledArrays

### A beneficial type piracy, todo: pull-request in SparseArrays
using SparseArrays
import SparseArrays: nnz
using LinearAlgebra
import LinearAlgebra: Adjoint
"""
    nnz(x::LinearAlgebra.Adjoint)

Call `SparseArrays.nnz` of transpose for performance.
"""
nnz(x::LinearAlgebra.Adjoint) =
    nnz(x.parent)

############################################################
export NBijection, EnumerationDict, VectorDict
"""
    NBijection(x::Vector{T})

KEY <-> Int, bijection map. Synonym to possibly better name `EnumerationDict`, or `VectorDict`.

Store inverse indices in a `Dict{T,I}` as well as original `Vector{T}`.
- `copy`: set to false to store `x` without copying.
  (make sure that Vector will be mutated only through `NBijection` API.
"""
struct NBijection{KEY, I<:Integer}
    enumeration::Vector{KEY}
    index::Dict{KEY,I}
end
EnumerationDict{KEY, I<:Integer} = NBijection{KEY, I}
VectorDict{KEY, I<:Integer} = NBijection{KEY, I}
NBijection{KEY, I}() where {I,KEY} =
    NBijection{KEY,I}(KEY[], Dict{KEY,I}())
NBijection(x) = NBijection(collect(x), copy=false)
NBijection(x::Vector{KEY}; copy=true) where KEY =
    NBijection{KEY,Int}(
        copy ? Base.copy(x) : x,
        Dict((t => i for (t,i) = zip(x, 1:length(x)))...))

@deprecate terms(x::NBijection) keys(x)
"""
    keys(b::NBijection)

terms are keys in sequence that share position index with some other
data structure, like a count matrix in Count.
"""
Base.keys(b::NBijection) = b.enumeration
Base.values(b::NBijection) = 1:length(b)
Base.keytype(x::NBijection{KEY,I}) where {KEY,I} = KEY
Base.valtype(x::NBijection{KEY,I}) where {KEY,I} = I
Base.eltype(x::NBijection{KEY,I}) where {KEY,I} = KEY
Base.show(io::IO,x::NBijection) =
    print(io,"NBijection{",valtype(x),"}: ",x.index)
Base.length(b::NBijection) = length(b.enumeration)
Base.iterate(b::NBijection,a...) =
    iterate(b.enumeration,a...)
"""
    Base.getindex(b::NBijection,i::Integer)

`b.enumeration[i]`.
"""
@inline Base.@propagate_inbounds Base.getindex(b::NBijection,i::Integer) =
    b.enumeration[i]

import Base: match
Base.match(x::NBijection, y) = false

import Base: indexin
Base.indexin(i,b::NBijection{KEY,I}) where {KEY,I} =
    Union{I,Nothing}[ get(b.index,e,nothing) for e in i]

import Base: get!, push!
"""
    Base.get!(f::Function,b::NBijection, key)

Calls `f(length(b),key)` hook function if `key !in b`, push `key` to `b.enumeration` and store `key=>index`.
Similar to `get!(f,b::Dict, key)`.

TODO: remove f(n,key)
```jldoctest
julia> x = NBijection(['a','b','c'])
julia> get!(x,'b')
2

julia> get!((i,k)->println("neu"), x, 'd')
neu
4

julia> get!((i,k)->println("neu"), x, 'd')
4
```
"""
Base.get!(f::Function,b::NBijection, key) =
    get!(b.index, key) do
        push!(b.enumeration, key)
        f(length(b),key)
        length(b)
    end

"""
    Base.get!(b::NBijection, key)

When `key !in b.index`, push `nothing` to `b.enumeration` and store `key=>index`.
```jldoctest
julia> x = NBijection(['a','b','c'])
julia> get!(x,'b')
2

julia> get!(x,'d')
4

julia> get!(x,'d')
4

```
"""
Base.get!(b::NBijection, key) =
    get!((_,_)->nothing,b, key)
@deprecate get_index!(b::NBijection, key) get!(b,key)

"""
    Base.push!(b::NBijection, key)

Push `key` to `b.enumeration` and store `key=>index`.
Return last index.

@asserts that key is not in index, fails otherwise.
Consider [`get!`](@ref).
"""
function Base.push!(b::NBijection, key)
    @assert !haskey(b.index,key)
    push!(b.enumeration,key)
    b.index[key] = length(b.enumeration)
end

"""
    sparse_vector(x::NBijection,sel)

Create a sparse vector representation of `sel`ected values in `x::`[`NBijection`](@ref).
"""
function SparseArrays.sparsevec(x::NBijection,sel)
    idxs = filter(
        x->x!==nothing,
        indexin(sel,x))
    sparsevec(
        idxs,
        ConstArray(1,length(idxs)),
        length(x))
end

export ConstArray
struct ConstArray{T,N,l} <: AbstractArray{T,N}
    value::T
    ConstArray(x::T, l::Integer...) where T =
        new{T,length(l),l}(x)
end
Base.size(::ConstArray{T,N,l}) where {T,N,l}	= l
Base.getindex(A::ConstArray, i::Int...) = A.value
Base.setindex!(A::ConstArray, v, i::Int) = error("create a normal array and convert")
Base.sum(A::ConstArray) = A.value * (*(size(A)...))

export LabeledMatrix
"""
A matrix, `AbstractArray{T,2}`, with a layer labeled rows and cols.
"""
struct LabeledMatrix{T,A <: AbstractArray{T,2},R,C} <: AbstractArray{T,2}
    row::R
    col::C
    values::A
    function LabeledMatrix(row,col,values::A) where {A<:AbstractArray}
        new{eltype(values),A,typeof(row),typeof(col)}(row,col,values)
    end
end
export collapse_observations
function collapse_observations(f::Function,x,T=Any)
   nb = NBijection(T[])
    cols = Int32[  get!(nb,f(e))
                   for e in values(x)
                   ]
    LabeledMatrix(x,nb,
                 sparse(1:length(cols),cols,
                        ConstArray(1,length(cols)),
                        length(cols),length(nb)))
end
export sparsevec
SparseArrays.sparsevec(x::LabeledMatrix,y) = sparsevec(x.col,y)
@deprecate row(x,y) sparsevec(x,y)
# ? col, row ?
# SparseArrays.sparsecol(x::LabeledMatrix,y) = sparsevec(x.row,y)

Base.getindex(x::LabeledMatrix, a...) = x.values[a...]
Base.getindex(x::LabeledMatrix, row, a...) =
    if match(x.row,row)
        x.values[get(row),a...]
    else
        error("invalid index: $row for $(x.row)")
    end

import Base: +, *, adjoint, size
Base.size(x::LabeledMatrix) =
    Base.size(x.values)
adjoint(x::LabeledMatrix) =
    LabeledMatrix(x.col, x.row, x.values')
@deprecate reverse(x::LabeledMatrix) adjoint(x)

function *(x::LabeledMatrix, y::LabeledMatrix)
    @assert x.col == y.row
    LabeledMatrix(
        x.row, y.col,
        x.values * y.values)
end
*(x::LabeledMatrix, y::AbstractArray) = x.values*y
*(x::AbstractArray, y::LabeledMatrix) = x*y.values
*(x::LabeledMatrix, y::Number) =
    LabeledMatrix(
        x.row, x.col,
        x.values*y)
*(x::Number, y::LabeledMatrix) =
    LabeledMatrix(
        y.row, y.col,
        x*y.values)

function +(x::LabeledMatrix, y::LabeledMatrix)
    @assert x.row == y.row
    @assert x.col == y.col
    LabeledMatrix(
        x.row, y.col,
        x.values + y.values)
end
+(x::LabeledMatrix, y) = x.values+y
+(x, y::LabeledMatrix) = x+y.values
+(x::LabeledMatrix, y::Number) =
    LabeledMatrix(
        x.row, x.col,
        x.values+y)
+(x::Number, y::LabeledMatrix) =
    LabeledMatrix(
        y.row, y.col,
        x+y.values)

import Base: isless
Base.isless(than::Int, A::LabeledMatrix) =
    LabeledMatrix(A.row, A.col, A.values.>than)

Base.sum(x::LabeledMatrix; kw...) =
    sum(x.values; kw...)
using LinearAlgebra
import LinearAlgebra: Adjoint
Base.sum(x::Adjoint; dims=:) =
    if dims isa Colon
        sum(x.parent,dims=:)'
    elseif dims==1
        sum(x.parent,dims=2)'
    else
        sum(x.parent,dims=1)'
    end

import Printf: @printf
function Base.show(io::IO, x::LabeledMatrix)
    Base.show(io, MIME"text/plain"(), x)
end
compactstring(io::IO, x) =
    printstyled(IOContext(io, :compact => true),x,color=:red)

function show_stats(x::LabeledMatrix)
    nzs = x.values.>0
    S=sum(x.values,dims=2)
    (nnz=nnz(x.values)
     , relnz = (nnz(x.values)/*(size(x.values)...))
     , mean_row_sum = sum(S[S.>0])/length(S) #?
     , unique = sum(sum(nzs,dims=2))/size(nzs)[2] #?
     , zeros_1 = sum(iszero.(sum(nzs,dims=1)))
     , zeros_2 = sum(iszero.(sum(nzs,dims=2)))
     )
end

function Base.show(io::IO, ::MIME"text/plain", x::LabeledMatrix)
    rel="root-word"
    s = show_stats(x)
    print(io, size(x.values)[1]," * ", size(x.values)[2])
    printstyled(io," ", x.row, color=:red)
    printstyled(io," with $rel ", color=:yellow)
    printstyled(io," ", x.col, color=:red)
    @printf(io, "\n%d non-zero (~%2.4f%% sparse). ",
            s.nnz, (1.0-s.relnz)*100.0)
    @printf(io, "\nAverage %s row-sum %.2f (unique %.2f) ",
            rel,
            s.mean_row_sum,
            s.unique)
    compactstring(io,valtype(x.col))
    print(io,"s / ")
    compactstring(io,valtype(x.row))
    zeros=s.zeros_2
    print(io,". ")
    if zeros > 0
        print(io,"\n$zeros ")
        compactstring(io,valtype(x.row))
        print(io, "s have all-zero ")
        compactstring(io,valtype(x.col))
        print(io,". ")
    end
    zeros=s.zeros_1
    if zeros > 0
        print(io,"\n$zeros ")
        compactstring(io,valtype(x.col))
        print(io, "s occur in no ")
        compactstring(io,valtype(x.row))
        print(io,". ")
    end
end

end # module


text:
@info "Starting precompile exec."

#import Pluto

#session = Pluto.ServerSession()

nb_path = joinpath(@__DIR__(), "notebook.jl")
include(nb_path)

#@info "Loading sample notebook"
#nb = Pluto.load_notebook(Pluto.tamepath(nb_path));
#session.notebooks[nb.notebook_id] = nb;

#@info "Running sample notebook"
#Pluto.update_save_run!(session, nb, nb.cells; run_async=false, prerender_text=true)

# nice! we ran the notebook, so we already precompiled a lot

@info "Precompile exec finished."

max_stars_repo_name: greimel/pluto-binder

max_issues_repo_name: greimel/pluto-binder

max_forks_repo_name: greimel/pluto-binder


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

@show gethostname()

rep   = parse(Int,ARGS[1])
ip    = parse(Int,ARGS[2])
ialpha = parse(Int,ARGS[3])
dir = ARGS[4]

alpha = 0.1:0.1:1
pArr = [50, 100, 150]
n = 300

bnm = Binomial(1, 0.05)

p = pArr[ip]
Random.seed!(54298)

# generate X
Ωx = Matrix{Float64}(I, p, p)
for i = 1:(p-1)
  for j = (i+1):p
    Ωx[i, j] = rand(bnm) * 0.8
    Ωx[j, i] = Ωx[i, j]
  end
end
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Ωx = D * Ωx * D

level = maximum( diag(Ωx) ) * sqrt(2. * log(p) / n)

bc = Bernoulli()
Ωy = copy(Ωx)
rd = sample(1:p-1, 8; replace=false)
sort!(rd)
Ωy[rd[1], rd[5]] += (2 * rand(bc) - 1) * level
Ωy[rd[2], rd[6]] += (2 * rand(bc) - 1) * level
Ωy[rd[3], rd[7]] += (2 * rand(bc) - 1) * level
Ωy[rd[4], rd[8]] += (2 * rand(bc) - 1) * level

Ωy[rd[5], rd[1]] = Ωy[rd[1], rd[5]]
Ωy[rd[6], rd[2]] = Ωy[rd[2], rd[6]]
Ωy[rd[7], rd[3]] = Ωy[rd[3], rd[7]]
Ωy[rd[8], rd[4]] = Ωy[rd[4], rd[8]]

de = max( abs(minimum(eigvals(Ωx))), abs(minimum(eigvals(Ωy))) ) + 0.05
Ωx = (Ωx + de*I)/(1+de)
Ωy = (Ωy + de*I)/(1+de)

Ωy = (1. - alpha[ialpha]) * Ωx + alpha[ialpha] * Ωy
Δ = Ωx - Ωy

Σx = inv(Symmetric(Ωx))
Σy = inv(Symmetric(Ωy))
dist_X = MvNormal(convert(Matrix, Σx))
dist_Y = MvNormal(convert(Matrix, Σy))

# generate data
Random.seed!(12346 + rep*100 + ialpha)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

Sx = Symmetric( X'X / n )
Sy = Symmetric( Y'Y / n )

S = BitArray( abs(Δ[r, c]) > 1e-4 for r = 1:p, c = 1:p )

eS_init = Array{BitArray}(undef, div((p + 1)*p, 2))
for j=1:div((p + 1)*p, 2)
    eS_init[j] = S
end

@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init)

@save "$(dir)/res_$(rep).jld" boot_res


text:
function read_request(channel::gRPCChannel, controller::gRPCController, services)
    @debug("before connection")
    connection = channel.session
    @debug("before evt")
    evt = Session.take_evt!(connection)

    @debug("before ifisa")
    if isa(evt, EvtGoaway)
        @debug("received EvtGoaway")
        close(channel)
        return (nothing, nothing, nothing)
    end

    @debug("before notifisa")
    if !isa(evt, EvtRecvHeaders)
        @debug("unexpected event while reading request (closing channel)", evt)
        close(channel)
        return (nothing, nothing, nothing)
    end

    channel.stream_id = evt.stream_identifier

    headers = evt.headers
    method = headers[":method"]
    path = headers[":path"]
    pathcomps = split(path, "/"; keep=false)
    if length(pathcomps) != 2
        @debug("unexpected path while reading request (closing channel)", path)
        close(channel)
        return (nothing, nothing, nothing)
    end

    servicename, methodname = split(path, "/"; keep=false)

    if evt.is_end_stream
        data = UInt8[]
    else
        data_evt = Session.take_evt!(connection)
        data = data_evt.data
    end

    @debug("received request", method, path, stream_id=channel.stream_id, servicename, methodname, nbytes=length(data))

    service = services[servicename]
    method = find_method(service, methodname)
    request_type = get_request_type(service, method)
    request = request_type()
    from_delimited_message_bytes(data, request)

    service, method, request
end

function write_response(channel::gRPCChannel, controller::gRPCController, response)
    @debug("write_response")
    sending_headers = [(":status", "200")]
    data_buff = to_delimited_message_bytes(response)

    Session.put_act!(channel.session, Session.ActSendHeaders(channel.stream_id, sending_headers, false))
    Session.put_act!(channel.session, Session.ActSendData(channel.stream_id, data_buff, true))
    nothing
end

function call_method(channel::gRPCChannel, service::ServiceDescriptor, method::MethodDescriptor, controller::gRPCController, request)
    @debug("call_method")
    write_request(channel, controller, service, method, request)
    response_type = get_response_type(method)
    response = response_type()
    read_response(channel, controller, response)
end

# gRPC server implementation
#
mutable struct gRPCServer
    sock::TCPServer
    services::Dict{String, ProtoService}
    run::Bool

    gRPCServer(services::Tuple{ProtoService}, ip::IPv4, port::Integer) = gRPCServer(services, listen(ip, port))
    gRPCServer(services::Tuple{ProtoService}, port::Integer) = gRPCServer(services, listen(port))
    function gRPCServer(services::Tuple{ProtoService}, sock::TCPServer)
        svcdict = Dict{String,ProtoService}()
        for svc in services
            svcdict[svc.desc.name] = svc
        end
        new(sock, svcdict, true)
    end
end

# TODO: close all channels, wait for/interrupt processors
close(srvr::gRPCServer) = close(srvr.sock)

function process(controller::gRPCController, srvr::gRPCServer, channel::gRPCChannel)
    @info("start processing channel")
    try
        while(!channel.session.closed)
            @info("start processing channel 0")
            service, method, request = read_request(channel, controller, srvr.services)
            @info("start processing channel 1")
            (service === nothing) && continue
            @info("start processing channel 2")

            response = call_method(service, method, controller, request)
            @info("start processing channel 3")
            #@debug("response from method", response)
            write_response(channel, controller, response)
            @info("start processing channel 4")
        end
    catch ex
        @warn("channel stopped with exception", ex)
    end
    # TODO: close channel if not closed, remove reference from srvr
    @info("stopped processing channel")
end

function run(srvr::gRPCServer)
    controller = gRPCController()
    try
        while(srvr.run)
            buffer = accept(srvr.sock)
            connection = Session.new_connection(buffer; isclient=false)

            # TODO: keep channel reference in server
            channel = gRPCChannel(connection)
            @async process(controller, srvr, channel)
        end
    catch ex
        @warn("server stopped with exception", ex)
    end
    @info("stopped server")
end


text:
function avg(x)
    sum(x)/length(x)
end


text:
# This file was generated, do not modify it. # hide
mode.(ŷ)


text:
# Copyright (c) 2013: Cbc.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

const MOI = MathOptInterface

MOI.Utilities.@product_of_sets(
    _LPProductOfSets,
    MOI.EqualTo{T},
    MOI.LessThan{T},
    MOI.GreaterThan{T},
    MOI.Interval{T},
)

MOI.Utilities.@struct_of_constraints_by_set_types(
    _CbcConstraints,
    Union{MOI.EqualTo{T},MOI.LessThan{T},MOI.GreaterThan{T},MOI.Interval{T}},
    MOI.SOS1{T},
    MOI.SOS2{T},
)

const OptimizerCache = MOI.Utilities.GenericModel{
    Float64,
    MOI.Utilities.ObjectiveContainer{Float64},
    MOI.Utilities.VariablesContainer{Float64},
    _CbcConstraints{Float64}{
        MOI.Utilities.MatrixOfConstraints{
            Float64,
            MOI.Utilities.MutableSparseMatrixCSC{
                Float64,
                Cint,
                MOI.Utilities.ZeroBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{Float64},
            _LPProductOfSets{Float64},
        },
        MOI.Utilities.VectorOfConstraints{
            MOI.VectorOfVariables,
            MOI.SOS1{Float64},
        },
        MOI.Utilities.VectorOfConstraints{
            MOI.VectorOfVariables,
            MOI.SOS2{Float64},
        },
    },
}

"""
    Optimizer()

Create a new Cbc Optimizer.
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Ptr{Cvoid}
    silent::Bool
    params::Dict{String,String}
    variable_start::Dict{MOI.VariableIndex,Float64}
    objective_constant::Float64
    solve_time::Float64
    termination_status::Cint
    has_solution::Bool
    variable_primal::Union{Nothing,Vector{Float64}}
    constraint_primal::Union{Nothing,Vector{Float64}}

    function Optimizer()
        model = new(
            Cbc_newModel(),
            false,
            Dict{String,String}(),
            Dict{MOI.VariableIndex,Float64}(),
            0.0,
            0.0,
            Cint(-1),
            false,
            nothing,
            nothing,
        )
        finalizer(Cbc_deleteModel, model)
        return model
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, model::Optimizer) = model
Base.unsafe_convert(::Type{Ptr{Cvoid}}, model::Optimizer) = model.inner

function MOI.default_cache(::Optimizer, ::Type{Float64})
    return MOI.Utilities.UniversalFallback(OptimizerCache())
end

function MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute)
    # TODO(odow): There is no programatical way throught the C API to check if a
    # parameter name (or value) is valid. Fix this upstream.
    return true
end

function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return MOI.set(model, param, string(value))
end

function MOI.set(
    model::Optimizer,
    param::MOI.RawOptimizerAttribute,
    value::String,
)
    if !MOI.supports(model, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    model.params[param.name] = value
    if param.name == "threads" && Sys.iswindows()
        @warn(
            "Ignoring threads parameter due to known bugs in CBC. Read " *
            "https://github.com/jump-dev/Cbc.jl/issues/186 for more details.",
        )
        return
    end
    if !(model.silent && param.name == "logLevel")
        Cbc_setParameter(model, param.name, value)
    end
    return
end

function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    # TODO: This gives a poor error message if the name of the parameter is
    # invalid.
    return model.params[param.name]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    if value == model.silent
        return
    end
    log_level = value ? "0" : get(model.params, "logLevel", "1")
    Cbc_setParameter(model, "logLevel", log_level)
    model.silent = value
    return
end
MOI.get(model::Optimizer, ::MOI.Silent) = model.silent

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    if value === nothing
        delete!(model.params, "seconds")
        Cbc_setParameter(model, "seconds", "??")
    else
        MOI.set(model, MOI.RawOptimizerAttribute("seconds"), value)
    end
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    value = get(model.params, "seconds", nothing)
    return value === nothing ? value : parse(Float64, value)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "COIN Branch-and-Cut (Cbc)"

MOI.get(::Optimizer, ::MOI.SolverVersion) = unsafe_string(Cbc_getVersion())

function MOI.empty!(model::Optimizer)
    Cbc_deleteModel(model)
    model.inner = Cbc_newModel()
    model.objective_constant = 0.0
    model.termination_status = Cint(-1)
    model.solve_time = 0.0
    for (name, value) in model.params
        MOI.set(model, MOI.RawOptimizerAttribute(name), value)
    end
    if model.silent
        Cbc_setParameter(model, "logLevel", "0")
    end
    empty!(model.variable_start)
    model.has_solution = false
    model.variable_primal = nothing
    model.constraint_primal = nothing
    return
end

function MOI.is_empty(model::Optimizer)
    return Cbc_getNumCols(model) == 0 && Cbc_getNumRows(model) == 0
end

function _index_map(
    src::OptimizerCache,
    ::MOI.IndexMap,
    ci::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction},
)
    return MOI.Utilities.rows(src.constraints.moi_equalto, ci)
end

function _index_map(
    ::OptimizerCache,
    index_map::MOI.IndexMap,
    ci::MOI.ConstraintIndex{<:MOI.VariableIndex},
)
    return index_map[MOI.VariableIndex(ci.value)].value
end

function _index_map(
    ::OptimizerCache,
    ::MOI.IndexMap,
    ci::MOI.ConstraintIndex{<:MOI.VectorOfVariables},
)
    return ci.value
end

function _index_map(
    src::OptimizerCache,
    index_map::MOI.IndexMap,
    ::Type{F},
    ::Type{S},
) where {F,S}
    inner = index_map.con_map[F, S]
    for ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        inner[ci] = MOI.ConstraintIndex{F,S}(_index_map(src, index_map, ci))
    end
    return
end

"""
    _index_map(src::OptimizerCache)

Create an `IndexMap` mapping the variables and constraints in `OptimizerCache`
to their corresponding 1-based columns and rows.
"""
function _index_map(src::OptimizerCache)
    index_map = MOI.IndexMap()
    for (i, x) in enumerate(MOI.get(src, MOI.ListOfVariableIndices()))
        index_map[x] = MOI.VariableIndex(i)
    end
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        _index_map(src, index_map, F, S)
    end
    return index_map
end

function _constraint_matrix(constraints, n)
    @assert n == constraints.coefficients.n
    return (
        m = constraints.coefficients.m,
        n = constraints.coefficients.n,
        colptr = constraints.coefficients.colptr,
        rowval = constraints.coefficients.rowval,
        nzval = constraints.coefficients.nzval,
        lower = constraints.constants.lower,
        upper = constraints.constants.upper,
    )
end

function _constraint_matrix(::Nothing, n)
    return (
        m = Cint(0),
        n = Cint(n),
        colptr = fill(Cint(0), n + 1),
        rowval = Cint[],
        nzval = Float64[],
        lower = Float64[],
        upper = Float64[],
    )
end

function MOI.copy_to(dest::Optimizer, src::OptimizerCache)
    matrix = _constraint_matrix(
        src.constraints.moi_equalto,
        MOI.get(src, MOI.NumberOfVariables()),
    )
    c = zeros(matrix.n)
    obj =
        MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    for term in obj.terms
        c[term.variable.value] += term.coefficient
    end
    dest.objective_constant = obj.constant
    zeroone_attr = MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}()
    binaries = Cint[Cint(ci.value - 1) for ci in MOI.get(src, zeroone_attr)]
    variable_lower = copy(src.variables.lower)
    variable_upper = copy(src.variables.upper)
    for b in binaries
        variable_lower[b+1] = max(variable_lower[b+1], 0.0)
        variable_upper[b+1] = min(variable_upper[b+1], 1.0)
    end
    Cbc_loadProblem(
        dest,
        matrix.n,
        matrix.m,
        matrix.colptr,
        matrix.rowval,
        matrix.nzval,
        variable_lower,
        variable_upper,
        c,
        matrix.lower,
        matrix.upper,
    )
    sense = MOI.get(src, MOI.ObjectiveSense())
    if sense == MOI.MIN_SENSE
        Cbc_setObjSense(dest, 1)
    elseif sense == MOI.MAX_SENSE
        Cbc_setObjSense(dest, -1)
    else
        @assert sense == MOI.FEASIBILITY_SENSE
        Cbc_setObjSense(dest, 0)
    end
    Cbc_setInteger.(dest, binaries)
    attr = MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}()
    for ci in MOI.get(src, attr)
        Cbc_setInteger(dest, Cint(ci.value - 1))
    end
    any_sos = false
    for (S, type) in ((MOI.SOS1{Float64}, 1), (MOI.SOS2{Float64}, 2))
        starts, indices, weights = Cint[], Cint[], Float64[]
        attr = MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S}()
        for ci in MOI.get(src, attr)
            any_sos = true
            push!(starts, Cint(length(weights)))
            f = MOI.get(src, MOI.ConstraintFunction(), ci)
            for x in f.variables
                push!(indices, Cint(x.value - 1))
            end
            s = MOI.get(src, MOI.ConstraintSet(), ci)
            append!(weights, s.weights)
        end
        N = Cint(length(starts))
        if N > 0
            push!(starts, length(weights))
            Cbc_addSOS(dest, N, starts, indices, weights, Cint(type))
        end
    end
    if any_sos && Cbc_getNumIntegers(dest) == 0
        @warn(
            "There are known correctness issues using Cbc with SOS " *
            "constraints and no binary variables.",
        )
    end
    return _index_map(src)
end

function MOI.copy_to(
    dest::Optimizer,
    src::MOI.Utilities.UniversalFallback{OptimizerCache},
)
    attr = MOI.VariablePrimalStart()
    MOI.Utilities.throw_unsupported(
        src;
        excluded_attributes = Any[MOI.VariablePrimalStart()],
    )
    index_map = MOI.copy_to(dest, src.model)
    if attr in MOI.get(src, MOI.ListOfVariableAttributesSet())
        for (x_src, x_dest) in index_map.var_map
            value = MOI.get(src, attr, x_src)
            MOI.set(dest, attr, x_dest, value)
        end
    end
    return index_map
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
    cache = MOI.default_cache(dest, Float64)
    src_cache = MOI.copy_to(cache, src)
    cache_dest = MOI.copy_to(dest, cache)
    index_map = MOI.IndexMap()
    for (src_x, cache_x) in src_cache.var_map
        index_map[src_x] = cache_dest[cache_x]
    end
    for (src_ci, cache_ci) in src_cache.con_map
        index_map[src_ci] = cache_dest[cache_ci]
    end
    return index_map
end

###
### supports and supports_constraint
###

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.EqualTo{Float64},
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.Interval{Float64},
            MOI.ZeroOne,
            MOI.Integer,
        },
    },
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{<:Union{MOI.SOS1{Float64},MOI.SOS2{Float64}}},
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{
        <:Union{
            MOI.EqualTo{Float64},
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.Interval{Float64},
        },
    },
)
    return true
end

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

###
### Variable starting values
###

function MOI.supports(
    ::Optimizer,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
    ::Nothing,
)
    delete!(model.variable_start, x)
    return
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
    value,
)
    model.variable_start[x] = value
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
)
    return get(model.variable_start, x, nothing)
end

###
### Optimize and post-optimize functions
###

function MOI.optimize!(model::Optimizer)
    if !isempty(model.variable_start)
        columns = fill(Cint(0), length(model.variable_start))
        values = fill(0.0, length(model.variable_start))
        for (i, (variable, value)) in enumerate(model.variable_start)
            columns[i] = Cint(variable.value - 1)
            values[i] = value
        end
        Cbc_setMIPStartI(model, length(columns), columns, values)
    end
    t = time()
    model.variable_primal = nothing
    model.constraint_primal = nothing
    model.termination_status = Cbc_solve(model)
    model.has_solution = _result_count(model)
    model.solve_time = time() - t
    return
end

MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return Cbc_getNumCols(model)
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    if Cbc_getNumIntegers(model) == 0
        return MOI.get(model, MOI.ObjectiveValue())
    end
    return Cbc_getBestPossibleObjValue(model) + model.objective_constant
end

function MOI.get(model::Optimizer, ::MOI.NodeCount)::Int64
    return Cbc_getNodeCount(model)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return Cbc_getObjValue(model) + model.objective_constant
end

# Cbc does not provide a native way of accessing the relative gap,
# use the Gurobi convention instead.
function MOI.get(model::Optimizer, ::MOI.RelativeGap)
    incumbent = MOI.get(model, MOI.ObjectiveValue())
    bound = MOI.get(model, MOI.ObjectiveBound())
    gap = abs(bound - incumbent) / abs(incumbent)
    return isnan(gap) ? Inf : gap
end

_update_cache(::Optimizer, data::Vector{Float64}, ::Any, ::Any) = data

function _update_cache(model::Optimizer, ::Nothing, f_p::F, f_n::G) where {F,G}
    p = f_p(model)
    n = f_n(model)
    if p == C_NULL
        return fill(NaN, n)
    end
    return unsafe_wrap(Array, p, (n,))
end

_get_cached_solution(data::Vector{Float64}, x) = data[x.value]

function _get_cached_solution(data::Vector{Float64}, x::Vector)
    return [data[xi.value] for xi in x]
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::Union{MOI.VariableIndex,Vector{MOI.VariableIndex}},
)
    MOI.check_result_index_bounds(model, attr)
    model.variable_primal = _update_cache(
        model,
        model.variable_primal,
        Cbc_getColSolution,
        Cbc_getNumCols,
    )
    return _get_cached_solution(model.variable_primal, x)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:Any},
)
    MOI.check_result_index_bounds(model, attr)
    model.constraint_primal = _update_cache(
        model,
        model.constraint_primal,
        Cbc_getRowActivity,
        Cbc_getNumRows,
    )
    return _get_cached_solution(model.constraint_primal, c)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    index::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(index.value))
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    sense = Cbc_getObjSense(model)
    if sense == 1.0
        return MOI.MIN_SENSE
    elseif sense == -1.0
        return MOI.MAX_SENSE
    else
        @assert sense == 0.0
        return MOI.FEASIBILITY_SENSE
    end
end

struct Status <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::Status) = true

struct SecondaryStatus <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::SecondaryStatus) = true

MOI.get(model::Optimizer, ::Status) = Cbc_status(model)
MOI.get(model::Optimizer, ::SecondaryStatus) = Cbc_secondaryStatus(model)

const _STATUS = Dict{Cint,String}(
    Cint(-1) => "before branchAndBound",
    Cint(0) =>
        "finished - check isProvenOptimal or isProvenInfeasible to see if " *
        "solution found (or check value of best solution)",
    Cint(1) => "stopped - on maxnodes, maxsols, maxtime",
    Cint(2) => "execution abandoned due to numerical dificulties",
    Cint(5) => "user programmed interruption",
)

const _SECONDARY_STATUS = Dict{Cint,String}(
    Cint(-1) => "unset (status_ will also be -1)",
    Cint(0) => "search completed with solution",
    Cint(1) => "linear relaxation not feasible (or worse than cutoff)",
    Cint(2) => "stopped on gap",
    Cint(3) => "stopped on nodes",
    Cint(4) => "stopped on time",
    Cint(5) => "stopped on user event",
    Cint(6) => "stopped on solutions",
    Cint(7) => "linear relaxation unbounded",
    Cint(8) => "stopped on iteration limit",
)

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return """
    Cbc_status          = $(_STATUS[model.termination_status])
    Cbc_secondaryStatus = $(_SECONDARY_STATUS[Cbc_secondaryStatus(model)])
    """
end

MOI.get(model::Optimizer, ::MOI.ResultCount) = model.has_solution ? 1 : 0

function _result_count(model::Optimizer)
    if Cbc_getNumIntegers(model) == 0
        # Cbc forwards the solve to the LP solver if there are no integers, so
        # check the termination status for the result count.
        return model.termination_status == 0 ? 1 : 0
    end
    return Cbc_numberSavedSolutions(model) > 0 ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    status = model.termination_status
    if status == -1
        return MOI.OPTIMIZE_NOT_CALLED
    elseif Cbc_isProvenOptimal(model) != 0
        return MOI.OPTIMAL
    elseif Cbc_isProvenInfeasible(model) != 0
        if Cbc_getNumIntegers(model) == 0
            # Why Cbc. For LPs, this could mean dual infeasible.
            return MOI.INFEASIBLE_OR_UNBOUNDED
        else
            return MOI.INFEASIBLE
        end
    elseif Cbc_isContinuousUnbounded(model) != 0
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif Cbc_isNodeLimitReached(model) != 0
        return MOI.NODE_LIMIT
    elseif Cbc_isSecondsLimitReached(model) != 0
        return MOI.TIME_LIMIT
    elseif Cbc_isSolutionLimitReached(model) != 0
        return MOI.SOLUTION_LIMIT
    elseif status == 0
        return MOI.OTHER_ERROR
    elseif status == 1
        return MOI.OTHER_LIMIT
    elseif status == 2
        return MOI.NUMERICAL_ERROR
    else
        @assert status == 5
        return MOI.INTERRUPTED
    end
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    elseif MOI.get(model, MOI.ResultCount()) == 1
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function MOI.get(::Optimizer, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end


text:
module HyperNEAT

using IniFile

export CPPN

type Configuration
  prob_synapse_change_weight::Float64 # 1 - prob_synapse_change_weight = new random value
  prob_synapse_new::Float64
  prob_neuron_new::Float64
  prob_neuron_change_function::Float64
  #= prob_neuron_change_parameter::Float64 =#
  prob_deactivate_synapse::Float64
  prob_activate_synapse::Float64
  value_new_synapse::Float64
  delta_weight::Float64
  #= delta_neuron_parameter::Float64 =#
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  r::Float64 # mathing percentage
  threshold::Float64 # speciation threshold
  innovation::Int64 # global innovation number
  id::Int64
  life_time::Int64
  max_individuals::Int64
  min_individuals::Int64
  yars_port::Int64
  log_dir::String
  generations::Int64
  population_size::Int64
  use_gui::Bool
  use_single_gui::Bool
  use_plot::Bool
  xml::String
  ports::Vector{Int64}
  wd::String
  structure::String
  normalise_population_size::Bool
  max_weight::Float64
end

function read_cfg(filename::String)
  ini = Inifile()
  read(ini, filename)
  pscw = float(get(ini, "Mutation",     "probability synapse change weight",   -1.0))
  psn  = float(get(ini, "Mutation",     "probability synapse new",             -1.0))
  pnn  = float(get(ini, "Mutation",     "probability neuron new",              -1.0))
  pnc  = float(get(ini, "Mutation",     "probability neuron change function",  -1.0))
  #= pcp  = float(get(ini, "Mutation",     "probability neuron change parameter", -1.0)) =#
  pds  = float(get(ini, "Mutation",     "probability deactivate synapse",      -1.0))
  pas  = float(get(ini, "Mutation",     "probability activate synapse",        -1.0))
  vns  = float(get(ini, "Mutation",     "value new synapse",                   -1.0))
  dw   = float(get(ini, "Mutation",     "delta weight",                        -1.0))
  #= dnp  = float(get(ini, "Mutation",     "neuron parameter delta",              -1.0)) =#
  msw  = float(get(ini, "Mutation",     "maximal synaptic weight",             10.0))
  c1   = float(get(ini, "Speciation",   "c1",                                  -1.0))
  c2   = float(get(ini, "Speciation",   "c2",                                  -1.0))
  c3   = float(get(ini, "Speciation",   "c3",                                  -1.0))
  c4   = float(get(ini, "Speciation",   "c4",                                  -1.0))
  th   = float(get(ini, "Speciation",   "threshold",                           -1.0))
  r    = float(get(ini, "Reproduction", "r",                                   -1.0))
  mai  = int(get(ini,   "Reproduction", "maximum number of individuals",        10))
  mii  = int(get(ini,   "Reproduction", "minimum number of individuals",        10))
  np   = (get(ini,      "Reproduction", "normalise population size", "false") == "true")
  lt   = int(get(ini,   "General",      "life time",                           -1.0))
  yp   = int(get(ini,   "General",      "yars port",                           -1.0))
  ld   = get(ini,       "General",      "log directory",                       ".")
  g    = int(get(ini,   "General",      "generations",                         -1))
  ps   = int(get(ini,   "General",      "population size",                     -1))
  p    = int(get(ini,   "General",      "port",                                []))
  ug   = (get(ini,      "General",      "use gui",  "false")          == "true")
  usg  = (get(ini,      "General",      "use only one gui",  "false") == "true")
  up   = (get(ini,      "General",      "use plot", "false")          == "true")
  xml  = get(ini,       "General",      "xml",                    "")
  wd   = get(ini,       "General",      "working directory",      "")
  st   = get(ini,       "General",      "structure", "neat")

  #ffc  = get(ini, "fitness function coefficients", "c",                                  "-1.0")
  cfg = Configuration(
  pscw,
  psn,
  pnn,
  pnc,
  #= pcp, =#
  pds,
  pas,
  vns,
  dw,
  #= dnp, =#
  c1,
  c2,
  c3,
  c4,
  r,
  th,
  10, # innovation number
  0, # global id
  lt,
  mai,
  mii,
  yp,
  ld,
  g,
  ps,
  ug,
  usg,
  up,
  xml,
  [], # todo ports
  wd,
  st,
  np,
  msw
  )

  return cfg
end

#= function Base.show(io::IO, cfg::Configuration) =#
#= println("[Mutation]") =#
#= println("probability synapse change weight  = ", cfg.prob_synapse_change_weight) =#
#= println("probability synapse new            = ", cfg.prob_synapse_new) =#
#= println("probability neuron new             = ", cfg.prob_neuron_new) =#
#= println("probability neuron change function = ", cfg.prob_neuron_change_function) =#
#= println("probability deactivate synapse     = ", cfg.prob_deactivate_synapse) =#
#= println("probability activate synapse       = ", cfg.prob_activate_synapse) =#
#= println("value new synapse                  = ", cfg.value_new_synapse) =#
#= println("delta weight                       = ", cfg.delta_weight) =#
#= println("") =#
#= println("[Speciation]") =#
#= println("c1                                 = ", cfg.c1) =#
#= println("c2                                 = ", cfg.c2) =#
#= println("c3                                 = ", cfg.c3) =#
#= println("c4                                 = ", cfg.c4) =#
#= println("threshold                          = ", cfg.threshold) =#
#= println("") =#
#= println("[Reproduction]") =#
#= println("r                                  = ", cfg.r) =#
#= println("") =#
#= end =#

function write_cfg(filename::ASCIIString, cfg::Configuration)
  ini = Inifile()
  ini.sections["Mutation"]     = IniFile.HTSS()
  ini.sections["Speciation"]   = IniFile.HTSS()
  ini.sections["Reproduction"] = IniFile.HTSS()

  sec                                       = ini.sections["Mutation"]
  sec["probability synapse change weight"]  = "$(cfg.prob_synapse_change_weight)"
  sec["probability synapse new"]            = "$(cfg.prob_synapse_new)"
  sec["probability neuron new"]             = "$(cfg.prob_neuron_new)"
  sec["probability neuron change function"] = "$(cfg.prob_neuron_change_function)"
  sec["probability deactivate synapse"]     = "$(cfg.prob_deactivate_synapse)"
  sec["probability activate synapse"]       = "$(cfg.prob_activate_synapse)"
  sec["value new synapse"]                  = "$(cfg.value_new_synapse)"
  sec["delta weight"]                       = "$(cfg.delta_weight)"

  sec                                       = ini.sections["Speciation"]
  sec["c1"]                                 = "$(cfg.c1)"
  sec["c2"]                                 = "$(cfg.c2)"
  sec["c3"]                                 = "$(cfg.c3)"
  sec["c4"]                                 = "$(cfg.c4)"
  sec["threshold"]                          = "$(cfg.threshold)"
  sec                                       = ini.sections["Reproduction"]
  sec["r"]                                  = "$(cfg.r)"

  fd = open(filename, "w")
  write(fd, ini)
  close(fd)
end

#c = read_cfg("cfg.ini")
#write_cfg("/Users/zahedi/Desktop/test.ini", c)

type HyperNEATError <: Exception
  var::String
end

type NeuronGene
  innovation::Int64
  ntype::Symbol # 1 = Input Neuron, 2 = Output Neuron, 3 = Hidden Neuron
  func::Symbol
  #= parameters::Vector{Float64} =#
end

#= new_neuron(;innovation=-1, ntype=:input, func=:id, parameters = [1.0, 0.0]) = NeuronGene(innovation, ntype, func, parameters) =#

new_neuron(;innovation=-1, ntype=:input, func=:id) = NeuronGene(innovation, ntype, func)

type SynapseGene
  src::Int64
  dest::Int64
  weight::Float64
  innovation::Int64
  active::Bool
end

new_synapse(;src=-1, dest=-1, weight=0.0, innovation=-1, active=true) = SynapseGene(src, dest, weight, innovation, active)

type Genome
  neurons::Vector{NeuronGene}
  synapses::Vector{SynapseGene}
end

new_genome(;neurons=[], synapses=[]) = Genome(neurons, synapses)

type Individual
  genomes::Vector{Genome}
  fitness::Float64
  id::Int64
  parents::(Int64, Int64)
  species_id::Int64
  age::Int64
  stats::Vector{Float64}
end

new_individual(;genomes=[], fitness=0.0, id=-1, parents=(-1,-1), species_id = -1, age = 0, stats = []) = Individual(genomes, fitness, id, parents, species_id, age, stats)

type Species
  individuals::Vector{Individual}
  size::Int64
end

new_species(;individuals=[], size=-1) = Species(individuals, size)

type Population
  species::Vector{Species}
  generation::Int64
end

new_population(;species=[], generation=-1) = Population(species, generation)

function population_write(filename::String, p::Population)
  f = open(filename, "w")
  serialize(f, p)
  close(f)
end

function population_read(filename::String)
  f = open(filename, "r")
  p = deserialize(f)
  close(f)
  return p
end

function mutation_add_neuron!(genome::Genome, cfg::Configuration)
  if length(genome.synapses) == 0
    return
  end
  if rand() < cfg.prob_neuron_new
    active_synapses = filter(s->s.active, genome.synapses)
    synapse         = active_synapses[ceil(rand() * length(active_synapses))]
    synapse.active  = false
    function_keys   = collect(keys(functions))
    func            = function_keys[ceil(rand() * length(function_keys))]
    cfg.innovation  = cfg.innovation + 1
    #= genome.neurons  = [genome.neurons, NeuronGene(cfg.innovation, :hidden, func, [1.0, 0.0])] =#
    genome.neurons  = [genome.neurons, NeuronGene(cfg.innovation, :hidden, func)]
    neuron_index    = length(genome.neurons)
    cfg.innovation  = cfg.innovation + 1
    genome.synapses = [genome.synapses, SynapseGene(neuron_index, synapse.dest, synapse.weight, cfg.innovation, true)]
    cfg.innovation  = cfg.innovation + 1
    genome.synapses = [genome.synapses, SynapseGene(synapse.src,  neuron_index, 1.0,            cfg.innovation, true)]
  end
end

function mutation_change_neuron!(genome::Genome, cfg::Configuration)
  if rand() < cfg.prob_neuron_change_function
    neurons       = filter(n->(n.ntype == :hidden), genome.neurons)
    if length(neurons) == 0
      return
    end
    neuron        = neurons[int(ceil(rand() * length(neurons)))]
    function_keys = collect(keys(functions))
    func          = neuron.func
    while func == neuron.func
      func          = function_keys[ceil(rand() * length(function_keys))]
    end
    neuron.func   = func
  end
end

#= function mutation_change_neuron_parameter!(genome::Genome, cfg::Configuration) =#
  #= if rand() < cfg.prob_neuron_change_parameter =#
    #= neurons           = filter(n->(n.ntype == :hidden || n.ntype == :output), genome.neurons) =#
    #= neuron            = neurons[int(ceil(rand() * length(neurons)))] =#
    #= neuron.parameters = map(p->p + (2.0 * rand() - 1.0) * cfg.delta_neuron_parameter, neuron.parameters) =#
  #= end =#
#= end =#

function list_of_possible_synapses(genome::Genome)
  last_input_neuron       = length(filter(n->n.ntype == :input, genome.neurons))
  nr_of_non_input_neurons = length(filter(n->n.ntype != :input, genome.neurons))

  all_possible_synapses = []
  for src = 1:length(genome.neurons)
    for dest = last_input_neuron+1:last_input_neuron+nr_of_non_input_neurons
      if length(filter(s->(s.src == src && s.dest == dest), genome.synapses)) == 0
        all_possible_synapses = [all_possible_synapses, (src, dest)]
      end
    end
  end
  all_possible_synapses
end

function mutation_add_synapse!(genome::Genome, cfg::Configuration)
  if rand() < cfg.prob_synapse_new
    l = list_of_possible_synapses(genome)
    if length(l) == 0
      return
    end
    cfg.innovation = cfg.innovation + 1
    (src,dest) = l[int(ceil(rand() * length(l)))]
    genome.synapses = [genome.synapses, SynapseGene(src, dest, (2.0 * rand() - 1.0) * cfg.value_new_synapse, cfg.innovation, true)]
  end
end

function mutation_change_synapse!(genome::Genome, cfg::Configuration)
  if length(genome.synapses) == 0
    return
  end
  synapse_index = int(ceil(rand() * length(genome.synapses)))
  synapse = genome.synapses[synapse_index]
  if rand() < cfg.prob_synapse_change_weight
    synapse.weight = synapse.weight + (2.0 * rand() - 1.0) * cfg.delta_weight
  else
    synapse.weight =                  (2.0 * rand() - 1.0) * cfg.value_new_synapse
  end

  genome.synapses[synapse_index] = synapse
end

function mutation_deactivate_synapse!(genome::Genome, cfg::Configuration)
  if length(genome.synapses) == 0
    return
  end
  if rand() < cfg.prob_deactivate_synapse
    if length(genome.synapses) != sum(map(s->s.active?0:1, genome.synapses))
      return # there no active synapses
    end

    active_synapses = filter(s->s.active, genome.synapse)
    active_synapses[ceil(rand() * length(active_synapses))].active = false

  end
end

function mutation_activate_synapse!(genome::Genome, cfg::Configuration)
  if length(genome.synapses) == 0
    return
  end
  if rand() < cfg.prob_activate_synapse
    if length(genome.synapses) != sum(map(s->s.active?1:0, genome.synapses))
      return # there no deactivated synapses
    end
    deactivated_synapses = filter(s->s.active == false, genome.synapses)
    if length(deactivated_synapses) == 0
      return
    end
    deactivated_synapses[ceil(rand() * length(deactivated_synapses))].active = true
  end
end

function crossover(mother::Genome, father::Genome)
  child = new_genome()
  max_innovation_number = 1
  if length(mother.synapses) > 0 && length(father.synapses) > 0
    max_innovation_number = maximum([map(s->s.innovation, mother.synapses),map(s->s.innovation, father.synapses)])
  elseif length(mother.synapses) > 0
    max_innovation_number = maximum([map(s->s.innovation, mother.synapses)])
  elseif length(father.synapses) > 0
    max_innovation_number = maximum([map(s->s.innovation, father.synapses)])
  end
  for i = 1:max_innovation_number
    # synapses
    ms = filter(s->s.innovation == i, mother.synapses)
    fs = filter(s->s.innovation == i, father.synapses)
    cs = []
    if length(ms) > 0 && length(fs) > 0 && length(ms) != length(fs)
      throw(HyperNEATError("crossover: Number of synapses don't match $(length(ms)) vs $(length(fs)) for innovation $i"))
    end
    if length(ms) > 0 && length(fs) > 0
      cs = [rand() < 0.5?ms[i]:fs[i] for i=1:length(ms)]
    elseif length(ms) > 0
      cs = ms
    else
      cs = fs
    end

    for s in cs
      child.synapses = [child.synapses, SynapseGene(s.src, s.dest, s.weight, s.innovation, s.active)]
    end

    # neurons
    mn = filter(n->n.innovation == i, mother.neurons)
    fn = filter(n->n.innovation == i, father.neurons)
    cn = []
    if (length(mn) > 0 && length(fn) > 0) && (length(mn) != length(fn))
      throw(HyperNEATError("crossover: Number of neurons don't match: father: $(length(fn)) vs mother: $(length(mn))"))
    end
    if length(mn) > 0 && length(fn) > 0
      cn = [rand() < 0.5? mn[i] : fn[i] for i = 1:length(mn)]
    elseif length(mn) > 0
      cn = mn
    else
      cn = fn
    end

    for n in cn
      #= child.neurons = [child.neurons, NeuronGene(n.innovation, n.ntype, n.func, n.parameters)] =#
      child.neurons = [child.neurons, NeuronGene(n.innovation, n.ntype, n.func)]
    end
  end
  child
end

function speciation_distance(g1::Genome, g2::Genome, cfg::Configuration)
  N                        = float(maximum([length(g1.synapses), length(g2.synapses)]))
  max_innovation_number_g1 = length(g1.synapses)>0?maximum(map(s->s.innovation, g1.synapses)):0
  max_innovation_number_g2 = length(g2.synapses)>0?maximum(map(s->s.innovation, g2.synapses)):0
  excess                   = float(abs(max_innovation_number_g1 - max_innovation_number_g2))
  disjoint                 = 0.0
  matching                 = 0.0
  weight_difference        = 0.0
  for i = 1:maximum([max_innovation_number_g1, max_innovation_number_g2])
    found_g1 = filter(s->s.innovation == i, g1.synapses)
    found_g2 = filter(s->s.innovation == i, g2.synapses)
    if length(found_g1) != length(found_g2)
      disjoint = disjoint + 1.0
    elseif length(found_g1) > 0
      weight_difference = weight_difference + abs(found_g1[1].weight - found_g2[1].weight)
      matching = matching + 1.0
    end
  end

  neuron_difference = 0.0
  n_matching        = 0.0
  for i = 1:maximum([max_innovation_number_g1, max_innovation_number_g2])
    found_g1 = filter(n->n.innovation == i, g1.neurons)
    found_g2 = filter(n->n.innovation == i, g2.neurons)
    if length(found_g1) > 0 && length(found_g2) > 0
      if found_g1[1].func != found_g2[1].func
        neuron_difference = neuron_difference + 1.0
        n_matching        = n_matching + 1.0
      end
    else
      neuron_difference = neuron_difference + 1.0
      n_matching        = n_matching + 1.0
    end
  end
  a = (N>0)?((cfg.c1 * excess + cfg.c2 * disjoint) / N):0
  b = (matching>0)?(cfg.c3 * weight_difference / matching):0
  c = (n_matching>0)?(cfg.c4 * neuron_difference / n_matching):0

  a + b + c
end

function calculate_species_sizes!(population::Population, cfg::Configuration)

  fitness_values = []
  for s in population.species
    for individual in s.individuals
      fitness_values = [fitness_values, individual.fitness]
    end
  end

  minimal_fitness = minimum(fitness_values)

  fitness_values = []
  for s in population.species
    l = float(length(s.individuals))
    for individual in s.individuals
      individual.fitness = (individual.fitness + abs(minimal_fitness)) / l
      fitness_values = [fitness_values, individual.fitness]
    end
  end

  average_fitness = mean(fitness_values)

  for s in population.species # new size depends
    ss     = float(sum(map(i->i.fitness, s.individuals)))
    av     = float(average_fitness)
    s.size = maximum([cfg.max_individuals, int(round(sum(map(i->i.fitness, s.individuals)) / average_fitness))])
  end

  for s in population.species
    if s.size == 1
      s.size = cfg.min_individuals
    end
  end
end

function mutation_scale_weights!(genome::Genome, cfg::Configuration)
  for s in genome.synapses
    if abs(s.weight) > cfg.max_weight
      s.weight = (s.weight<0:-1:1) * cfg.max_weight
    end
  end
end

function mutate!(g::Genome, cfg::Configuration)
  mutation_deactivate_synapse!(g, cfg)
  mutation_add_neuron!(g, cfg)
  #= mutation_change_neuron_parameter!(g, cfg) =#
  mutation_change_synapse!(g, cfg)
  mutation_add_synapse!(g, cfg)
  mutation_change_neuron!(g, cfg)
  mutation_activate_synapse!(g, cfg)
  mutation_scale_weights!(g, cfg)
  g
end

function reproduce!(population::Population, cfg::Configuration)
  for s in population.species
    s.individuals = sort(s.individuals, by=x->x.fitness, rev=true)
  end

  calculate_species_sizes!(population, cfg)

  if cfg.normalise_population_size == true
    sum_sizes = sum(map(s->s.size, population.species))

    for s in population.species
      s.size = int(round(float(s.size) / float(sum_sizes) * cfg.max_individuals))
    end

    if length(population.species) == 1
      if population.species[1].size != cfg.max_individuals
        population.species[1].size = cfg.max_individuals
      end
    end
  else
    for s in population.species
      if s.size > cfg.max_individuals
        s.size = cfg.max_individuals
      elseif s.size < cfg.min_individuals
        s.size = cfg.min_individuals
      end
    end

    #= sum_sizes = sum(map(s->s.size, population.species)) =#
    #= if sum_sizes < cfg.min_individuals =#
      #= factor = float64(cfg.min_individuals / sum_sizes) =#
      #= for s in population.species =#
        #= s.size = int(ceil(s.size * factor)) =#
      #= end =#
    #= end =#
  end

  individuals = []

  for s in population.species
    N               = ceil(cfg.r*length(s.individuals))
    parents         = s.individuals[1:N]
    s.individuals   = filter(p->p.age < 5, parents)
    nr_of_offspring = s.size - length(s.individuals)
    for i in s.individuals
      i.age = i.age + 1
    end
    if s.size > length(s.individuals)
      for i = 1:nr_of_offspring
        mother        = parents[ceil(rand() * N)]
        father        = parents[ceil(rand() * N)]

        nr_of_genomes = length(mother.genomes)

        child         = new_individual()
        child.genomes = map(i->crossover(mother.genomes[i], father.genomes[i]), [1:nr_of_genomes])
        child.parents = (mother.id, father.id)
        child.id      = cfg.id
        cfg.id        = cfg.id + 1

        for g in child.genomes
          mutate!(g, cfg)
        end
        individuals = [individuals, child]
      end
    else
      s.size = 0
    end
  end

  population.generation = population.generation + 1

  population.species = filter(s->s.size > 0, population.species)

  first = new_species(individuals=[individuals[1]])
  population.species = [first]

  for i in individuals[2:end]
    found = false
    for s in population.species
      random_individual = s.individuals[int(ceil(rand() * length(s.individuals)))]
      dist = maximum([speciation_distance(i.genomes[j], random_individual.genomes[j], cfg) for j=1:length(i.genomes)])
      if dist < cfg.threshold
        s.individuals = [s.individuals, i]
        found = true
        break
      end
    end
    if found == false
      population.species = [population.species, new_species(individuals=[i])]
    end
  end
  population
end

type Neuron
  inputs::Vector{Neuron}
  weights::Vector{Float64}
  func::Function
  output::Float64
  #= parameters::Vector{Float64} =#
end

type CPPN
  inputs::Vector{Neuron}
  hidden::Vector{Neuron}
  outputs::Vector{Neuron}
  neurons::Vector{Neuron}
end

function update_neuron!(n)
  act = sum([n.inputs[i].output * n.weights[i] for i=1:length(n.inputs)])
  if act == NaN
    println("Inputs: ",  [n.inputs[i].output for i=1:length(n.inputs)])
    println("Weights: ", [n.weights[i]       for i=1:length(n.inputs)])
  end
  if abs(act) > 1000
    act = (act<0)?-1000.0:1000.0
  end
  #= println(n.parameters) =#
  #= n.output = n.func(n.parameters, act) =#
  n.output = n.func(act)
end

function update_cppn!(cppn::CPPN, inputs::Vector{Float64})
  for i in 1:minimum([length(cppn.inputs), length(inputs)])
    cppn.inputs[i].output = inputs[i]
  end

  for i = 1:10
    for n in cppn.hidden
      update_neuron!(n)
    end
    for n in cppn.outputs
      update_neuron!(n)
    end
  end
  map(o->o.output, cppn.outputs)
end

create_cppn()     = CPPN([], [], [])
gauss(x,x0,sigma) = exp(-(x-x0)^2 / (2.0 * sigma^2))
sigmoid(x)        = 1.0 / (1.0 + exp(-x))
id                = x->x

#= functions = { =#
#= :gauss => (p::Vector{Float64}, v::Float64)->gauss(v, p[2], p[1]), =#
#= :sin   => (p::Vector{Float64}, v::Float64)->sin(     p[1] * v + p[2]), =#
#= :cos   => (p::Vector{Float64}, v::Float64)->cos(     p[1] * v + p[2]), =#
#= :ncos  => (p::Vector{Float64}, v::Float64)->(-cos(   p[1] * v + p[2])), =#
#= :tanh  => (p::Vector{Float64}, v::Float64)->tanh(    p[1] * v + p[2]), =#
#= :sigm  => (p::Vector{Float64}, v::Float64)->sigmoid( p[1] * v + p[2]), =#
#= :id    => (p::Vector{Float64}, v::Float64)->id(      p[1] * v + p[2])} =#

#= function add_input_neuron!(cppn::CPPN,  f::Symbol, parameters::Vector{Float64}) =#
  #= n = Neuron([], [], functions[f], 0.0, parameters) =#
  #= cppn.inputs  = [cppn.inputs, n] =#
  #= cppn.neurons = [cppn.neurons, n] =#
#= end =#

#= function add_output_neuron!(cppn::CPPN, f::Symbol, parameters::Vector{Float64}) =#
  #= n = Neuron([], [], functions[f], 0.0, parameters) =#
  #= cppn.outputs = [cppn.outputs, n] =#
  #= cppn.neurons = [cppn.neurons, n] =#
#= end =#

#= function add_hidden_neuron!(cppn::CPPN, f::Symbol, parameters::Vector{Float64}) =#
  #= n = Neuron([], [], functions[f], 0.0, parameters) =#
  #= cppn.hidden  = [cppn.hidden,  n] =#
  #= cppn.neurons = [cppn.neurons, n] =#
#= end =#

functions = {
:gauss => (v::Float64)->gauss(v, 0.0, 1.0),
:sin   => (v::Float64)->  sin(v),
:nsin  => (v::Float64)->(-sin(v)),
:cos   => (v::Float64)->  cos(v),
:ncos  => (v::Float64)->(-cos(v)),
:tanh  => (v::Float64)->tanh(v),
:sigm  => (v::Float64)->sigmoid(v),
:id    => (v::Float64)->id(v)}

function add_input_neuron!(cppn::CPPN,  f::Symbol)
  n = Neuron([], [], functions[f], 0.0)
  cppn.inputs  = [cppn.inputs, n]
  cppn.neurons = [cppn.neurons, n]
end

function add_output_neuron!(cppn::CPPN, f::Symbol)
  n = Neuron([], [], functions[f], 0.0)
  cppn.outputs = [cppn.outputs, n]
  cppn.neurons = [cppn.neurons, n]
end

function add_hidden_neuron!(cppn::CPPN, f::Symbol)
  n = Neuron([], [], functions[f], 0.0)
  cppn.hidden  = [cppn.hidden,  n]
  cppn.neurons = [cppn.neurons, n]
end

function add_synapse!(cppn::CPPN, src::Neuron, dest::Neuron, weight::Float64)
  dest.inputs  = [dest.inputs,  src]
  dest.weights = [dest.weights, weight]
end

function phenotype(genome::Genome)
  cppn = CPPN([], [], [], [])
  for n in genome.neurons
    if n.ntype == :input
      #= add_input_neuron!(cppn, n.func, n.parameters) =#
      add_input_neuron!(cppn, n.func)
    elseif n.ntype == :output
      #= add_output_neuron!(cppn, n.func, n.parameters) =#
      add_output_neuron!(cppn, n.func)
    elseif n.ntype == :hidden
      #= add_hidden_neuron!(cppn, n.func, n.parameters) =#
      add_hidden_neuron!(cppn, n.func)
    else
      throw(HyperNEATError("cppn.jl: unknown neuron type given in function phenotype(...)"))
    end
  end
  for s in genome.synapses
    if s.active == true
      add_synapse!(cppn, cppn.neurons[s.src], cppn.neurons[s.dest], s.weight)
    end
  end
  cppn
end

end


text:
Maybe(T::Type) = Union{T,Nothing}

###################
### DATASETTYPE ###
###################

abstract type DatasetType end

### PonctualDataset ###
struct PonctualDataset <: DatasetType
    data::DataFrame
    target::Vector{Symbol}
end
function PonctualDataset(data, target ; keys::Vector{tK}=Symbol[]) where tK<:Union{String,Symbol}
    check_keys(names(data),keys)
    PonctualDataset(data, target)
end
PonctualDataset(data, target::Symbol) = PonctualDataset(data, [target])

### TimeseriesDataset ###
struct TimeseriesDataset <: DatasetType
    data::DataFrame
    target::Vector{Symbol}
    var::Symbol
end
function TimeseriesDataset(data, target, var ; keys::Vector{tK}=Symbol[]) where tK<:Union{String,Symbol}
    check_keys(names(data),keys)
    TimeseriesDataset(data, target, var)
end
TimeseriesDataset(data, target::Symbol, var) = TimeseriesDataset(data, [target], var)

# helper functions :
function group(ds::TimeseriesDataset)
    groupby(ds.data,ds.var)
end

function get_timeseries(tsds::TimeseriesDataset,val)
    df = tsds.data
    return TimeseriesDataset(df[df[:,tsds.var] .== val,:], tsds.target, tsds.var)
end

# ALL DATASETS :
Base.@kwdef struct Dataset
    ponctual::Vector{PonctualDataset} = PonctualDataset[]
    timeseries::Vector{TimeseriesDataset} = TimeseriesDataset[]
end
Dataset(tsd::Vector{TimeseriesDataset},pd::Vector{PonctualDataset}) = Dataset(pd,tsd)
Dataset(pd::PonctualDataset,tsd::TimeseriesDataset) = Dataset([pd],[tsd])
Dataset(pd::Vector{PonctualDataset},tsd::TimeseriesDataset) = Dataset(pd,[tsd])
Dataset(pd::PonctualDataset,tsd::Vector{TimeseriesDataset}) = Dataset([pd],tsd)

function check_keys(colnames::Vector{String},keys::Vector{tK}) where tK<:Union{String,Symbol}
    if !isempty(keys)
        isin = colnames .∈ Ref(string.(keys))
        id_notin = findall(.!isin)
        !isempty(id_notin) && throw(ArgumentError("column names $(colnames[id_notin]) are in keys argument"))
    end
end

export DatasetType, PonctualDataset, TimeseriesDataset, Dataset, group, get_timeseries

max_stars_repo_name: Leooop/DataFormatter

max_issues_repo_name: Leooop/DataFormatter

max_forks_repo_name: Leooop/DataFormatter


text:
"""
    struct FEMBasis

Provides a Julia-based implementation of a finite element basis, but using HelFEM-provided
primitive polynomials. Unlike [`RadialBasis`](@ref), it does not scale the basis functions
with ``1/r``, and is therefore valid for any ``x`` interval on ``\\mathbb{R}``.

# Constructors

    FEMBasis(pb::PolynomialBasis, boundaries; [nquad])

Constructs a FEM basis with elements defined by the boundaries specified by `boundaries`
using the set of polynomials defined by `pb` as the primitive polynomials basis. Optionally,
`nquad` can be specified to change the number of quadrature points used in the Gauss-Chebyshev
quadrature when evaluating matrix elements.
"""
struct FEMBasis
    pb :: PolynomialBasis
    boundaries :: Vector{Float64}
    # Option whether to impose Dirichet 0 boundary conditions?
    bcdropleft :: Bool
    bcdropright :: Bool
    # Cache for numerical quadrature
    qxs :: Vector{Float64} # quadrature points on [-1, 1]
    qws :: Vector{Float64} # quadrature weights corresponding to qxs
    # Caches of basis polynomials and their derivatives evaluated at quadrature points qxs
    # Resulting in length(qxs) × length(pb) matrices
    pbf :: Matrix{Float64}
    pbdf :: Matrix{Float64}

    function FEMBasis(
            pb::PolynomialBasis, boundaries::AbstractVector;
            nquad=5*length(pb),
            # TODO: Implement actual support for not dropping the edge functions
            # bcdropleft = true, bcdropright = true,
        )
        @assert length(boundaries) >= 2
        # Check and sort the boundaries
        @assert all(isfinite.(boundaries))
        boundaries_sorted = collect(boundaries)
        sort!(boundaries_sorted)
        # Quadrature and primitive basis functions
        qxs, qws = chebyshev(nquad)
        pbf, pbdf = pb_eval(pb, qxs)
        # Construct the FEMBasis object
        b = new(pb, boundaries_sorted, true, true, #=TODO: bcdropleft, bcdropright,=# qxs, qws, pbf, pbdf)
        return b
    end
end

function Base.show(io::IO, b::FEMBasis)
    pbname = primbas_name(b.pb)
    pb_nnodes = nnodes(b.pb)
    nelem = nelements(b)
    x1, xn = first(b.boundaries), last(b.boundaries)
    nbf = length(b)
    # TODO: add introspection for BCs
    write(io, "FEMBasis($pbname($pb_nnodes nodes), $nelem elements, on [$x1, $xn]) with $nbf basis function$(nbf>1 ? "s" : "")")
end

"""
    length(b::FEMBasis) -> Int

Return the number of basis functions represented by `b`.
"""
function Base.length(b::FEMBasis)
    nelem = nelements(b)
    pb_length = length(b.pb)
    nbf = (pb_length - 1)*nelem + 1
    # Imposing a boundary condition on either left or right end of the interval will means
    # removing the left- or right-most function which is non-zero at the boundary.
    b.bcdropleft && (nbf -= 1)
    b.bcdropright && (nbf -= 1)
    return nbf
end

"""
    boundaries(b::FEMBasis)

Return an array of element boundaries as a sorted read-only array.
"""
boundaries(b::FEMBasis) = ReadOnlyArray(b.boundaries)

"""
    nelements(b::FEMBasis)

Return the number of finite elements defining the basis.
"""
nelements(b::FEMBasis) = length(b.boundaries) - 1

"""
    HelFEM.nquad(b::FEMBasis)

Return the number of quadrature points used when evaluating the matrix elements.
"""
nquad(b::FEMBasis) = length(b.qxs)

"""
    (b::FEMBasis)(qs; derivative=0) -> Matrix

Functor syntax to evaluate the basis functions at specified points. Returns a matrix where
each column corresponds to a basis function and row corresponds to a point in the input
vector `qs`.
"""
function (b::FEMBasis)(qs; derivative=0)
    @assert derivative in [0, 1]
    @assert first(b.boundaries) <= minimum(qs)
    @assert last(b.boundaries) >= maximum(qs)
    bf = zeros(Float64, (length(qs), length(b)))
    pb_nbf = length(b.pb)
    # TODO: support _not_ dropping the edge functions
    for iel = 1:nelements(b)
        qmin, qmax = b.boundaries[iel], b.boundaries[iel+1]
        q0, qλ = (qmax + qmin) / 2, (qmax - qmin) / 2
        # We assign the r values on the boundaries to the element on the left. In other words,
        # we'll define the element to be r ∈ (rmin, rmax], except when iel == 1.
        idxs = (iel == 1) ? findall(q -> qmin <= q <= qmax, qs) : findall(q -> qmin < q <= qmax, qs)
        isempty(idxs) && continue
        # Scale the rs values in this element down to the [-1, 1] range
        xs = (qs[idxs] .- q0) ./ qλ
        # Evaluate basis functions (or their derivative, if dy is true) at xs
        ys = (derivative == 0) ? b.pb(xs) : pb_eval(b.pb, xs).dfs ./ qλ
        # Assign the calculated basis values to the correct place in the basis matrix
        bfrange_start = (iel == 1) ? 1 : (pb_nbf - 1) * (iel - 1)
        bfrange_end = (iel == nelements(b)) ? length(b) : (pb_nbf - 1) * iel
        ysrange_start = (iel == 1) ? 2 : 1
        ysrange_end = (iel == nelements(b)) ? (pb_nbf - 1) : pb_nbf
        # Assign them into the output matrix
        bf[idxs, bfrange_start:bfrange_end] .= ys[:, ysrange_start:ysrange_end]
    end
    return bf
end

"""
    radial_integral(b::FEMBasis, f; lderivative=false, rderivative=false) -> Matrix

Evaluate the matrix elements of the function `f` in the basis

```math
f_{ij} = \\int_{x_1}^{x_{N+1}} b^*_{i}(x) f(x) b_{j}(x) ~ dx
```

`f` is expected to be a callable object that takes a single argument.

If the keyword arguments `lderivative` and `rderivative`, the derivative of the corresponding
basis function is used instead.
"""
function radial_integral(b::FEMBasis, f; lderivative=false, rderivative=false)
    Fij = zeros(Float64, (length(b), length(b)))
    pb_nbf = length(b.pb)
    # TODO: support _not_ dropping the edge functions
    for iel = 1:nelements(b)
        qmin, qmax = b.boundaries[iel], b.boundaries[iel+1]
        q0, qλ = (qmax + qmin) / 2, (qmax - qmin) / 2
        qs = b.qxs .* qλ .+ q0 # scale the quadrature points to q coordinate in element
        # Evaluate the function at quadrature points in this element
        fs = f.(qs)
        # Left and right function values
        Bi = lderivative ? b.pbdf : b.pbf
        Bj = rderivative ? b.pbdf : b.pbf
        # The integrals get a qλ^(1 - n - m) factor, where n and m are the orders of the
        # derivatives of the left and right basis functions. The reason is because we scale
        # the x axis in each element (by qλ), which means the derivative of the basis
        # function is 1/qλ times the derivative of the primitive polynomial. The additional
        # 1 comes from the integral itself (dq = qλ * dx).
        qλfact = qλ^(1 - (lderivative ? 1 : 0) - (rderivative ? 1 : 0))
        # Operator matrix in the element
        Fij_el = Bi' * (b.qws .* fs .* Bj) .* qλfact
        # Assign the calculated basis values to the correct place in the basis matrix
        bfrange_start = (iel == 1) ? 1 : (pb_nbf - 1) * (iel - 1)
        bfrange_end = (iel == nelements(b)) ? length(b) : (pb_nbf - 1) * iel
        ysrange_start = (iel == 1) ? 2 : 1
        ysrange_end = (iel == nelements(b)) ? (pb_nbf - 1) : pb_nbf
        # Assign them into the output matrix
        bfrange, ysrange = bfrange_start:bfrange_end, ysrange_start:ysrange_end
        Fij[bfrange, bfrange] .+= Fij_el[ysrange, ysrange]
    end
    return Fij
end

"""
    elementrange(b::FEMBasis, k::Integer) -> (r_k, r_{k+1})

Returns a tuple with the start and end boundary of the `k`-th element.
"""
function elementrange(b::FEMBasis, k::Integer)
    element_boundaries = boundaries(b)
    1 <= k < length(element_boundaries) || throw(DomainError(k, "Element index k out of range (1 <= k <= $(nelements(b)))"))
    element_boundaries[k], element_boundaries[k+1]
end

"""
    scale_to_element(b::RadialBasis, k::Integer, xs)

Scales the ``x`` values within the element `k` to the corresponding ``r`` coordinates.
"""
function scale_to_element(b::FEMBasis, k::Integer, xs)
    rmin, rmax = elementrange(b, k)
    r0, rλ = (rmax + rmin) / 2, (rmax - rmin) / 2
    xs .* rλ .+ r0
end

function controlpoints(b::FEMBasis)
    # We only know how to calculate control points for the LIP basis currently
    b.pb.primbas == 4 || throw(ArgumentError("Can only calculate control points for LIP basis"))
    # The LIPBasis uses Gauss-Lobatto nodes, so the following only applies to LIP
    pb_nnodes = nnodes(b.pb)
    xs, _ = lobatto(pb_nnodes) # we only need the nodes, can discard the weights
    rs = Vector{Float64}(undef, (pb_nnodes - 1)*nelements(b) + 1)
    for k in 1:nelements(b)
        rs[(pb_nnodes-1)*(k-1)+1:(pb_nnodes-1)*k+1] .= scale_to_element(b, k, xs)
    end
    return rs
end

function functionoverlap(b::FEMBasis, f; derivative=false)
    fi = zeros(Float64, length(b))
    pb_nbf = length(b.pb)
    # TODO: support _not_ dropping the edge functions
    for iel = 1:nelements(b)
        qmin, qmax = b.boundaries[iel], b.boundaries[iel+1]
        q0, qλ = (qmax + qmin) / 2, (qmax - qmin) / 2
        qs = b.qxs .* qλ .+ q0 # scale the quadrature points to q coordinate in element
        # Evaluate the function at quadrature points in this element
        fs = f.(qs)
        # Left and right function values
        Bi = derivative ? b.pbdf : b.pbf
        # The integrals get a qλ^(1 - n - m) factor, where n and m are the orders of the
        # derivatives of the left and right basis functions. The reason is because we scale
        # the x axis in each element (by qλ), which means the derivative of the basis
        # function is 1/qλ times the derivative of the primitive polynomial. The additional
        # 1 comes from the integral itself (dq = qλ * dx).
        qλfact = qλ^(1 - (derivative ? 1 : 0))
        # Operator matrix in the element
        fi_el = sum(b.qws .* fs .* Bi, dims=1) .* qλfact
        # Assign the calculated basis values to the correct place in the basis matrix
        bfrange_start = (iel == 1) ? 1 : (pb_nbf - 1) * (iel - 1)
        bfrange_end = (iel == nelements(b)) ? length(b) : (pb_nbf - 1) * iel
        ysrange_start = (iel == 1) ? 2 : 1
        ysrange_end = (iel == nelements(b)) ? (pb_nbf - 1) : pb_nbf
        # Assign them into the output matrix
        bfrange, ysrange = bfrange_start:bfrange_end, ysrange_start:ysrange_end
        fi[bfrange] .+= fi_el[ysrange]
    end
    return fi
end


text:
# Date functions

### Core query functions

# Monday = 1....Sunday = 7
dayofweek(days) = mod1(days,7)

# Number of days in year
daysinyear(y) = 365 + isleapyear(y)

# Day of the year
const MONTHDAYS = [0,31,59,90,120,151,181,212,243,273,304,334]
dayofyear(y,m,d) = MONTHDAYS[m] + d + (m > 2 && isleapyear(y))

### Days of the Week
dayofweek(dt::TimeType) = dayofweek(days(dt))

const Monday,Tuesday,Wednesday,Thursday,Friday,Saturday,Sunday = 1,2,3,4,5,6,7
const Mon,Tue,Wed,Thu,Fri,Sat,Sun = 1,2,3,4,5,6,7
const english_daysofweek = [1=>"Monday",2=>"Tuesday",3=>"Wednesday",
                            4=>"Thursday",5=>"Friday",6=>"Saturday",7=>"Sunday"]
const VALUETODAYOFWEEK = (UTF8String=>Dict{Int,UTF8String})["english"=>english_daysofweek]
const english_daysofweekabbr = [1=>"Mon",2=>"Tue",3=>"Wed",
                                4=>"Thu",5=>"Fri",6=>"Sat",7=>"Sun"]
const VALUETODAYOFWEEKABBR = (UTF8String=>Dict{Int,UTF8String})["english"=>english_daysofweekabbr]
dayname(dt::Integer;locale::String="english") = VALUETODAYOFWEEK[locale][dt]
dayabbr(dt::Integer;locale::String="english") = VALUETODAYOFWEEKABBR[locale][dt]
dayname(dt::TimeType;locale::String="english") = VALUETODAYOFWEEK[locale][dayofweek(dt)]
dayabbr(dt::TimeType;locale::String="english") = VALUETODAYOFWEEKABBR[locale][dayofweek(dt)]

# Convenience methods for each day
ismonday(dt::TimeType) = dayofweek(dt) == Mon
istuesday(dt::TimeType) = dayofweek(dt) == Tue
iswednesday(dt::TimeType) = dayofweek(dt) == Wed
isthursday(dt::TimeType) = dayofweek(dt) == Thu
isfriday(dt::TimeType) = dayofweek(dt) == Fri
issaturday(dt::TimeType) = dayofweek(dt) == Sat
issunday(dt::TimeType) = dayofweek(dt) == Sun

# i.e. 1st Monday? 2nd Monday? 3rd Wednesday? 5th Sunday?
dayofweekofmonth(dt::TimeType) = (d = day(dt); return d < 8 ? 1 :
    d < 15 ? 2 : d < 22 ? 3 : d < 29 ? 4 : 5)

# Total number of a day of week in the month
# e.g. are there 4 or 5 Mondays in this month?
const TWENTYNINE = IntSet(1,8,15,22,29)
const THIRTY = IntSet(1,2,8,9,15,16,22,23,29,30)
const THIRTYONE = IntSet(1,2,3,8,9,10,15,16,17,22,23,24,29,30,31)
function daysofweekinmonth(dt::TimeType)
    y,m,d = yearmonthday(dt)
    ld = daysinmonth(y,m)
    return ld == 28 ? 4 : ld == 29 ? ((d in TWENTYNINE) ? 5 : 4) :
           ld == 30 ? ((d in THIRTY) ? 5 : 4) :
           (d in THIRTYONE) ? 5 : 4
end

@vectorize_1arg TimeType dayname
@vectorize_1arg TimeType dayabbr
@vectorize_1arg TimeType dayofweek
@vectorize_1arg TimeType dayofweekofmonth
@vectorize_1arg TimeType daysofweekinmonth

### Months
const January,February,March,April,May,June = 1,2,3,4,5,6
const July,August,September,October,November,December = 7,8,9,10,11,12
const Jan,Feb,Mar,Apr,May,Jun,Jul,Aug,Sep,Oct,Nov,Dec = 1,2,3,4,5,6,7,8,9,10,11,12
const english_months = [1=>"January",2=>"February",3=>"March",4=>"April",
                        5=>"May",6=>"June",7=>"July",8=>"August",9=>"September",
                        10=>"October",11=>"November",12=>"December"]
const VALUETOMONTH = (UTF8String=>Dict{Int,UTF8String})["english"=>english_months]
const englishabbr_months = [1=>"Jan",2=>"Feb",3=>"Mar",4=>"Apr",
                            5=>"May",6=>"Jun",7=>"Jul",8=>"Aug",9=>"Sep",
                            10=>"Oct",11=>"Nov",12=>"Dec"]
const VALUETOMONTHABBR = (UTF8String=>Dict{Int,UTF8String})["english"=>englishabbr_months]
monthname(dt::Integer;locale::String="english") = VALUETOMONTH[locale][dt]
monthabbr(dt::Integer;locale::String="english") = VALUETOMONTHABBR[locale][dt]
monthname(dt::TimeType;locale::String="english") = VALUETOMONTH[locale][month(dt)]
monthabbr(dt::TimeType;locale::String="english") = VALUETOMONTHABBR[locale][month(dt)]

daysinmonth(dt::TimeType) = ((y,m) = yearmonth(dt); return daysinmonth(y,m))

@vectorize_1arg TimeType monthname
@vectorize_1arg TimeType monthabbr
@vectorize_1arg TimeType daysinmonth

### Years
isleapyear(dt::TimeType) = isleapyear(year(dt))

dayofyear(dt::TimeType) = ((y,m,d) = yearmonthday(dt); return dayofyear(y,m,d))

daysinyear(dt::TimeType) = 365 + isleapyear(dt)

@vectorize_1arg TimeType isleapyear
@vectorize_1arg TimeType dayofyear
@vectorize_1arg TimeType daysinyear

### Quarters
function quarterofyear(dt::TimeType)
    m = month(dt)
    return m < 4 ? 1 : m < 7 ? 2 : m < 10 ? 3 : 4
end
const QUARTERDAYS = [0,90,181,273]
dayofquarter(dt::TimeType) = dayofyear(dt) - QUARTERDAYS[quarterofyear(dt)]

@vectorize_1arg TimeType quarterofyear
@vectorize_1arg TimeType dayofquarter


text:
using BCVUMPS
using Test

@testset "BCVUMPS.jl" begin
    @testset "fixedpoint" begin
        println("fixedpoint tests running...")
        include("fixedpoint.jl")
    end

    @testset "environment" begin
        println("environment tests running...")
        include("environment.jl")
    end

    @testset "bcvumpsruntime" begin
        println("bcvumpsruntime tests running...")
        include("bcvumpsruntime.jl")
    end

    @testset "example tensors" begin
        println("exampletensors tests running...")
        include("exampletensors.jl")
    end

    @testset "MCMC" begin
        println("MCMC tests running...")
        include("MCMC.jl")
    end

    @testset "example obs" begin
        println("example obs tests running...")
        include("exampleobs.jl")
    end
end


text:
"""
    NoRecorder()

Records nothing.
"""
struct NoRecorder end
export NoRecorder
recordspike!(recorder, net, spike) = Nothing # Void

mutable struct SpikeRecorder
    neuronindices::Array{Int64, 1}
    timeintervals::Array{Tuple{Float64, Float64}, 1}
    intervalcounter::Int64
    spikes::Array{Tuple{Int64, Float64}, 1}
end
export SpikeRecorder
"""
    SpikeRecorder(neuronindices, timeintervals)

Records all spikes in `neuronindices` and `timeintervals`.
`neuronindices` is an array of indices.
`timeintervals` should be given as `[(start1, end1); (start2, end2); ...]`.
"""
SpikeRecorder(neuronindices, timeintervals) =
    SpikeRecorder(neuronindices, timeintervals, 1, Tuple{Int64, Float64}[], Tuple{Int64, Float64, Float64}[])
function recordspike!(recorder::SpikeRecorder, net, spike)
    !(spike[1] in recorder.neuronindices) && return
    recorder.intervalcounter > length(recorder.timeintervals) && return
    while spike[2] > recorder.timeintervals[recorder.intervalcounter][2]
        recorder.intervalcounter += 1
        recorder.intervalcounter > length(recorder.timeintervals) && return
    end
    if spike[2] > recorder.timeintervals[recorder.intervalcounter][1]
        push!(recorder.spikes, spike)
    end
    return
end
struct AllSpikesRecorder
    spikes::Array{Tuple{Int64, Float64}, 1}
end
export AllSpikesRecorder
"""
    AllSpikesRecorder() = AllSpikesRecorder(Tuple{Int64, Float64}[])

Records all spikes.
"""
AllSpikesRecorder() = AllSpikesRecorder(Tuple{Int64, Float64}[])
recordspike!(recorder::AllSpikesRecorder, net, spike) = push!(recorder.spikes, spike)

struct AllTracesRecorder
    spikes::Array{Tuple{Int64, Float64, Float64}, 1}
end
export AllTracesRecorder
"""
    AllTracesRecorder() = AllTracesRecorder(Tuple{Int64, Float64, Float64}[])

Records all traces at respective spike times.
"""
AllTracesRecorder() = AllTracesRecorder(Tuple{Int64, Float64, Float64}[])
recordspike!(recorder::AllTracesRecorder, net, spike) = push!(recorder.spikes, (spike[1], spike[2], net.spiketrace[spike[1]]))

struct AllStateRecorder
    spiketimes::Array{Float64, 1}
    spikes::Array{Tuple{Int64, Float64}, 1}
    traces::Array{Array{Float64, 1}, 1}
    membranepotentials::Array{Array{Float64, 1}, 1}
    weights::Array{Array{Float64, 1}, 1}
end
export AllStateRecorder
"""
AllStateRecorder() = AllStateRecorder(Array{Float64, 1}(), Tuple{Int64, Float64}[], Array{Array{Float64, 1}, 1}(),
                                      Array{Array{Float64, 1}, 1}(), Array{Array{Float64, 1}, 1}())

Records the state (spiketime, spike (neuronindex, time), traces, membrane potential, weights) at respective spike times.
"""
AllStateRecorder() = AllStateRecorder(Array{Float64, 1}(), Tuple{Int64, Float64}[], Array{Array{Float64, 1}, 1}(),
                                      Array{Array{Float64, 1}, 1}(), Array{Array{Float64, 1}, 1}())
function recordspike!(recorder::AllStateRecorder, net, spike)
  push!(recorder.spiketimes, deepcopy(net.t))
  push!(recorder.spikes, (spike[1], spike[2]))
  push!(recorder.traces, deepcopy(net.spiketrace))
  push!(recorder.membranepotentials, deepcopy(net.v))
  push!(recorder.weights, deepcopy(net.weights))
end


text:
struct EnsembleBatchQs <: BatchPQs
    model::SVDD.OCClassifier
    k::Int
    sequential_strategy::SequentialPQs
    solver::JuMP.MOI.OptimizerWithAttributes
    model_indices::Vector{Vector{Int}}
    batch_models::Vector{SVDD.OCClassifier}

    function EnsembleBatchQs(model::SVDD.OCClassifier, sequential_strategy::SequentialPQs; k::Int, solver::JuMP.MOI.OptimizerWithAttributes)::EnsembleBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))

        labels = MLLabelUtils.labelmap2vec(model.pools)
        num_observations = size(model.data, 2)
        samples_per_batch = ceil(Int, num_observations/k)

        model_indices = Vector{Vector{Int}}()
        batch_models = Vector{SVDD.OCClassifier}()
        for i in 1:k
            new_model = deepcopy(model)
            push!(batch_models, new_model)
            # sample 1 / k indices
            indices = sample(1:num_observations, samples_per_batch, replace=false)
            push!(model_indices, indices)
        end
        return new(model, k, sequential_strategy, solver, model_indices, batch_models)
    end
end

"""
Select best batch by selecting best observation from each ensemble member.
"""
function select_batch(qs::EnsembleBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    label_list = MLLabelUtils.labelmap2vec(labels)

    batch_samples = []
    for i in 1:qs.k
        # set ensemble model parameters (data + labels)
        indices = qs.model_indices[i]
        model = qs.batch_models[i]
        SVDD.set_data!(model, qs.model.data[:, indices])
        SVDD.set_pools!(model, MLLabelUtils.labelmap(label_list[indices]))

        STDOUT_orig, STDERR_orig = stdout, stderr
        redirect_stdout(); redirect_stderr()
        status = SVDD.fit!(model, qs.solver)
        redirect_stdout(STDOUT_orig); redirect_stderr(STDERR_orig)

        if status !== JuMP.MathOptInterface.OPTIMAL
            throw(ErrorException("Ensemble model could not find optimal solution."))
        end

        candidate_scores = qs_score(qs.sequential_strategy, x, labels)[candidate_indices]
        descending_indices = sortperm(candidate_scores; rev=true)
        push!(batch_samples, first(x for x in candidate_indices[descending_indices] if x ∉ batch_samples))
    end
    return batch_samples
end


text:
# abstract type AbstractBoxPartition{B <: Box} end

"""
A generalized box with
`center`:   vector where the box's center is located
`radius`:   vector of radii, length of the box in each dimension

"""
struct Box{N,T <: AbstractFloat}
    center::SVector{N,T}
    radius::SVector{N,T}

    function Box(center, radius)
        N = length(center)

        if length(radius) != N
            throw(DimensionMismatch("Center vector and radius vector must have same length ($N)"))
        end

        if any(x -> x <= 0, radius)
            throw(DomainError(radius, "radius must be positive in every component"))
        end

        T = promote_type(eltype(center), eltype(radius))
        if !(T <: AbstractFloat)
            T = Float64
        end

        return new{N,T}(SVector{N,T}(center), SVector{N,T}(radius))
    end
end

Base.in(point, box::Box) = all(box.center - box.radius  .<=  point  .<  box.center + box.radius)

volume(box::Box) = prod(2 * box.radius)

function Base.show(io::IO, box::Box)
    print(io, "Box: center = $(box.center), radii = $(box.radius)")
end


text:
using SafeTestsets

# --- TESTS ---
@safetestset "SPICE Utils" begin include("spice_utils.jl") end
@safetestset "NBP" begin include("nbp.jl") end
@safetestset "3BP" begin include("3bp.jl") end
@safetestset "Document tests" begin include("doctests.jl") end
@safetestset "Sensitivity" begin include("sensitivity.jl") end


text:
using Distributions
using FileIO
using POMDPGym
using POMDPs
using Random

include("utils.jl")

# Define SimpleMDP.
struct SimpleMDP <: MDP{Float64, Float64} end

POMDPs.initialstate(mdp::SimpleMDP) = MersenneTwister()
POMDPs.actions(mdp::SimpleMDP) = DiscreteNonParametric([1.0, 2.0], [0.5, 0.5])
POMDPs.isterminal(mdp::SimpleMDP, s) = s > 5
POMDPs.discount(mdp::SimpleMDP) = 1.0

function POMDPs.gen(mdp::SimpleMDP, s, a, rng::AbstractRNG=Random.GLOBAL_RNG; kwargs...)
    sp = s + a
    r = sp > 2 ? 1 : 2
    return (sp=sp, r=r)
end

mdp = SimpleMDP()

fixed_s = rand(initialstate(mdp))

N = 1000
base_N = 1000
c = 0.3
vloss = 0.0
α = 0.1
β = 1.0
γ = 0.0
is_baseline = false

path = "data"

tree_mdp = create_tree_amdp(mdp, actions; reduction="sum")

alpha_list = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

function baseline()
    baseline_out = run_baseline(mdp, fixed_s, actions; N=base_N)
save(joinpath(path, "simple_baseline_$(N).jld2"),
     Dict("risks:" => baseline_out[1], "states:" => baseline_out[2]))
    println("Baseline metrics")

    for alpha in alpha_list
        m = eval_metrics(baseline_out[1]; alpha)
        println("[Alpha=$(alpha)] Mean: $(m.mean), VaR: $(m.var), CVaR: $(m.cvar), Worst: $(m.worst)")
    end
end

function mcts()
    mcts_out, planner = run_mcts(tree_mdp, fixed_s; N=N, c=c, vloss=vloss, α=α, β=β, γ=γ)
    mcts_info = mcts_out[4]
    save(joinpath(path, "simple_mcts_$(N).jld2"),
         Dict("risks:" => mcts_out[1], "states:" => mcts_out[2], "IS_weights:" => mcts_out[3], "tree:" => mcts_info[:tree]))

    println("TIS metrics: N=$(N), c=$(c), α=$(α), β=$(β), γ=$(γ)")
    for alpha in alpha_list
        m = eval_metrics(mcts_out[1]; weights=exp.(mcts_out[3]), alpha=alpha)
        println("[Alpha=$(alpha)] Mean: $(m.mean), VaR: $(m.var), CVaR: $(m.cvar), Worst: $(m.worst)")
    end
end

is_baseline ? baseline() : mcts()


text:
using ProjectX
using Test

@testset "ProjectX.jl" begin
    # Write your own tests here.
end