text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

"""
    module RRRMC

This module implements methods for reduced-rejection-rate Monte Carlo on Ising spin models.

See [`standardMC`](@ref), [`rrrMC`](@ref), [`bklMC`](@ref), [`wtmMC`](@ref) and [`extremal_opt`](@ref).
"""
module RRRMC

export standardMC, rrrMC, bklMC, wtmMC, extremal_opt

using LinearAlgebra, Random
using ExtractMacro

include("DFloats.jl")
include("Common.jl")
include("Interface.jl")
include("DeltaE.jl")
include("WaitingTimes.jl")

using .Common
using .Interface
using .DeltaE
using .WaitingTimes

include("load_graphs.jl")

include("QAliases.jl")
include("REAliases.jl")
include("LEAliases.jl")
include("TLEAliases.jl")
using .QAliases
using .REAliases
using .LEAliases
using .TLEAliases

@inline accept(x) = x ≥ 0 || rand() < exp(x)
@inline function accept(c, x)
    c ≥ 1 && x ≥ 0 && return true
    a = c * exp(x)
    return a ≥ 1 || rand() < a
end

"""
    standardMC(X::AbstractGraph, β::Real, iters::Integer; keywords...)

Runs `iters` iterations of a standard Metropolis Monte Carlo algorithm on the given Ising spin model `X`, at inverse temperature `β`.
Each spin flip attempt counts as an iteration.

Returns two objects: a vector of energies and the last configuration (see [`Config`](@ref)).

Possible keyord arguments are:

* `step`: the interval of iterations with which to collect energies (for the returned result) and to call the `hook` function (see below).
  The default is `1`, which is good for debugging but otherwise generally a bad idea.
* `seed`: the random seed. The default is some arbitrary number.
* `C0`: the initial configuration. The default is `nothing`, in which case it is initialized at random. Otherwise it can be a [`Config`](@ref) object.
  Passing the result of a previous run can be useful e.g. when implementing a simulated annealing protocol, or if the system has not equilibrated yet.
* `hook`: a function to be executed after every `step` number of iterations (see above). It must take five arguments: the current iteration, the graph `X`,
  the current configuration, the number of accepted moves so far, and the current energy. Useful to collect data other than the energy, write to files ecc;
  you'd probably want to use a closure, see the example below. The return value must be a `Bool`: return `false` to interrupt the simulation, `true` otherwise.
  The default does nothing and just returns `true`.

Basic example:

```
julia> Random.seed!(76543); X = RRRMC.GraphPSpin3(3999, 5); β = 1.0;
julia> Es, C = standardMC(X, β, 100_000, step = 1_000);
```

Example of using the `hook` for collecting samples as the columns of a `BitMatrix`:

```
julia> iters = 100_000; step = 1_000; l = iters ÷ step; N = RRRMC.getN(X);
julia> Cs = BitArray(undef, N, l); hook = (it, X, C, acc, E) -> (Cs[:,it÷step]=C.s; true);
julia> Es, C = standardMC(X, β, iters, step = step, hook = hook);
```
"""
function standardMC(X::AbstractGraph{ET}, β::Real, iters::Integer;
                    seed = 167432777111,
                    step::Integer = 1,
                    hook = (x...)->true,
                    C0::Union{Config,Nothing} = nothing,
                    pp = nothing,
                    quiet::Bool = false
                   ) where {ET}
    seed > 0 && Random.seed!(seed)
    Es = empty!(Array{ET}(undef, min(10^8, iters ÷ step)))

    N = getN(X)
    C::Config = C0 ≡ nothing ? Config(N) : C0
    C.N == N || throw(ArgumentError("Invalid C0, wrong N, expected $N, given: $(C.N)"))
    E = energy(X, C)
    accepted = 0
    # pp ≡ nothing && (pp = randperm(N))
    # j = 0
    it = 0
    while it < iters
        it += 1
        # println("it=$it")
        # @assert isapprox(E , energy(X, C), atol=1e-8) (E, energy(X, C))
        if (it % step == 0)
            # println("it=$it")
            # @assert isapprox(E , energy(X, C), atol=1e-8) (E, energy(X, C))
            push!(Es, E)
            hook(it, X, C, accepted, E) || break
        end
        # j += 1
        # j > N && (j -= N)
        # i = pp[j]
        i = rand(1:N)
        ΔE = delta_energy(X, C, i)
        accept(-β * ΔE) || continue
        spinflip!(X, C, i)
        E += ΔE
        accepted += 1
    end
    if !quiet
        println("samples = ", length(Es))
        println("iters = ", it)
        println("accept rate = ", accepted / it)
    end
    # @show E, energy(X, C), abs(E - energy(X, C))
    return Es, C
end

### Begin RRR-related functions

function step_rrr(X, C::Config, ΔEcache)
    z = get_z(ΔEcache)
    move, ΔE = rand_move(ΔEcache)
    compute_staged!(X, C, move, ΔEcache)
    z′ = compute_reverse_probabilities!(ΔEcache)
    c = z / z′
    return move, c, ΔE
end

"""
    rrrMC(X::AbstractGraph, β::Real, iters::Integer; keywords...)

Same as [`standardMC`](@ref), but uses the reduced-rejection-rate method. Each iteration takes more time, but has a higher chance of being accepted,
so fewer iterations overall should be needed normally. Whether this trade-off is convenient depends on the parameters and the details of the model.
This function has specialized versions for [`DiscrGraph`](@ref) and [`DoubleGraph`](@ref) models.

The return values and the keyword arguments are the same as [`standardMC`](@ref), see the usage examples for that function.
"""
function rrrMC(X::SingleGraph{ET}, β::Real, iters::Integer;
               seed = 167432777111,
               step::Integer = 1,
               hook = (x...)->true,
               C0::Union{Config,Nothing} = nothing,
               staged_thr::Real = NaN,
               staged_thr_fact::Real = 5.0,
               quiet::Bool = false
              ) where {ET}

    isfinite(β) || throw(ArgumentError("β must be finite, given: $β"))
    seed > 0 && Random.seed!(seed)
    Es = empty!(Array{ET}(undef, min(10^8, iters ÷ step)))

    if staged_thr ≡ NaN
        staged_thr = X isa SimpleGraph ? 0.8 : 0.5
    end

    N = getN(X)
    C::Config = C0 ≡ nothing ? Config(N) : C0
    C.N == N || throw(ArgumentError("Invalid C0, wrong N, expected $N, given: $(C.N)"))
    E = energy(X, C)
    ΔEcache = gen_ΔEcache(X, C, β)
    #check_consistency(ΔEcache)

    λ = staged_thr_fact / N
    staged_its = 0

    it = 0
    accepted = 0
    acc_rate = 0.5
    while it < iters
        #@assert E == energy(X, C)
        #check_consistency(ΔEcache)
        it += 1
        if (it % step == 0)
            push!(Es, E)
            hook(it, X, C, accepted, E) || break
        end
        acc = false
        if acc_rate < staged_thr
            staged_its += 1
            move, c, ΔE = step_rrr(X, C, ΔEcache)
            if rand() < c
                spinflip!(X, C, move)
                apply_staged!(ΔEcache)
                E += ΔE
                accepted += 1
                acc = true
            end
        else
            move, ΔE = rand_move(ΔEcache)
            c = apply_move!(X, C, move, ΔEcache)
            if rand() < c
                E += ΔE
                accepted += 1
                acc = true
            else
                apply_move!(X, C, move, ΔEcache)
            end
        end
        acc_rate = acc_rate * (1 - λ) + acc * λ
    end
    if !quiet
        println("samples = ", length(Es))
        println("iters = ", it)
        println("accept rate = ", accepted / it)
        println("frac. staged iters = ", staged_its / it)
    end
    return Es, C
end

function rrrMC(X::DoubleGraph{GT,ET}, β::Real, iters::Integer;
               seed = 167432777111,
               step::Integer = 1,
               hook = (x...)->true,
               C0::Union{Config,Nothing} = nothing,
               staged_thr::Real = 0.5,
               staged_thr_fact::Real = 5.0,
               quiet::Bool = false
              ) where {GT,ET}
    isfinite(β) || throw(ArgumentError("β must be finite, given: $β"))
    seed > 0 && Random.seed!(seed)
    Es = empty!(Array{ET}(undef, min(10^8, iters ÷ step)))

    N = getN(X)
    C::Config = C0 ≡ nothing ? Config(N) : C0
    C.N == N || throw(ArgumentError("Invalid C0, wrong N, expected $N, given: $(C.N)"))
    E = energy(X, C)
    @assert isfinite(E)
    X0 = inner_graph(X)
    ΔEcache = gen_ΔEcache(X0, C, β)
    #check_consistency(ΔEcache)

    λ = staged_thr_fact / N
    staged_its = 0

    it = 0
    accepted = 0
    acc_rate = 0.5
    while it < iters
        #@assert abs(E - energy(X, C)) < 1e-10 (E, energy(X, C), abs(E - energy(X,C)))
        #DeltaE.check_consistency(ΔEcache)
        it += 1
        if (it % step == 0)
            push!(Es, E)
            hook(it, X, C, accepted, E) || break
        end
        acc = false
        if acc_rate < staged_thr
            staged_its += 1
            move, c, ΔE0 = step_rrr(X0, C, ΔEcache)
            ΔE1 = delta_energy_residual(X, C, move)
            if accept(c, -β * ΔE1)
                spinflip!(X, C, move)
                apply_staged!(ΔEcache)
                E += ΔE0 + ΔE1
                accepted += 1
                acc = true
            end
        else
            move, ΔE0 = rand_move(ΔEcache)
            ΔE1 = delta_energy_residual(X, C, move)
            c = apply_move!(X, C, move, ΔEcache)
            if accept(c, -β * ΔE1)
                E += ΔE0 + ΔE1
                accepted += 1
                acc = true
            else
                apply_move!(X, C, move, ΔEcache)
            end
        end
        acc_rate = acc_rate * (1 - λ) + acc * λ
    end
    if !quiet
        println("samples = ", length(Es))
        println("iters = ", it)
        println("accept rate = ", accepted / it)
        println("frac. staged iters = ", staged_its / it)
    end
    return Es, C
end

### Begin BKL-related functions

apply_step_bkl!(X::AbstractGraph, C::Config, move::Int, ΔEcache) =
    apply_move!(X, C, move, ΔEcache, Val{false})

apply_step_bkl!(X::DiscrGraph, C::Config, move::Int, ΔEcache) =
    apply_move!(X, C, move, ΔEcache)

"""
    bklMC(X::AbstractGraph, β::Real, iters::Integer; keywords...)

Same as [`standardMC`](@ref), but uses the rejection-free method by Bortz, Kalos and Lebowitz. Each step takes more time, but rejected moves
are almost free, since they are skipped entirely. This function has a specialized version for [`DiscrGraph`](@ref) models.

The return values and the keyword arguments are the same as [`standardMC`](@ref), see the usage examples for that function.

Note that the number of iterations includes the rejected moves. This makes the results directly comparable with those of `standardMC`. It also
means that increasing `β` at fixed `iters` will result in fewer steps being actually computed.
"""
function bklMC(X::AbstractGraph{ET}, β::Real, iters::Integer;
               seed = 167432777111,
               step::Integer = 1,
               hook = (x...)->true,
               C0::Union{Config,Nothing} = nothing,
               quiet::Bool = false
              ) where {ET}
    seed > 0 && Random.seed!(seed)
    Es = empty!(Array{ET}(undef, min(10^8, iters ÷ step)))

    N = getN(X)
    C::Config = C0 ≡ nothing ? Config(N) : C0
    C.N == N || throw(ArgumentError("Invalid C0, wrong N, expected $N, given: $(C.N)"))
    E = energy(X, C)
    ΔEcache = gen_ΔEcache(X, C, β, false)
    #check_consistency(ΔEcache)

    it = 0
    accepted = 0
    nextstep = step
    stop = false
    while it < iters
        #@assert E == energy(X, C)
        #check_consistency(ΔEcache)

        skip = rand_skip(ΔEcache)
        move, ΔE = rand_move(ΔEcache)

        while it + skip + 1 ≥ nextstep
            push!(Es, E)
            hook(nextstep, X, C, accepted, E) || @goto out
            nextstep += step
            nextstep > iters && @goto out
        end

        apply_step_bkl!(X, C, move, ΔEcache)
        it += skip + 1
        E += ΔE
        accepted += 1
    end
    @label out
    if !quiet
        println("samples = ", length(Es))
        println("iters = ", it)
        println("accept rate = ", accepted / iters)
        println("true it = ", accepted)
    end
    return Es, C
end

"""
    wtmMC(X::AbstractGraph, β::Real, samples::Integer; keywords...)

Same as [`standardMC`](@ref), but uses the rejection-free waiting-time method by Dall and Sibani. It is similar to [`bklMC`](@ref).

The return values and the keyword arguments are *almost* the same as [`standardMC`](@ref), see the usage examples for that function.
However, the waiting time method uses an internal "global time" variable, which takes the place of the "iterations" counter of [`standardMC`](@ref)
and of [`bklMC`](@ref). Thus, this function has two differences with respect to the other samplers in the module:

* the function takes a `samples` integer argument, with the maximum number of samples which will be collected.
* the `step` keyword argument is of type `Float64` instead of integer (default value = `1.0`; you'll probably want to change this).
  The `step` is measured in terms of the global time variable and is scaled with the size of the problem `N`.

The total number of samples actually collected can still be less than `samples` if the `hook` function from the keyword arguments returns `false` earlier.
"""
function wtmMC(X::AbstractGraph{ET}, β::Real, samples::Integer;
               seed = 167432777111,
               step::Float64 = 1.0,
               hook = (x...)->true,
               C0::Union{Config,Nothing} = nothing,
               quiet::Bool = false
              ) where {ET}
    seed > 0 && Random.seed!(seed)
    Es = empty!(Array{ET}(undef, min(10^8, samples)))

    N = getN(X)
    C::Config = C0 ≡ nothing ? Config(N) : C0
    C.N == N || throw(ArgumentError("Invalid C0, wrong N, expected $N, given: $(C.N)"))
    E = energy(X, C)
    theap = THeap(X, C, β)

    step /= N

    num_moves = 0
    tmax = step * samples
    t = 0.0
    nextstep = step
    while t < tmax
        #@assert abs(E - energy(X, C)) < 1e-10 E-energy(X,C)

        t′, move = pick_next(theap)
        while t′ ≥ nextstep
            push!(Es, E)
            hook(nextstep, X, C, num_moves, E) || @goto out
            nextstep += step
            nextstep > tmax + 1e-10 && @goto out
        end
        t = t′
        ΔE = update_heap!(theap, X, C, β, move, t)
        E += ΔE

        num_moves += 1
    end
    @label out
    if !quiet
        println("samples = ", length(Es))
        println("num_moves = ", num_moves)
        println("global time = ", t)
        println("ratio = ", t / num_moves)
    end
    return Es, C
end

# τ-Extremal Optimization

"""
    extremal_opt(X::AbstractGraph, τ::Real, iters::Integer; keywords...)

Extremal optimization algorithm: it seeks the lowest energy state of a given Ising spin model `X` by performing
a random walk biased towards low-energy states, but with a heavy tail which allows it to avoid getting trapped.

The interface is very similar to that of [`standardMC`](@ref) and the other RRRMC Monte Carlo functions, but it
takes a parameter `τ` instead of `β`, the return values are different, and the `hook` keyword argument has a
different signature, see below.

The parameter `τ` controls the shape of the tail and should be larger than 1 (a reasonable value could be 1.3);
`iters` is the total number of spin flips performed.

This function has a specialized version for [`DiscrGraph`](@ref) models; otherwise, it works but it is not
implemented very efficiently at the moment (each spin flip takes O(N) time even on diluted graphs).

Returns 4 objects: the final configuration, the minimum energy found, the configuration of minimum energy, and the
iteration at which such configuration was found (see [`Config`](@ref)).

Possible keyord arguments are:
* `step`: the interval of iterations with which to call the `hook` function (see below).
  The default is `1`, which is good for debugging but otherwise generally a bad idea, probably.
* `seed`: the random seed. The default is some arbitrary number.
* `C0`: the initial configuration. The default is `nothing`, in which case it is initialized at random. Otherwise it can be a [`Config`](@ref) object.
* `hook`: a function to be executed after every `step` number of iterations (see above). It must take five arguments: the current iteration, the graph `X`,
  the current configuration, the current energy, and the minimum energy found so far. Useful to collect data other than the energy, write to files ecc;
  you'd probably want to use a closure, see the example below. The return value must be a `Bool`: return `false` to interrupt the simulation, `true` otherwise.
  The default does nothing and just returns `true`. Note that the signature is similar to that of the corresponding `hook` argument in [`standardMC`](@ref),
  but different.

Basic example:

```
julia> Random.seed!(76543); X = RRRMC.GraphPSpin3(3999, 5); τ = 1.3;
julia> C, Emin, Cmin, itmin = extremal_opt(X, τ, 100_000, step = 1_000);
```

Example of using the `hook` for collecting samples as the columns of a `BitMatrix`:

```
julia> iters = 100_000; step = 1_000; l = iters ÷ step; N = RRRMC.getN(X);
julia> Cs = BitArray(undef, N, l); hook = (it, X, C, E, Emin) -> (Cs[:,it÷step]=C.s; true);
julia> C, Emin, Cmin, itmin = extremal_opt(X, τ, iters, step = step, hook = hook);
```

"""
function extremal_opt(X::AbstractGraph{ET}, τ::Real, iters::Integer;
                      seed = 167432777111,
                      step::Integer = 1,
                      hook = (x...)->true,
                      C0::Union{Config,Nothing} = nothing,
                      quiet::Bool = false
                     ) where {ET}
    seed > 0 && Random.seed!(seed)

    N = getN(X)
    C::Config = C0 ≡ nothing ? Config(N) : C0
    C.N == N || throw(ArgumentError("Invalid C0, wrong N, expected $N, given: $(C.N)"))
    E = energy(X, C)
    Emin = E
    Cmin = copy(C)
    itmin = 0
    ΔEcache = gen_EOcache(X, C, τ)
    # check_consistency(ΔEcache)

    it = 0
    stop = false
    while it < iters
        it += 1
        # @assert E == energy(X, C)
        # check_consistency(ΔEcache)

        if (it % step == 0)
            hook(it, X, C, E, Emin) || break
            # E = energy(X, C)
        end

        move, ΔE = rand_move(ΔEcache)
        apply_move!(X, C, move, ΔEcache)
        E += ΔE

        if E < Emin
            Emin = E
            copy!(Cmin, C)
            itmin = it
        end
    end
    if !quiet
        # println("samples = ", length(Es))
        println("iters = ", it)
        println("min [it = $itmin] = $Emin")
    end
    return C, Emin, Cmin, itmin
end

### auxiliary/miscellanaous/ugly stuff. Hic sunt leones.

ba2int(b::BitArray) = Int(b.chunks[1])
int2ba(i::Int, N) = (@assert N≤64; b = BitVector(undef, N); b.chunks[1]=i; b)

function truep(X::DiscrGraph, β::Real)
    N = X.N
    @assert N ≤ 64
    S = 1 << N
    C = Config(N)
    p = zeros(S)
    Z = 0.0
    for i = 1:S
        C.s.chunks[1] = i - 1
        E = energy(X, C)
        p[i] = exp(-β * E)
        Z += p[i]
    end
    p ./= Z
    return p
end

function samplep(ws::Vector{BitVector})
    N = length(ws[1])
    S = 1 << N
    C = MCNew.Config(N)
    p = zeros(S)
    for w in ws
        copy!(C.s, w)
        i = ba2int(w) + 1
        p[i] += 1
    end
    p ./= length(ws)
    return p
end

function tm(Es::Vector{T}, step = 1, skip0 = 0.1, skip1 = 0.05) where {T<:Real}
    N = length(Es)
    i0 = floor(Int, N * skip0)
    N0 = N - i0
    n = N0 ÷ step
    m = zeros(n)
    m[1] = mean(Es[i0 + (1:step)])
    for j = 2:n
        m[j] = (mean(Es[i0 + (j-1)*step + (1:step)]) + m[j-1] * (j-1)) / j;
    end
    return m[1 + floor(Int, skip1 * n):end]
end

function ravg(Es::Vector{T}, step = 1_000, skip0 = 0.0) where {T<:Real}
    N = length(Es)
    i0 = floor(Int, N * skip0)
    N0 = N - i0
    n = N0 ÷ step
    m = zeros(n)
    m[1] = mean(Es[i0 + (1:step)])
    for j = 2:n
        m[j] = mean(Es[i0 + (j-1)*step + (1:step)])
    end
    return m
end

function second_eigenvalue(Q::AbstractMatrix{Float64})
    ev = eigvals(Q)
    @assert all(x->(abs(imag(x)) ≤ 1e-10), ev)
    evr = Float64[real(x) for x in ev]

    τ = -1 / log(sort!(evr)[end-1])
end

function second_eigenvalue_standard(X::DiscrGraph, β::Float64)
    N = getN(X)
    @assert N ≤ 64
    S = 1 << N
    Q = zeros(S, S)
    C = Config(N)
    p = zeros(S)
    for ks = 1:S
        C.s.chunks[1] = ks - 1
        pT = 0.0
        for j = 1:N
            kd = Int(((ks-1) ⊻ (UInt64(1) << (j-1))) + 1)

            #=
            C.s[j] $= 1
            @assert C.s.chunks[1] == kd
            C.s[j] $= 1
            kd += 1
            =#

            ΔE = delta_energy(X, C, j)
            p = min(1.0, exp(-β * ΔE)) / N
            pT += p
            Q[kd, ks] = p
        end
        Q[ks, ks] = 1 - pT
    end

    τ = second_eigenvalue(Q)

    return Q, τ
end

function second_eigenvalue_bkl(Q::Matrix{Float64})
    pr = diag(Q)
    rfQ = (Q - diagm(pr)) ./ (1 .- pr')
    τ = second_eigenvalue(rfQ)

    return rfQ, τ
end

function second_eigenvalue_bkl(X::DiscrGraph, β::Float64)
    Q, _ = second_eigenvalue_standard(X, β)
    return second_eigenvalue_bkl(Q)
end

function second_eigenvalue_rrr(X::DiscrGraph, β::Float64)
    N = getN(X)
    @assert N ≤ 64
    S = 1 << N

    Q = zeros(S, S)
    C = Config(N)

    for ks = 1:S
        C.s.chunks[1] = ks - 1

        energy(X, C)
        ΔEcache = gen_ΔEcache(X, C, β)
        #DeltaE.check_consistency(ΔEcache)

        pchg = 0.0
        for j = 1:N
            @extract ΔEcache : k=pos[j] z
            compute_staged!(X, C, j, ΔEcache)
            z′ = compute_reverse_probabilities!(ΔEcache)
            w = get_class_f(ΔEcache, k)

            # pp = min(1.0, z / z′) * w / z # this is error-prone when z is small
            pp = w / max(z, z′)
            @assert isfinite(pp)
            pchg += pp

            kd = Int(((ks-1) ⊻ (UInt64(1) << (j-1))) + 1)
            Q[kd, ks] = pp
        end
        @assert -1e-15 ≤ 1 - pchg ≤ 1 + 1e-15
        Q[ks, ks] = clamp(1 - pchg, 0.0, 1.0)
    end

    τ = second_eigenvalue(Q)

    return Q, τ
end

function second_eigenvalue_stats(;seed::Integer = 86823, graph = GraphRRG, args = (10, 5), βr = 0.5:0.5:3.0, n = 10)
    τs = Matrix{Float64}[]
    rrs = Matrix{Float64}[]
    for j = 1:n
        Random.seed!(seed + j)
        println("seed = $(seed + j)")
        X = graph(args...)
        τ = Array{Float64}(undef, length(βr), 3)
        rr = Array{Float64}(undef, length(βr), 3)
        for (l,β) in enumerate(βr)
            println("  β=$β")
            p = truep(X, β);
            Q, τ[l,1] = second_eigenvalue_standard(X, β)
            @assert maximum(abs.(p - Q * p)) < 1e-13
            rr[l,1] = sum(diag(Q) .* p)
            rfQ, τ[l,2] = second_eigenvalue_bkl(Q)
            pr = diag(Q)
            @assert maximum(abs.(p .* (1 .- pr) - rfQ * (p .* (1 .- pr)))) < 1e-13
            rr[l,2] = 0.0
            Q, τ[l,3] = second_eigenvalue_rrr(X, β)
            @assert maximum(abs(p - Q * p)) < 1e-13
            rr[l,3] = sum(diag(Q) .* p)
            println("    τ:  ", τ[l,:])
            println("    τr: ", τ[l,1] ./ τ[l,:])
            println("    rr: ", rr[l,:])
        end
        push!(τs, τ)
        push!(rrs, rr)
    end
    println("---------\n")
    println("τ stats1")
    τm = mean(τs)
    for (l,β) in enumerate(βr)
        println("  β=$β")
        println("    ", τm[l,:])
        println("    ", τm[l,1] ./ τm[l,:])
    end
    println("---------\n")
    println("τ stats2")
    τr = mean(Matrix{Float64}[τ[:,1] ./ τ for τ in τs])
    for (l,β) in enumerate(βr)
        println("  β=$β")
        println("    ", τr[l,:])
    end
    println("---------\n")
    println("rr stats")
    rrm = mean(rrs)
    for (l,β) in enumerate(βr)
        println("  β=$β")
        println("    ", rrm[l,:])
    end
    println("---------\n")
    return τs, rrs
end

function runtest(;N = 1_000, D = 4, β = 1.0, seedlist = 101:110, iters = 10_000_000, stepfact = 10, step = 100)
    Random.seed!(111)
    X = GraphRRG(N, D)

    nsamples = iters ÷ step

    nseeds = length(seedlist)
    EsBKL = zeros(nsamples, nseeds)
    Es5 = zeros(nsamples, nseeds)
    elbkl = 0.0
    elsmrt = 0.0
    for (l,seed) in enumerate(seedlist)
        info("SEED = $seed")
        #info("  bkl")
        #el = @elapsed EsBKL[:,l] = bkl(X, β, iters * stepfact, seed=seed, step=step*stepfact)
        #info("    time: $el")
        #elbkl += el
        info("  rrr")
        el = @elapsed Es5[:,l] = rrrMC(X, β, iters, seed=seed, step=step)
        info("    time: $el")
        elsmrt += el
        println()
    end
    info("average elapsed time bkl:  ", elbkl / nseeds)
    info("average elapsed time smrt: ", elsmrt / nseeds)
    info("ratio: ", elsmrt / elbkl)
    return EsBKL, Es5
end

end


text:
module DummyRecalls

@noinline notsupported() = error("tracing requiers Recalls.jl")

macro note(_...)
    :(notsupported())
end

end  # module

ext: jl


text:
# ---
# title: 630. Course Schedule III
# id: problem630
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy
# link: <https://leetcode.com/problems/course-schedule-iii/description/>
# hidden: true
# ---
#
# There are `n` different online courses numbered from `1` to `n`. Each course
# has some duration(course length) `t` and closed on `dth` day. A course should
# be taken **continuously** for `t` days and must be finished before or on the
# `dth` day. You will start at the `1st` day.
#
# Given `n` online courses represented by pairs `(t,d)`, your task is to find
# the maximal number of courses that can be taken.
#
# **Example:**
#
#
#
#     Input: [[100, 200], [200, 1300], [1000, 1250], [2000, 3200]]
#     Output: 3
#     Explanation:
#     There're totally 4 courses, but you can take 3 courses at most:
#     First, take the 1st course, it costs 100 days so you will finish it on the 100th day, and ready to take the next course on the 101st day.
#     Second, take the 3rd course, it costs 1000 days so you will finish it on the 1100th day, and ready to take the next course on the 1101st day.
#     Third, take the 2nd course, it costs 200 days so you will finish it on the 1300th day.
#     The 4th course cannot be taken now, since you will finish it on the 3300th day, which exceeds the closed date.
#
#
#
#
# **Note:**
#
#   1. The integer 1 <= d, t, n <= 10,000.
#   2. You can't take two courses simultaneously.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "Evolution" begin

    struct MockCand <: AbstractCandidate
        val::Real
    end
    NaiveGAflux.fitness(c::MockCand) = c.val

    struct DummyCand <: AbstractCandidate end

    @testset "EliteSelection" begin
        pop = MockCand.([3, 7, 4, 5, 9, 0])
        @test fitness.(evolve(EliteSelection(3), pop)) == [9, 7, 5]
    end

    @testset "SusSelection" begin
        pop = MockCand.(0.1:0.1:0.8)
        @test fitness.(evolve(SusSelection(4, NoOpEvolution(), MockRng([0.05])), pop)) == [0.1, 0.4, 0.6, 0.7]
    end

    @testset "TournamentSelection" begin
        pop = MockCand.(1:10)
        Random.shuffle(::MockRng, a::AbstractVector) = reverse(a)
        Random.rand(r::MockRng, n::Int) = map(i -> rand(r), 1:n)

        @test fitness.(evolve(TournamentSelection(4, 2, 1, NoOpEvolution(), MockRng([0.05])), pop)) == [10, 8, 6, 4]
        @test fitness.(evolve(TournamentSelection(8, 3, 0.7, NoOpEvolution(), MockRng([0.01, 0.9, 0.01])), pop)) == [9, 9, 7, 4, 2, 8, 5, 2]

        @test fitness.(evolve(TournamentSelection(20, 2, 1, NoOpEvolution(), MockRng([0.01])), pop)) == repeat(10:-2:2, 4)

        # Edge case, nselect*k % popsize == 0
        @test fitness.(evolve(TournamentSelection(5, 2, 1, NoOpEvolution(), MockRng([0.01])), pop)) == 10:-2:2
        # Edge case, nselect*k % popsize == 1
        @test fitness.(evolve(TournamentSelection(2, 2, 1, NoOpEvolution(), MockRng([0.01])), MockCand.(1:3))) == [3, 2]
    end

    @testset "CombinedEvolution" begin
        pop = [DummyCand() for i in 1:5]
        @test evolve(CombinedEvolution(NoOpEvolution(), NoOpEvolution()), pop) == vcat(pop,pop)
    end

    @testset "EvolutionChain" begin
        dfitness = EvolveCandidates(mc -> MockCand(2 * fitness(mc)))
        sqrfitness = EvolveCandidates(mc -> MockCand(fitness(mc)^2))
        @test evolve(EvolutionChain(), MockCand.(1:5)) == MockCand.(1:5)
        @test evolve(EvolutionChain(dfitness, dfitness, sqrfitness), MockCand.(1:10)) == MockCand.((4:4:40).^2)
    end

    @testset "PairCandidates" begin
        struct PairConsumer <: AbstractEvolution end
        nseen = 0
        function NaiveGAflux._evolve(::PairConsumer, pop)
            nseen = length(pop)
            @test typeof(pop) == Vector{Tuple{MockCand, MockCand}}
            return pop
        end

        @test evolve(PairCandidates(NoOpEvolution()), []) == []
        pop = MockCand.(1:9)
        @test evolve(PairCandidates(PairConsumer()), pop) == pop
        @test nseen == 5
        pop = MockCand.(1:4)
        @test evolve(PairCandidates(PairConsumer()), pop) == pop
        @test nseen == 2
    end

    @testset "ShuffleCandidates" begin
        pop = MockCand.(1:10)
        @test evolve(ShuffleCandidates(MersenneTwister(1)), pop) != pop
    end

    @testset "EvolveCandidates" begin
        pop = MockCand.([2,3,5,9])
        @test fitness.(evolve(EvolveCandidates(mc -> MockCand(2mc.val)), pop)) == 2fitness.(pop)
    end

end


text:
#!/usr/bin/env julia
using DelimitedFiles: readdlm

function solve(case, R, C, grid)
    numl = 0

    function is_l(R1, R2, C)
        0 < R1 <= size(grid, 1) && 0 < R2 <= size(grid, 1) && @views all(grid[R1:R2, C]) && (numl += 1; @info "" R1 R2 C)
    end

    for r in 1:R
        len = 0
        runr = 0
        runc = 0
        for c in 1:C
            if grid[r, c]
                len += 1
                len == 1 && (runr = r; runc = c)

                # looking for long
                for d in 2:len
                    is_l(runr, runr + 2d - 1, runc)
                    is_l(runr - 2d + 1, runr, runc)
                    is_l(r, r + 2d - 1, c)
                    is_l(r - 2d + 1, r, c)
                end

                # looking for short
                for d in 4:2:len
                    is_l(r, r + d ÷ 2 - 1, c)
                    is_l(r - d ÷ 2 + 1, r, c)
                    is_l(runr, runr + len ÷ 2 - 1, runc)
                    is_l(runr - len ÷ 2 + 1, runr, runc)
                end
            else
                len = 0
            end
        end
    end

    println("Case #$case: $numl")
end

function main()
    T = parse(Int, readline(stdin))

    for case in 1:T
        R, C = map(x -> parse(Int, x), split(readline(stdin)))
        grid = BitMatrix(undef, R, C)
        for r in 1:R
            grid[r,:] .= map(x -> parse(Bool, x), split(readline(stdin)))
        end
        @info "" grid
        solve(case, R, C, grid)
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
module TSClassifiers

using Random
using CSV
using DataFrames: DataFrame, nrow
using Dates
using Serialization

using ..Pipelines
using ..DecisionTreeLearners: RandomForest
using ..Statifiers
using ..ValDateFilters
using ..AbsTypes
using ..Utils
import ..AbsTypes: fit, fit!, transform, transform!

export fit, fit!, transform, transform!
export TSClassifier, getstats

"""
    TSClassifier(
       Dict(
          # training directory
          :trdirectory => "",
          :tstdirectory => "",
          :modeldirectory => "",
          :feature_range => 7:20,
          :juliarfmodelname => "juliarfmodel.serialized",
          # Output to train against
          # (:class).
          :output => :class,
          # Options specific to this implementation.
          :impl_args => Dict(
             # Merge leaves having >= purity_threshold CombineMLd purity.
             :purity_threshold => 1.0,
             # Maximum depth of the decision tree (default: no maximum).
             :max_depth => -1,
             # Minimum number of samples each leaf needs to have.
             :min_samples_leaf => 1,
             # Minimum number of samples in needed for a split.
             :min_samples_split => 2,
             # Minimum purity needed for a split.
             :min_purity_increase => 0.0
          )
       )
    )

Given a bunch of time-series with specific types. Get the statistical features of each,
use these as inputs to RF classifier with output as the TS type, train and test. Another
option is to use these stat features for clustering and check cluster quality. If
accuracy is poor, add more stat features and repeat same process as outlined for training
and testing. Assume that each time-series is named based on their type which will be
used as target output. For example, temperature time series will be named as temperature?.csv
where ? is an integer. Loop over each file in a directory, get stat and
record in a dictionary/dataframe, train/test. Default to using RandomForest
for classification of data types.
"""
mutable struct TSClassifier <: Learner
   name::String
   model::Dict{Symbol,Any}

   function TSClassifier(args=Dict())
      default_args = Dict{Symbol,Any}(
          :name => "tscl",
          # training directory
          :trdirectory => "",
          :tstdirectory => "",
          :modeldirectory => "",
          :feature_range => 7:20,
          :juliarfmodelname => "juliarfmodel.serialized",
          # Output to train against
          # (:class).
          :output => :class,
          # Options specific to this implementation.
          :impl_args => Dict{Symbol,Any}(
              # Merge leaves having >= purity_threshold CombineMLd purity.
              :purity_threshold => 1.0,
              # Maximum depth of the decision tree (default: no maximum).
              :max_depth => -1,
              # Minimum number of samples each leaf needs to have.
              :min_samples_leaf => 1,
              # Minimum number of samples in needed for a split.
              :min_samples_split => 2,
              # Minimum purity needed for a split.
              :min_purity_increase => 0.0
          )
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

"""
    fit!(tsc::TSClassifier, features::T=[], labels::Vector=[])

Get the stats of each file, collect as dataframe, and train.
"""
function fit!(tsc::TSClassifier, features::DataFrame=DataFrame(), labels::Vector=[])::Nothing
   ispathnotempty(tsc.model) || throw(ArgumentError("empty training/testing/modeling directory"))
   ldirname = tsc.model[:trdirectory]
   mdirname = tsc.model[:modeldirectory]
   modelfname=tsc.model[:juliarfmodelname]
   trdata = getstats(ldirname)
   rfmodel = RandomForest(tsc.model)
   xfeatures = tsc.model[:feature_range]
   X=trdata[:,xfeatures]
   Y=trdata[:,:dtype]
   fit!(rfmodel,X,Y)
   mkpath(mdirname)
   serializedmodel = joinpath(mdirname,modelfname)
   open(serializedmodel,"w") do file
      serialize(file,rfmodel)
   end
   trstatfname = joinpath(mdirname,modelfname*".csv")
   trdata |> CSV.write(trstatfname)
   tsc.model[:features] = names(X)
   tsc.model[:fmodel] = rfmodel
   return nothing
end

function fit(tsc::TSClassifier, features::DataFrame=DataFrame(), labels::Vector=[])::TSClassifier
   fit!(tsc,features,labels)
   return deepcopy(tsc)
end

"""
    transform!(tsc::TSClassifier, features::T=[])

Apply the learned parameters to the new data.
"""
function transform!(tsc::TSClassifier, features::DataFrame=DataFrame())::DataFrame
   ldirname = tsc.model[:tstdirectory]
   mdirname = tsc.model[:modeldirectory]
   modelfname=tsc.model[:juliarfmodelname]
   trdata = getstats(ldirname)
   xfeatures = tsc.model[:feature_range]
   X=trdata[:,xfeatures]
   mfeatures=tsc.model[:features]
   (sum(names(X) .== mfeatures ) == length(mfeatures)) || error("features mismatch")
   serializedmodel = joinpath(mdirname,modelfname)
   if isfile(serializedmodel)
      println("loading model from file: "*serializedmodel)
      model=open(serializedmodel,"r") do file
         deserialize(file)
      end
   else
      model= tsc.model[:fmodel]
   end
   mpred = transform!(model,X)
   return DataFrame(fname=trdata.fname,predtype=mpred)
end

function transform(tsc::TSClassifier, features::DataFrame=DataFrame())::DataFrame
   return transform!(tsc,features)
end

@enum TSType begin
  temperature = 1
  weather = 2
  footfall = 3
  AirOffTemp = 4
  Energy = 5
  Pressure = 6
  RetTemp = 7
end

function ispathnotempty(margs::Dict)
  return margs[:trdirectory] != "" && margs[:tstdirectory] != "" && margs[:modeldirectory] != ""
end

# return stat of a file
function getfilestat(ldirname::AbstractString,lfname::AbstractString)
   myregex = r"(?<dtype>[A-Z _ - a-z]+)(?<number>\d*).(?<ext>\w+)"
   m=match(myregex,lfname)
   ext = m[:ext]; dtype=m[:dtype];num = m[:number]
   (dtype != "" && ext != "")  || error("wrong filename format: dtype[n].csv")
   dtype in string.(instances(TSType)) || error(dtype * ", filename does not indicate known data type.")
   # create a pipeline to get stat
   fname = joinpath(ldirname,lfname)
   csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
   valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
   valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
   stfier = Statifier(Dict(:processmissing=>false))
   mpipeline = @pipeline csvfilter |> valgator |> valnner |> stfier
   df = fit_transform!(mpipeline)
   df.dtype = repeat([dtype],nrow(df))
   df.fname = repeat([lfname],nrow(df))
   return (df)
end

function serialloop(ldirname,mfiles)
   trdata = DataFrame()
   for file in mfiles
      try
         df=getfilestat(ldirname,file)
         trdata = vcat(trdata,df)
         println("getting stats of "*file)
      catch errormsg
         println("skipping "*file*": "*string(errormsg))
         continue
      end
   end
   return trdata
end

function threadloop(ldirname,mfiles)
   @eval using Base.Threads
   trdata = DataFrame()
   mutex = SpinLock()
   Base.Threads.@threads for file in mfiles
      try
         df=getfilestat(ldirname,file)
         lock(mutex)
         trdata = vcat(trdata,df)
         println("getting stats of "*file*" on thread:"*string(Base.Threads.threadid()))
         unlock(mutex)
      catch errormsg
         println("skipping "*file*": "*string(errormsg))
      end
   end
   return trdata
end

# loop over the directory and get stats of each file
# return a dataframe containing stat features and ts type for target
function getstats(ldirname::AbstractString)
   ldirname != "" || error("directory name empty")
   mfiles = readdir(ldirname) |> x->filter(y->match(r".csv",y) != nothing,x)
   mfiles != [] || error("empty csv directory")
   #df = serialloop(ldirname,mfiles)
   # get julia version and run threads if julia 1.3
   jversion = string(Base.VERSION)
   df = DataFrame()
   if match(r"^1.5",jversion) === nothing
      df = serialloop(ldirname,mfiles)
   else
      df = threadloop(ldirname,mfiles)
   end
   return df
end

end


text:
function testAddVariable(client, context, variableLabels, variableTypes, variableTypeStrings)
    resultIds = Task[]
    for (index, label) in enumerate(variableLabels)
        resultId = addVariable(client, context, Variable(label, variableTypes[index]))
        @test resultId != "Error"
        push!(resultIds, resultId)
    end

    # Wait for them to be done before proceeding.
    waitForCompletion(client, resultIds, expectedStatuses=["Complete"], maxSeconds=180)

    return resultIds
end

function testLs(client, context, variableLabels, variableTypes, variableTypeStrings)
    @test length( setdiff(variableLabels, fetch( ls(client, context) )) ) == 0
end

function testGetVariable(client, context, variableLabels, variableTypes, variableTypeStrings)
    for i in 1:length(variableLabels)
        actualVariable = getVariable(client,context,variableLabels[i]) |> fetch
        @test actualVariable["label"] == variableLabels[i]
        @test actualVariable["variableType"] == variableTypeStrings[i]
    end
end

function testGetVariables(client, context, variableLabels, variableTypes, variableTypeStrings)
    # Make a quick dictionary of the expected variable Types
    varIdType = Dict(variableLabels .=> variableTypeStrings)

    variables = getVariables(client, context; detail=SUMMARY) |> fetch
    for v in variables
        @test v["variableType"] == varIdType[v["label"]]
    end
end

function runVariableTests(client, context)
    @testset "Testing Variables" begin
        @info "Running variable-testset"

        variableLabels = ["x0", "x1", "l1", "x3"]
        variableTypes = [:Pose2,"RoME.Pose2", :Point2, :ContinuousScalar]
        variableTypeStrings = ["RoME.Pose2","RoME.Pose2", "RoME.Point2", "IncrementalInference.Position{1}"]

        @testset "Adding" begin testAddVariable(client, context, variableLabels, variableTypes, variableTypeStrings); end
        @testset "Listing" begin testLs(client, context, variableLabels, variableTypes, variableTypeStrings); end
        @testset "Getting" begin testGetVariable(client, context, variableLabels, variableTypes, variableTypeStrings); end
        @testset "Getting Lists" begin testGetVariables(client, context, variableLabels, variableTypes, variableTypeStrings); end
    end
end


text:
function pdfmerge(output, args...)
    run(`pdfmerge $output $args`)
end

function ffmpeg_imgs(in_pattern, in_fps,
                     out_name, out_fps = 30,
                     vcodec = "h264", pix_fmt = "yuv420p",
                     assume_yes = true)
    run(`ffmpeg $(assume_yes ? "-y" : "") -framerate $in_fps -i $in_pattern -vcodec $vcodec -r $out_fps -pix_fmt $pix_fmt $out_name`)
end

export pdfmerge, ffmpeg_imgs


text:
classes = Dict("Introduction to Real Analysis" => 104, "Programming for Mathematical Applications" => 124 )

get(classes, "Programming for Mathematical Applications", "Class not found")

get(classes, "Introduction to Abstract Algebra", "Class not found")

classes["Programming for Mathematical Applications"]

push!(classes, "Introduction to Abstract Algebra" => 113)

classes["Introduction to Complex Analysis"] = 185

delete!(classes, "Introduction to Real Analysis")


text:
# Autogenerated wrapper script for LLVMBootstrap_jll for x86_64-linux-musl
export libLLVM, libLTO, libclang, llvm_config

using Zlib_jll
using XML2_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libLLVM`
const libLLVM_splitpath = ["lib", "libLLVM.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libLLVM_path = ""

# libLLVM-specific global declaration
# This will be filled out by __init__()
libLLVM_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libLLVM = "libLLVM-8.so"

# Relative path to `libLTO`
const libLTO_splitpath = ["lib", "libLTO.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libLTO_path = ""

# libLTO-specific global declaration
# This will be filled out by __init__()
libLTO_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libLTO = "libLTO.so.8"

# Relative path to `libclang`
const libclang_splitpath = ["lib", "libclang.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libclang_path = ""

# libclang-specific global declaration
# This will be filled out by __init__()
libclang_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libclang = "libclang.so.8"

# Relative path to `llvm_config`
const llvm_config_splitpath = ["bin", "llvm-config"]

# This will be filled out by __init__() for all products, as it must be done at runtime
llvm_config_path = ""

# llvm_config-specific global declaration
function llvm_config(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(llvm_config_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global prefix = abspath(joinpath(@__DIR__, ".."))

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    append!.(Ref(PATH_list), (Zlib_jll.PATH_list, XML2_jll.PATH_list,))
    append!.(Ref(LIBPATH_list), (Zlib_jll.LIBPATH_list, XML2_jll.LIBPATH_list,))

    global libLLVM_path = abspath(joinpath(artifact"LLVMBootstrap", libLLVM_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libLLVM_handle = dlopen(libLLVM_path)
    push!(LIBPATH_list, dirname(libLLVM_path))

    global libLTO_path = abspath(joinpath(artifact"LLVMBootstrap", libLTO_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libLTO_handle = dlopen(libLTO_path)
    push!(LIBPATH_list, dirname(libLTO_path))

    global libclang_path = abspath(joinpath(artifact"LLVMBootstrap", libclang_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libclang_handle = dlopen(libclang_path)
    push!(LIBPATH_list, dirname(libclang_path))

    global llvm_config_path = abspath(joinpath(artifact"LLVMBootstrap", llvm_config_splitpath...))

    push!(PATH_list, dirname(llvm_config_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
# CREATE LATTICES FROM UNITCELLS
# PUT TOGEHTER PATTERNS OF UNITCELLS PERIODICALLY

# Construction of periodic lattice (in 2D)
function getLatticePeriodic(
        :: Type{L},
        unitcell        :: U,
        extent          :: NTuple{2,Int64}
    ) :: L where {
        D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,2},
        U<:AbstractUnitcell{S,B},
        DL,LLS,LLB,SL<:AbstractSite{LLS,DL},BL<:AbstractBond{LLB,2},
        L<:AbstractLattice{SL,BL,U}
    }

    # extract the cardinal directions of the lattice from the array
    N_a1	= abs(extent[1])       :: Int64
    N_a2	= abs(extent[2])       :: Int64
	N_sites = numSites(unitcell)   :: Int64
	N_bonds = numBonds(unitcell)   :: Int64

    # GENERATE NEW POSITIONS
	site_list  = Array{SL,1}(undef, N_a1*N_a2*N_sites)

	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
	    for a in 1:N_sites
			# calculate the site index
			site_index             = index(i,j,a, N_a1,N_a2,N_sites)
			# calculate the position
			site_position          = point(site(unitcell,a)) .+ ((i-1).*a1(unitcell)) .+ ((j-1).*a2(unitcell))
			# set the site in the list
	        site_list[site_index]  = newSite(SL, site_position, label(site(unitcell,a)))
		end
	end
    end

    # GENERATE NEW CONNECTIONS
	bond_list  = Array{BL,1}(undef, N_a1*N_a2*N_bonds)

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
		# add all connections for unitcell (i,j)
		for b in 1:N_bonds
			# get the respective bond
			current_bond = bond(unitcell,b) :: BL
			# calculate the site index from where the connection goes
			index_from = index(i,j,from(current_bond), N_a1,N_a2,N_sites)
			# calculate the aimed unitcell
			i_to = i + wrap(current_bond)[1]
			j_to = j + wrap(current_bond)[2]
			# check if the connection goes around in a1 direction
            (i_to, offset_a1) = getConstrainedIndex(i_to, N_a1)
			# check if the connection goes around in a2 direction
            (j_to, offset_a2) = getConstrainedIndex(j_to, N_a2)
			# get the site index to where the connection goes
			index_to = index(i_to, j_to, to(current_bond), N_a1,N_a2,N_sites)
            # generate a new connection and set it in the list
			bond_list[index(i,j,b, N_a1,N_a2,N_bonds)] = newBond(BL, index_from, index_to, label(current_bond), (offset_a1, offset_a2))
		end
	end
	end

    # generate new lattice vectors, now spanning the WHOLE lattice
    lattice_vectors = Vector{Float64}[
			# multiply by extent of lattice
			a1(unitcell) .* N_a1,
			a2(unitcell) .* N_a2
		]

    # save everything to a Lattice object
    lattice = newLattice(
            L,
        	lattice_vectors,
        	site_list,
        	bond_list,
            unitcell
        )

    # return the lattice
    return lattice
end

# Construction of periodic lattice (in 3D)
function getLatticePeriodic(
        :: Type{L},
        unitcell        :: U,
        extent          :: NTuple{3,Int64}
    ) :: L where {
        D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,3},
        U<:AbstractUnitcell{S,B},
        DL,LLS,LLB,SL<:AbstractSite{LLS,DL},BL<:AbstractBond{LLB,3},
        L<:AbstractLattice{SL,BL,U}
    }

    # extract the cardinal directions of the lattice from the array
    N_a1	= abs(extent[1])       :: Int64
    N_a2	= abs(extent[2])       :: Int64
    N_a3	= abs(extent[3])       :: Int64
	N_sites = numSites(unitcell)   :: Int64
	N_bonds = numBonds(unitcell)   :: Int64

    # GENERATE NEW POSITIONS
	site_list  = Array{SL,1}(undef, N_a1*N_a2*N_a3*N_sites)

	# set all positions to their correct values
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
	    for a in 1:N_sites
			# calculate the site index
			site_index             = index(i,j,k,a, N_a1,N_a2,N_a3,N_sites)
			# calculate the position
			site_position          = point(site(unitcell,a)) .+ ((i-1).*a1(unitcell)) .+ ((j-1).*a2(unitcell)) .+ ((k-1).*a3(unitcell))
			# set the site in the list
	        site_list[site_index]  = newSite(SL, site_position, label(site(unitcell,a)))
		end
	end
    end
    end

    # GENERATE NEW CONNECTIONS
	bond_list  = Array{BL,1}(undef, N_a1*N_a2*N_a3*N_bonds)

	# iterate over all unit cells
	for i in 1:N_a1
	for j in 1:N_a2
	for k in 1:N_a3
		# add all connections for unitcell (i,j)
		for b in 1:N_bonds
			# get the respective bond
			current_bond = bond(unitcell, b) :: BL
			# calculate the site index from where the connection goes
			index_from = index(i,j,k, from(current_bond), N_a1,N_a2,N_a3,N_sites)
			# calculate the aimed unitcell
			i_to = i + wrap(current_bond)[1]
			j_to = j + wrap(current_bond)[2]
			k_to = k + wrap(current_bond)[3]
			# check if the connection goes around in a1 direction
            (i_to, offset_a1) = getConstrainedIndex(i_to, N_a1)
			# check if the connection goes around in a2 direction
            (j_to, offset_a2) = getConstrainedIndex(j_to, N_a2)
			# check if the connection goes around in a3 direction
            (k_to, offset_a3) = getConstrainedIndex(k_to, N_a3)
			# get the site index to where the connection goes
			index_to = index(i_to,j_to,k_to, to(current_bond), N_a1,N_a2,N_a3,N_sites)
            # generate a new connection and set it in the list
			bond_list[index(i,j,k,b, N_a1,N_a2,N_a3,N_bonds)] = newBond(BL, index_from, index_to, label(current_bond), (offset_a1, offset_a2, offset_a3))
		end
	end
	end
	end

    # generate new lattice vectors, now spanning the WHOLE lattice
    lattice_vectors = Vector{Float64}[
			# multiply by extent of lattice
			a1(unitcell) .* N_a1,
			a2(unitcell) .* N_a2,
			a3(unitcell) .* N_a3
		]

    # save everything to a Lattice object
    lattice = newLattice(
            L,
        	lattice_vectors,
        	site_list,
        	bond_list,
            unitcell
        )

    # return the lattice
    return lattice
end

# Construction wrapper for default type
function getLatticePeriodic(
        unitcell        :: U,
        extent          :: NTuple{N,Int64}
    ) :: Lattice{S,B,U} where {
        N,D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},
        U<:AbstractUnitcell{S,B}
    }

    # return the general function
    return getLatticePeriodic(Lattice{S,B,U}, unitcell, extent)
end

# short hand notation (only pass N instead of (N1,N2,N3))

# 1d (General)
function getLatticePeriodic(
        :: Type{L},
        unitcell        :: U,
        extent          :: Int64
    ) :: L where {
        D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,1},
        U<:AbstractUnitcell{S,B},
        DL,LLS,LLB,SL<:AbstractSite{LLS,DL},BL<:AbstractBond{LLB,1},
        L<:AbstractLattice{SL,BL,U}
    }

    # return the suitable function
    return getLatticePeriodic(L, unitcell, (extent, ))
end
# 1d (Wrapper)
function getLatticePeriodic(
        unitcell        :: U,
        extent          :: Int64
    ) :: Lattice{S,B,U} where {D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,1},U<:AbstractUnitcell{S,B}}

    # return the suitable function
    return getLatticePeriodic(Lattice{S,B,U}, unitcell, (extent, ))
end

# 2d
function getLatticePeriodic(
        :: Type{L},
        unitcell        :: U,
        extent          :: Int64
    ) :: L where {
        D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,2},
        U<:AbstractUnitcell{S,B},
        DL,LLS,LLB,SL<:AbstractSite{LLS,DL},BL<:AbstractBond{LLB,2},
        L<:AbstractLattice{SL,BL,U}
    }

    # return the suitable function
    return getLatticePeriodic(L, unitcell, (extent, extent))
end
# 2d (Wrapper)
function getLatticePeriodic(
        unitcell        :: U,
        extent          :: Int64
    ) :: Lattice{S,B,U} where {D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,2},U<:AbstractUnitcell{S,B}}

    # return the suitable function
    return getLatticePeriodic(Lattice{S,B,U}, unitcell, (extent, extent))
end

# 3d
function getLatticePeriodic(
        :: Type{L},
        unitcell        :: U,
        extent          :: Int64
    ) :: L where {
        D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,3},
        U<:AbstractUnitcell{S,B},
        DL,LLS,LLB,SL<:AbstractSite{LLS,DL},BL<:AbstractBond{LLB,3},
        L<:AbstractLattice{SL,BL,U}
    }

    # return the suitable function
    return getLatticePeriodic(L, unitcell, (extent, extent, extent))
end
# 3d (Wrapper)
function getLatticePeriodic(
        unitcell        :: U,
        extent          :: Int64
    ) :: Lattice{S,B,U} where {D,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,3},U<:AbstractUnitcell{S,B}}

    # return the suitable function
    return getLatticePeriodic(Lattice{S,B,U}, unitcell, (extent, extent, extent))
end

# export the functions
export getLatticePeriodic


text:
""" Functions indicating if a subproblem can be fathomed and giving information on the dominated nadir points of the lower bound and upper Bound solutions """

# Public functions :
# - SPS_partially_implicit : study sub problem using partial dichotomy

#-------------------------------------------------------------------#
#------------------------ PUBLIC FUNCTIONS -------------------------#
#-------------------------------------------------------------------#

"""
Optimized pruning test (for dichoSearch type UB computation)

Computes the UB of the subProblem and stops as soon as the sub problem cannot be fathomed.
"""
function SPS_partially_implicit(prob::BiOKP, assignment::Assignment, nadirsToStudy::Vector{PairOfSolution}; parentToChild = ParentToChild())

    CONFIG.debug && DEBUG_nadirs(prob,nadirsToStudy)
    CONFIG.debug && DEBUG_parentToChild(prob,parentToChild)

    UB = UpperBound() # init

    # get initial solutions
    UB.sols = copy(parentToChild.knownSols)

    # compute lexicographic optimal solutions
    leftSol, rightSol = getLexicoSolutions(prob, assignment, UB.sols, getSolver1OKP(COMPONENTS.methodUB), parentToChild = parentToChild)

    # In the event that the two solutions found are identical, this solution is optimal and thus added to the lower bound
    if !differentFloatSols(leftSol,rightSol)
        if UB.sols == nil(Sol) # if the UB.sols contains already one solution, we don't need to add this one because it is the same
            UB.sols = cons(leftSol, UB.sols)
        end
        # we don't have any pairs to study

        if leftSol.isBinary
            CONFIG.debug && DEBUG_UB(prob, UB)
            return OPTIMALITY, nadirsToStudy, Vector{PairOfSolution}(), UB
        else
            CONFIG.debug && DEBUG_UB(prob, UB)
            return NONE, nadirsToStudy, Vector{PairOfSolution}(), UB
        end

    else # the two lexicographic solutions are different

        # get the list of pairs of solutions to study
        toStudy, UB.sols = computesPairsToStudy(prob, UB.sols, parentToChild, leftSol, rightSol)

        # add the constraint extreme horizontal and vertical of the polytope
        A = Matrix{Float64}(undef,2,2)
        b = Vector{Float64}(undef,2)

        # vertical (on the right)
        A[1, 1] = 1
        A[1, 2] = 0

        b[1] = rightSol.y[1]

        # horizontal (on the left)
        A[2, 1] = 0
        A[2, 2] = 1

        b[2] = leftSol.y[2]

        subDominated, nadirsToStudy, nonDomNadirs = PRUNING_compareCstrToNadirs(prob, DualSet(A,b), nadirsToStudy)

        # while some pairs of solutions are found by the dichotomic search, we keep going
        while toStudy != [] && !subDominated

            currPair = pop!(toStudy)

            leftSol = currPair.solL
            rightSol = currPair.solR

            λ = [leftSol.y[2] - rightSol.y[2], rightSol.y[1] - leftSol.y[1]]

            solver1OKP = getSolver1OKP(COMPONENTS.methodUB)
            midSol     = solver1OKP(prob, λ, assignment)

            # if the solution is beyond the pair studied, it's added to the UB
            if round(sum(λ .* midSol.y), digits = 7) - round(sum(λ .* rightSol.y), digits = 7) > 0.00001
                # add the new solution found
                strictInsertInFront!(prob, UB.sols, midSol)
                # add the two newly found pairs
                push!(toStudy, PairOfSolution(leftSol, midSol))
                push!(toStudy, PairOfSolution(midSol, rightSol))
            end

            # we construct the constraint associated with the new solution
            b = (leftSol.y[2] - rightSol.y[2]) * midSol.y[1] + (rightSol.y[1] - leftSol.y[1]) * midSol.y[2]

            A = Matrix{Float64}(undef,1,2)
            A[1,1] = λ[1]
            A[1,2] = λ[2]

            subDominated, nadirsToStudy, nonDomNadirs = PRUNING_compareCstrToNadirs(prob, DualSet(A,[b]),nadirsToStudy, nonDomNadirs = nonDomNadirs)
        end

        if subDominated

            CONFIG.debug && DEBUG_nadirs(prob,nadirsToStudy)
            CONFIG.debug && DEBUG_UB(prob, UB)

            return DOMINANCE, Vector{PairOfSolution}(), nadirsToStudy, UB

        else

            CONFIG.debug && DEBUG_nadirs(prob, nadirsToStudy)
            CONFIG.debug && DEBUG_nadirs(prob, nonDomNadirs)
            CONFIG.debug && DEBUG_UB(prob, UB)

            return NONE, nadirsToStudy, nonDomNadirs, UB

        end
    end
end

#-------------------------------------------------------------------#
#----------------------- PRIVATE FUNCTIONS -------------------------#
#-------------------------------------------------------------------#

"""
Indicates if the subproblem is dominated, as well as a partition of nadir points (dominated vs non dominated)

[ATTENTION] For nadir points, "dominated" doesn't mean "discarded". It is the OPPOSITE ! If the subproblem isn't dominated, we will keep studying the dominated nadirs, not the nondominated ones. In the rest of the code, dominated nadirs are therefore called "nadirsToStudy".

Returns isSubDominated, dominatedNadirs, nonDomNadirs
"""
function PRUNING_compareCstrToNadirs(prob::BiOKP, constraint::DualSet, nadirsToStudy::Vector{PairOfSolution}; nonDomNadirs::Vector{PairOfSolution} = Vector{PairOfSolution}(), regionLimit::Vector{Float64} = [-1.,-1.])

    CONFIG.debug && DEBUG_nadirs(prob,nadirsToStudy)
    CONFIG.debug && DEBUG_nadirs(prob,nonDomNadirs)

    dominatedNadirs = Vector{PairOfSolution}()

    regionLimit == [-1,-1] ? checkLimits = false : checkLimits = true

    A = constraint.A
    b = constraint.b

    isSubDominated = true
    iterNadirs = 0
    # for each nadir point to study
    while iterNadirs < length(nadirsToStudy)

        iterNadirs += 1

        pairNadir = nadirsToStudy[iterNadirs]

        if prob.isInteger && COMPONENTS.nadirsShift
		    nadir = [pairNadir.solL.y[1] + 1, pairNadir.solR.y[2] + 1] # - 0.0...1 to avoid bad roundings
        else
            nadir = [pairNadir.solL.y[1], pairNadir.solR.y[2]]
        end

        # if the nadir is out of the bounds
        if checkLimits && (nadir[1] > regionLimit[1] || nadir[2] > regionLimit[2])
            # the nadir is non dominated
            push!(nonDomNadirs,pairNadir)

        else

            nadirA = A * nadir

            isNadirUnderUB = true
            iterConstraint = 1
            while iterConstraint <= length(b) && isNadirUnderUB
                # compute isNadirUnderUB : verify that the nadir point is lower that the upper bound on all objectives
                if prob.isInteger && COMPONENTS.nadirsShift
                    isNadirUnderUB = isNadirUnderUB && nadirA[iterConstraint] <= b[iterConstraint]
                else
                    isNadirUnderUB = isNadirUnderUB && nadirA[iterConstraint] < b[iterConstraint]
                end
                iterConstraint += 1
            end

            if isNadirUnderUB # the nadir is dominated but the UB
                push!(dominatedNadirs,pairNadir)
                isSubDominated = false
            else
                push!(nonDomNadirs,pairNadir)
            end
        end
	end

    while iterNadirs < length(nadirsToStudy)
        iterNadirs += 1
        push!(dominatedNadirs,nadirsToStudy[iterNadirs])
    end

    sort!(nonDomNadirs, lt = (x, y) -> isless(x.solL.y[1],y.solL.y[1]))

    CONFIG.debug && DEBUG_nadirs(prob, dominatedNadirs)
    CONFIG.debug && DEBUG_nadirs(prob, nonDomNadirs)

    return isSubDominated, dominatedNadirs, nonDomNadirs
end


text:
using Pkg
Pkg.instantiate()
Pkg.develop(PackageSpec(path=pwd()))
Pkg.build("Swalbe")


text:
using ScenTrees
using Test
using Statistics: std
using CSV, DataFrames, XLSX

@testset "ScenTrees.jl" begin
    @testset "Predefined tree - Tree 402" begin
        a = Tree(402)
        @test typeof(a) == Tree
        @test length(a.parent) == 15
        @test length(a.state) == length(a.probability) == length(a.parent) == 15
        @test sum(a.probability) == 8.0
        @test length(a.children) == 8
        @test length(root(a)) == 1
    end
    @testset "Initial Trees" begin
        init = Tree([1,2,2,2],1)
        @test typeof(init) == Tree
        @test length(init.parent) == 15
        @test length(init.state) == length(init.probability) == length(init.parent) == 15
        @test length(init.children) == 8
        @test length(stage(init)) == 15
        @test height(init) == 3
        @test length(leaves(init)) == 3
        @test nodes(init) == 1:15
        @test length(nodes(init)) == 15
        @test length(root(init)) == 1
    end
    @testset "A sample of a Scenario Tree 1D" begin
        x = Tree([1,3,3,3,3],1)
        @test typeof(x) == Tree
        @test length(x.parent) == 121
        @test length(x.state) == length(x.probability) == length(x.parent) == 121
        @test sum(x.probability) ≈ 41.0
        @test length(x.children) == 41
        @test length(root(x)) == 1
        @test length(leaves(x)) == 3
    end
    @testset "A sample of a Scenario Tree 2D" begin
        y = Tree([1,3,3,3,3],2)
        @test typeof(y) == Tree
        @test length(y.parent) == 121
        @test length(y.probability) == length(y.parent) == 121
        @test length(y.state) == length(y.parent)*2 == 242
        @test sum(y.probability) ≈ 41.0
        @test length(y.children) == 41
    end
    @testset "Sample stochastic functions" begin
        a = gaussian_path1D()
        b = running_maximum1D()
        c = path()
        d = gaussian_path2D()
        e = running_maximum2D()
        @test length(a) == 4
        @test length(b) == 4
        @test length(c) == 4
        @test size(d) == (4,2)
        @test size(e) == (4,2)
    end
    @testset "ScenTrees.jl - Tree Approximation 1D" begin
        paths = [gaussian_path1D,running_maximum1D]
        trees = [Tree([1,2,2,2]),Tree([1,3,3,3])]
        samplesize = 100000
        p = 2
        r = 2

        for path in paths
            for newtree in trees
                tree_approximation!(newtree,path,samplesize,p,r)
                @test length(newtree.parent) == length(newtree.state)
                @test length(newtree.parent) == length(newtree.probability)
                @test length(stage(newtree)) == length(newtree.parent)
                @test height(newtree) == maximum(stage(newtree))
                @test round(sum(leaves(newtree)[3]),digits=1) == 1.0   #sum of unconditional probabilities of the leaves
                @test length(root(newtree)) == 1
            end
        end
    end
    @testset "ScenTrees.jl - Tree Approximation 2D" begin
        twoD = tree_approximation!(Tree([1,3,3,3],2),gaussian_path2D,100000,2,2)
        @test size(twoD.state,2) == 2
        @test size(twoD.state,1) == length(twoD.parent) == length(twoD.probability)
    end

    @testset "ScenTrees.jl - Lattice Approximation" begin
        tstLat = lattice_approximation([1,2,3,4],gaussian_path1D,500000,2,1)
        @test length(tstLat.state) == length(tstLat.probability)
        @test round.(sum.(tstLat.probability), digits = 1)  == [1.0, 1.0, 1.0, 1.0] #sum of probs at every stage
    end

    @testset "ScenTrees.jl - Lattice Approximation 2D" begin
        lat2 = lattice_approximation([1,2,3,4],gaussian_path2D,500000,2,2)
        @test length(lat2) == 2 # resultant lattices are 2
        @test length(lat2[1].state) == length(lat2[1].probability)
        @test round.(sum.(lat2[1].probability), digits = 1)  == [1.0, 1.0, 1.0, 1.0]
        @test round.(sum.(lat2[2].probability), digits = 1)  == [1.0, 1.0, 1.0, 1.0]
    end
    @testset "ScenTrees.jl - Test Example Data" begin
        data = CSV.read("data5.csv",DataFrame) # check the new way of loading data
        gsData = Matrix(data)
        df = DataFrame(XLSX.readtable("Mappe1.xlsx", "Sheet1")...)
        df1 = Matrix(df)
        df1 = convert(Array{Float64,2},df1)
        df2 = DataFrame(XLSX.readtable("Mappe1.xlsx", "Sheet2")...)
        df22 = Matrix(df2)
        df22 = convert(Array{Float64,2},df22)
        RandomWalkData = CSV.read("RandomDataWalk.csv",DataFrame)
        RWData = Matrix(RandomWalkData)

        @test size(gsData) ==(100,10)
        @test size(df1) == (1000,7)
        @test size(df22) == (1000,4)
        @test size(RWData) == (1000,5)

        sd1 = zeros(5)
        sd2 = zeros(7)

        for t = 1 : size(RWData,2)
            sd1[t] = std(RWData[:,t])
        end

        for t = 1:size(df1,2)
            sd2[t] = std(df1[:,t])
        end
        @test (sd1 .< 5) == Bool[true, true, true, true, true]
        @test (sd2 .< 10) == Bool[true, true, true, true, true, true, true]

        LatFromKernel = lattice_approximation([1,3,4,5,6],kernel_scenarios(RWData),100000,2,1)
        @test round.(sum.(LatFromKernel.probability),digits=1) == [1.0, 1.0, 1.0, 1.0, 1.0]
        @test length(LatFromKernel.state) == length(LatFromKernel.probability)
    end
end


text:
using Test
using DivPopGa

include("testclusteredga.jl")
include("testsummary.jl")

max_stars_repo_name: akadal/divpopga

max_issues_repo_name: akadal/divpopga

max_forks_repo_name: akadal/divpopga


text:
using Pkg
add("https://github.com/snotskie/EpistemicNetworkAnalysis.jl")
add("DataFrames")
add("CSV")
add("Plots")
add("GLM")
add("Statistics")
add("MixedModels")
add("HypothesisTests")


text:
## A discrete distribution representing a point mass (DiracPM) or several point masses (MultvariateDiracPM)
## Like Distributions.Categorical but remembers the posiition of the point mass
## This is needed to specify the priors under the null and alternative hypotheses

struct DiracPM <: DiscreteUnivariateDistribution
    x::Float64
end

struct MultivariateDiracPM <: DiscreteMultivariateDistribution
    x::Vector{Float64}
end

import Statistics.var, Statistics.mean, Base.rand, Distributions.cdf, Base.length

mean(d::DiracPM) = d.x
var(d::DiracPM) = 0.0
rand(d::DiracPM) = d.x
rand(d::DiracPM, B::Int) = repeat([d.x], inner=B)
cdf(d::DiracPM, x::Float64) = Float64(d.x <= x)
length(d::DiracPM) = 1

mean(d::MultivariateDiracPM) = d.x
var(d::MultivariateDiracPM) = 0.0
rand(d::MultivariateDiracPM) = d.x
rand(d::MultivariateDiracPM, B::Int) = reshape(repeat(d.x, outer=B), length(d.x), B)
cdf(d::MultivariateDiracPM, x::Vector{Float64}) = Float64(all(d.x .<= x))
length(d::MultivariateDiracPM) = length(d.x)


text:
module Jaynes

# Yarrrr I'm a com-pirate!
using IRTools
using IRTools: @dynamo, IR, xcall, arguments, insertafter!, recurse!, isexpr, self, argument!, Variable
#using Mjolnir
#using Mjolnir: Basic, AType, Const, abstract, Multi, @abstract, Partial, trace
#using Mjolnir: Defaults
using MacroTools

using Reexport
@reexport using Distributions

# Differentiable.
using Zygote
using DistributionsAD
using Flux.Optimise: update!
@reexport using Flux.Optimise

# Toplevel importants :)
const Address = Union{Int, Symbol, Pair}

# ------------ Com-pirate fixes ------------ #

# TODO: This chunk below me is currently required to fix an unknown performance issue in Base. Don't be alarmed if this suddenly disappears in future versions.
unwrap(gr::GlobalRef) = gr.name
unwrap(gr) = gr

# Whitelist includes vectorized calls.
whitelist = [
             # Base
             :rand, :_apply_iterate,

             # Specialized call sites
             :learnable, :markov, :plate, :cond, :factor,

             # Foreign model interfaces
             :foreign
            ]

# Fix for specialized tracing.
function recur!(ir, to = self)
    for (x, st) in ir
        isexpr(st.expr, :call) && begin
            ref = unwrap(st.expr.args[1])
            ref in whitelist || continue
            ir[x] = Expr(:call, to, st.expr.args...)
        end
    end
    return ir
end

# Fix for _apply_iterate.
function f_push!(arr::Array, t::Tuple{}) end
f_push!(arr::Array, t::Array) = append!(arr, t)
f_push!(arr::Array, t::Tuple) = append!(arr, t)
f_push!(arr, t) = push!(arr, t)
function flatten(t::Tuple)
    arr = Any[]
    for sub in t
        f_push!(arr, sub)
    end
    return arr
end

# ------------ includes ------------ #

include("compiler/static.jl")
include("traces.jl")
include("selections.jl")
include("utils/numerical.jl")
include("utils/vectorized.jl")
include("compiler/diffs.jl")
include("contexts.jl")
include("inference.jl")
include("foreign_model_interfaces.jl")

# Contexts.
export Generate, generate
export Simulate, simulate
export Update, update
export Propose, propose
export Regenerate, regenerate
export Score, score
export Backpropagate, get_parameter_gradients, get_choice_gradients, train

# Tracer language features.
export learnable, plate, markov, cond, factor

# Diffs.
export NoChange, UndefinedChange, VectorDiff

# Selections and parameters.
export selection, array, parameters
export anywhere, intersection, union
export get_selection, compare, has_top, update_parameters, dump_queries, merge!, merge

# Inference.
export importance_sampling, is
export initialize_filter, filter_step!, check_ess_resample!, get_lmle, pf
export metropolis_hastings, mh
export automatic_differentiation_variational_inference, advi
export hamiltonian_monte_carlo, hmc

# Foreign model interfaces.
export @primitive
export @load_gen_fmi
export @load_soss_fmi
export @load_flux_fmi
export foreign

# Utilities.
export display, getindex, haskey, get_score, get_ret

end # module


text:
# Computer Number Guessing Game in Julia

using Random

function play_number_guess_computer()

    print("Please enter a number between 1 and 50 for the computer to try and guess: ")

    # Take in the user input and convert it to a number
    target_number = parse(Int64, readline())

    # Create an array of 50 numbers
    guess_order = collect(1:50)

    # Define our random seed
    rng = MersenneTwister(1234)

    # Shuffle the array randomly given ur seed
    shuffled_guess = shuffle(rng, guess_order)

    # Loop through each guess and see if it right
    for guess in shuffled_guess

        if guess == target_number
            print("\nThe computer cracked the code and guessed it right!")
            break # Stop the for loop if we get it right
        end

        print("\nComputer guessed: $guess")
    end
end


text:
function test_add_remove(storage::IS.TimeSeriesStorage)
    name = "component1"
    label = "val"
    component = IS.TestComponent(name, 5)
    ts = create_time_series_data()
    IS.add_time_series!(storage, IS.get_uuid(component), label, ts)

    ts_data = IS.get_time_series(storage, IS.get_uuid(ts))

    @test TimeSeries.timestamp(ts_data) == TimeSeries.timestamp(ts.data)
    @test TimeSeries.values(ts_data) == TimeSeries.values(ts.data)

    component2 = IS.TestComponent("component2", 6)
    IS.add_time_series!(storage, IS.get_uuid(component2), label, ts)

    @test IS.get_num_time_series(storage) == 1

    IS.remove_time_series!(storage, IS.get_uuid(ts), IS.get_uuid(component2), label)

    # There should still be one reference to the data.
    ts_data2 = IS.get_time_series(storage, IS.get_uuid(ts))
    @test ts_data2 isa TimeSeries.TimeArray

    IS.remove_time_series!(storage, IS.get_uuid(ts), IS.get_uuid(component), label)
    @test_throws ArgumentError IS.get_time_series(storage, IS.get_uuid(ts))
    IS.get_num_time_series(storage) == 0
end

function test_add_references(storage::IS.TimeSeriesStorage)
    label = "val"
    component1 = IS.TestComponent("component1", 5)
    component2 = IS.TestComponent("component2", 6)
    ts = create_time_series_data()
    ts_uuid = IS.get_uuid(ts)
    IS.add_time_series!(storage, IS.get_uuid(component1), label, ts)
    IS.add_time_series_reference!(storage, IS.get_uuid(component2), label, ts_uuid)

    @test IS.get_num_time_series(storage) == 1

    IS.remove_time_series!(storage, ts_uuid, IS.get_uuid(component1), label)

    # There should still be one reference to the data.
    @test IS.get_time_series(storage, ts_uuid) isa TimeSeries.TimeArray

    IS.remove_time_series!(storage, ts_uuid, IS.get_uuid(component2), label)
    @test_throws ArgumentError IS.get_time_series(storage, IS.get_uuid(ts))
    IS.get_num_time_series(storage) == 0
end

function test_get_subset(storage::IS.TimeSeriesStorage)
    name = "component1"
    label = "val"
    component = IS.TestComponent(name, 1)
    ts = create_time_series_data()
    IS.add_time_series!(storage, IS.get_uuid(component), label, ts)
    ts_data = IS.get_time_series(storage, IS.get_uuid(ts))

    @test TimeSeries.timestamp(ts_data) == TimeSeries.timestamp(ts.data)
    index = 3
    len = 5
    ts_subset = IS.get_time_series(storage, IS.get_uuid(ts); index = index, len = len)
    @test ts_subset[1] == ts_data[index]
    @test length(ts_subset) == len
end

function test_clear(storage::IS.TimeSeriesStorage)
    name = "component1"
    label = "val"
    component = IS.TestComponent(name, 5)
    ts = create_time_series_data()
    IS.add_time_series!(storage, IS.get_uuid(component), label, ts)

    ts_data = IS.get_time_series(storage, IS.get_uuid(ts))

    @test TimeSeries.timestamp(ts_data) == TimeSeries.timestamp(ts.data)
    @test TimeSeries.values(ts_data) == TimeSeries.values(ts.data)

    IS.clear_time_series!(storage)
    @test_throws ArgumentError IS.get_time_series(storage, IS.get_uuid(ts))
end

@testset "Test time series storage implementations" begin
    for in_memory in (true, false)
        test_add_remove(IS.make_time_series_storage(; in_memory = in_memory))
        test_get_subset(IS.make_time_series_storage(; in_memory = in_memory))
        test_clear(IS.make_time_series_storage(; in_memory = in_memory))
    end

    test_add_remove(IS.make_time_series_storage(; in_memory = false, directory = "."))
    test_get_subset(IS.make_time_series_storage(; in_memory = false, directory = "."))
    test_clear(IS.make_time_series_storage(; in_memory = false, directory = "."))
end

@testset "Test copy time series references" begin
    for in_memory in (true, false)
        test_add_remove(IS.make_time_series_storage(; in_memory = in_memory))
        test_add_references(IS.make_time_series_storage(; in_memory = in_memory))
        test_get_subset(IS.make_time_series_storage(; in_memory = in_memory))
        test_clear(IS.make_time_series_storage(; in_memory = in_memory))
    end
end


text:
using CoordinateTransformations
using Base.Test
using ForwardDiff: Dual, partials
using StaticArrays

@testset "CoordinateTransformations" begin

    include("core.jl")
    include("coordinatesystems.jl")
    include("affine.jl")

end


text:
using Documenter
using OffsetArrays

makedocs(
    sitename = "OffsetArrays",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    pages = ["index.md", "internals.md", "reference.md"],
    modules = [OffsetArrays],
    doctestfilters = [r"at \./.*", r"at /home.*", r"top-level scope.*", r"\[\d*\]\s*$"],   # for backtraces
)

deploydocs(
    repo = "github.com:JuliaArrays/OffsetArrays.jl.git"
)


text:
# Julia version : 1.3.1
# Flux version : v0.10.1

using Random
using Dates
using CuArrays
using CUDAdrv
using CUDAnative: device!
using Flux, Metalhead, Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Metalhead: trainimgs
using Images: channelview
using Statistics: mean
using Base.Iterators: partition

model_file = joinpath(dirname(@__FILE__),"cifar10_vgg16_model.bson")

# Get arguments

epochs = 100
batch_size = 128
gpu_device = 0

# Very important : this prevent loss NaN
ϵ = 1.0f-32

# use 1nd GPU
#CUDAnative.device!(0)
device!(gpu_device)
CuArrays.allowscalar(false)

# VGG16 and VGG19 models
vgg16() = Chain(
  Conv((3, 3), 3 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  MaxPool((2,2)),
  Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  MaxPool((2,2)),
  Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  MaxPool((2,2)),
  Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  MaxPool((2,2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  MaxPool((2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(512, 4096, relu),
  Dropout(0.5),
  Dense(4096, 4096, relu),
  Dropout(0.5),
  Dense(4096, 10),
  softmax)

vgg19() = Chain(
  Conv((3, 3), 3 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(64),
  MaxPool((2,2)),
  Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(128),
  MaxPool((2,2)),
  Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(256),
  Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  BatchNorm(512),
  Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
  MaxPool((2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(512, 4096, relu),
  Dropout(0.5),
  Dense(4096, 4096, relu),
  Dropout(0.5),
  Dense(4096, 10),
  softmax)

m = vgg16() |> gpu

# Function to convert the RGB image to Float32 Arrays
getarray(X) = Float32.(permutedims(channelview(X), (2, 3, 1)))

function make_minibatch(imgs,labels,batch_size)
  data_set = [(cat(imgs[i]..., dims = 4),
          labels[:,i])
          for i in partition(1:length(imgs), batch_size)]
  return data_set
end

X = trainimgs(CIFAR10)
train_idxs = 1:49000
train_imgs = [getarray(X[i].img) for i in train_idxs]
train_labels = float.(onehotbatch([X[i].ground_truth.class for i in train_idxs],1:10))
train_set = make_minibatch(train_imgs,train_labels,batch_size)

verify_idxs = 49001:50000
verify_imgs = cat([getarray(X[i].img) for i in verify_idxs]..., dims = 4)
verify_labels = float.(onehotbatch([X[i].ground_truth.class for i in verify_idxs],1:10))
verify_set = [(verify_imgs,verify_labels)]

# Fetch the test data from Metalhead and get it into proper shape.
# CIFAR-10 does not specify a verify set so valimgs fetch the testdata instead of testimgs
tX = valimgs(CIFAR10)
test_idxs = 1:10000
test_imgs = [getarray(tX[i].img) for i in test_idxs]
test_labels = float.(onehotbatch([tX[i].ground_truth.class for i in test_idxs], 1:10))
test_set = make_minibatch(test_imgs,test_labels,batch_size)

# Defining the loss and accuracy functions
loss(x, y) = crossentropy(m(x) .+ ϵ, y)

function accuracy(data_set)
  batch_size = size(data_set[1][1])[end]
  l = length(data_set)*batch_size
  s = 0f0
  for (x,y) in data_set
    s += sum((onecold(m(x|>gpu) |> cpu) .== onecold(y|>cpu)))
  end
  return s/l
end

# Make sure our is nicely precompiled befor starting our training loop
m(train_set[1][1] |> gpu)

# Defining the callback and the optimizer
opt = ADAM(0.001)

@info "Training model..."

for epoch_idx in 1:epochs
  accs = Array{Float32}(undef,0)

  train_set_len = length(train_set)
  shuffle_idxs = collect(1:train_set_len)
  shuffle!(shuffle_idxs)

  for (idx,data_idx) in enumerate(shuffle_idxs)
    (x,y) = train_set[data_idx]
    # We augment `x` a little bit here, adding in random noise
    x = (x .+ ϵ*randn(eltype(x),size(x))) |> gpu
    y = y|> gpu
    Flux.train!(loss,params(m),[(x,y)],opt)
    v_acc = accuracy(verify_set)
    @info "Epoch# $(epoch_idx)/$(epochs) - #$(idx)/$(train_set_len) loss: $(loss(x,y)), accuracy: $(v_acc)"
    push!(accs,v_acc)
  end
  m_acc = mean(accs)
  @info " -> Verify accuracy(mean) : $(m_acc)"
end
test_acc = accuracy(test_set)
@info "Test accuracy : $(test_acc)"

max_stars_repo_name: mrchaos/model-zoo

max_issues_repo_name: mrchaos/model-zoo

max_forks_repo_name: mrchaos/model-zoo


text:
using ADCME
using Libdl

#--------------------------------------------------------------
# install gmsh
const prefix = joinpath(ADCME.LIBDIR,"..")

download_info = Dict(
    "Linux64" => ("https://gmsh.info/bin/Linux/gmsh-4.5.6-Linux64-sdk.tgz", "gmsh.tgz"),
    "Windows64" => ("https://gmsh.info/bin/Windows/gmsh-4.5.6-Windows64-sdk.zip", "gmsh.zip"),
    "Windows32" => ("https://gmsh.info/bin/Windows/gmsh-4.5.6-Windows32-sdk.zip", "gmsh.zip"),
    "MacOSX" => ("https://gmsh.info/bin/MacOSX/gmsh-4.5.6-MacOSX-sdk.tgz", "gmsh.tgz"),
    "Linux32" => ("https://gmsh.info/bin/Linux/gmsh-4.5.6-Linux32-sdk.tgz", "gmsh.tgz")
)

key = ""
if Sys.iswindows()
    global key = "Windows$(Sys.WORD_SIZE)"
elseif Sys.isapple()
    global key = "MacOSX"
else
    global key = "Linux$(Sys.WORD_SIZE)"
end

if !isfile(joinpath(prefix, download_info[key][2]))
    download(download_info[key][1], joinpath(prefix, download_info[key][2]))
    if !isdir(joinpath(prefix, "gmsh"))
        mkdir(joinpath(prefix, "gmsh"))
    end
    if Sys.iswindows()
        run(`cmd /c unzip $(joinpath(prefix, download_info[key][2])) -d $(joinpath(prefix, "gmsh"))`)
    else
        run(`tar zxvf $(joinpath(prefix, download_info[key][2])) -C $(joinpath(prefix, "gmsh"))`)
    end
end

cp(joinpath(prefix, "gmsh", "gmsh-4.5.6-$key-sdk", "lib"), joinpath(@__DIR__, "Gmsh"), force=true)
cp(joinpath(@__DIR__, "gmsh.jl"), joinpath(@__DIR__, "Gmsh", "gmsh.jl"), force=true)

install_adept()

files = readdir("CustomOp")
all_files = []
for f in files
    if occursin("CMake", f) || occursin("build", f)
        continue
    end
    push!(all_files, f)
end

srcs = "add_library(DataLib SHARED "
for r in all_files
    fs = readdir("CustomOp/"*r)
    fs = filter(x->endswith(x, ".cpp"), fs)[1]
    global srcs *= "\${CMAKE_SOURCE_DIR}/$r/$fs "
end
srcs *= ")"

cnt = String(read("CustomOp/CMakeListsTemplate.txt"))
cnt = replace(cnt, "[REPLACE]"=>srcs)
open("CustomOp/CMakeLists.txt", "w") do io
    write(io, cnt)
end

change_directory("$(@__DIR__)/CustomOp/build")
ADCME.cmake()
ADCME.make()
cd("../../")