text:
using AmbulanceDeployment, DataFrames, JLD, Distributions

turnaround = Distributions.LogNormal(3.65, 0.3)
ncalls = 1000
namb = 30
lambda = 0

hourly_calls = DataFrames.readtable("data/processed/2-weekday_calls.csv")
# weekend_hourly_calls = DataFrames.readtable("data/processed/2-weekend_calls.csv")
adjacent_nbhd = DataFrames.readtable("data/processed/2-adjacent_nbhd.csv")
coverage = JLD.load("data/processed/3-coverage.jld", "stn_coverage")
hospitals = DataFrames.readtable("data/processed/3-hospitals.csv");
stations = DataFrames.readtable("data/processed/3-stations.csv");

solverstats = JLD.load("data/processed/4-solve-stats.jld")
amb_deployment = solverstats["amb_deployment"]
const model_names = (:Stochastic, :Robust01, :Robust005, :Robust001, :Robust0001, :Robust00001, :MEXCLP, :MALP)
name = model_names[1]
x = amb_deployment[name][namb]

p = DeploymentProblem(
    hourly_calls,
    adjacent_nbhd,
    coverage,
    namb = namb,
    train_filter = (hourly_calls[:year] .== 2012) .* (hourly_calls[:month] .<= 3)
)

# We focus on emergency calls during the "peak period" (8AM - 8PM),
# with the emergency calls from the first 3 month as our training set,
# and the subsequent emergency calls from the remaining months as our test set

# calls = DataFrames.readtable("data/processed/5-calls.csv");
# inc_test_filter  = !((calls[:year] .== 2012) .* (calls[:month] .<= 3))
# test_calls = calls[(1:nrow(calls))[inc_test_filter][1:ncalls],:];
test_calls = DataFrames.readtable("test_calls.csv")

problem = DispatchProblem(test_calls, hospitals, stations, p.coverage, x, turnaround=turnaround)
dispatch = ClosestDispatch(p, problem)
redeploy = AssignmentModel(p, x, hospitals, stations, lambda=Float64(lambda))

# id 145 dispatch to nbhd 88
srand(1234) # reset seed
@time df = simulate_events!(problem, dispatch, redeploy)
@show mean(df[:waittime]), maximum(df[:waittime])
@show mean(df[:waittime] + df[:responsetime])

# julia> mean(df[:responsetime])
# 8.34253333333334

test_calls = DataFrames.readtable("influx_calls.csv");

problem = DispatchProblem(test_calls, hospitals, stations, p.coverage, x, turnaround=turnaround);
dispatch = ClosestDispatch(p, problem);
redeploy = AssignmentModel(p, x, hospitals, stations, lambda=Float64(lambda));

# id 145 dispatch to nbhd 88
srand(1234); # reset seed
@time df = simulate_events!(problem, dispatch, redeploy);
@show mean(df[:waittime]), maximum(df[:waittime])
@show mean(df[:waittime] + df[:responsetime])


text:
"""
Module defining random play.

This module is the reference implementation for game decisions that should be
    overridden when improving strategy.
"""
module RandomPlay

import PlayingCards: Card

import Cribbage: Player
import Cribbage.Rules: validhands, validplays

export discardtocrib, pickcardtoplay

struct RandomPlayer <: Player
    name :: String
end

"""
    discardtocrib(player, cards)
Divides `cards` into `hand`, `crib`
"""
function discardtocrib(player::Player, cards::Tuple{Vararg{Card}})
    subsets = validhands(cards)
    hand = rand(subsets)
    crib = setdiff(cards, hand)
    @debug "keeping $(string(hand)), discarding $(string(crib))"
    return hand, crib
end

"""
    pickcardtoplay(player, prior, cards)
Returns either
- a card from `cards` to play
- `missing` if there are no cards can be played
"""
function pickcardtoplay(player::Player, prior::Tuple{Vararg{Card}}, cards::Tuple{Vararg{Card}})::Union{Missing,Card}
    possibleplays = validplays(prior, cards)
    if length(possibleplays) == 0
        return missing
    else
        return rand(possibleplays)
    end
end

end  # module RandomPlay

ext: jl


text:
#=
This shows basic data analysis using the Card data
on returns to education

See the github pages of the packages for documentation

=#

##
using CSV, DataFrames, DataFramesMeta, DelimitedFiles, StatPlots, Plots
pyplot(reuse=false, show=true)
card = CSV.read("../Data/card.csv") # how to read CSV data into DataFrame
size(card)
names(card)
head(card)
describe(card)

##
# select south and black (using DataFramesMeta)
sb = @where(card, :south .== 1, :black .== 1)
# make a scatter plot
@df sb scatter(:educ, :wage, leg=false)
# select north and white
nw = @where(card, :south .== 0, :black .== 0)
@df nw scatter!(:educ, :wage, leg=false)
##
# do the same thing using the all the data, with grouping
@df card scatter(:educ, :wage, group = (:black, :south), markersize=10)
##
# density plots by groups
@df card density(:wage, group=(:black), fill=true, fillalpha=0.3, legend=false)

##
# violin and box plots
@df card violin(:wage)
@df card boxplot!(:wage) # puts is beside
# violin plot of wage, by education
@df card violin(:educ, :wage)
# add a boxplot on top
@df card boxplot!(:educ, :wage)
##
# convert the dataframe to an ordinary array,
# as this is usually better for more advanced
# modeling
data = convert(Array{Float64}, card)
writedlm("card.ascii", data) # write to plain ASCII file


text:
flip( a::Dict ) = Dict( Iterators.flatten( [ v .=> k for (k,v) in a ] ) )

struct WorkerMetaData
    workertype::Type
    address::Union{Missing, String}
    CPU_speed::Int32
    RAM_available::UInt64
    #port::Int
end

"""
        worker_meta()

Returns a two dictionaries (local, and compute) containing metadata of all workers in scope.
This is for quick access to otherwise nested information and seperation of compute resources.

"""
function worker_meta()
    local_metadata, compute_metadata   = Dict(),  Dict()
    @sync for worker in Distributed.PGRP.workers
        @assert( !haskey(local_metadata, worker.id) )
        @assert( !haskey(compute_metadata, worker.id) )
        #Handle special case of master node
        if isa( worker, Distributed.LocalProcess )
            try
                ram = Sys.free_memory()
                cpuspd = Sys.cpu_info()[1].speed
                local_metadata[ worker.id ]      = WorkerMetaData( Distributed.LocalProcess,
                                                                worker.bind_addr, cpuspd, ram )
            catch err
                if isa(err, UndefRefError)
                    ram = Sys.free_memory()
                    cpuspd = Sys.cpu_info()[1].speed
                    local_metadata[ worker.id ]   = WorkerMetaData( Distributed.LocalProcess,
                                                                    missing, cpuspd, ram  )
                end
            end
        else
            cpuspd = fetch( @spawnat worker.id Sys.cpu_info()[1].speed )
            ram = fetch( @spawnat worker.id Sys.free_memory() )
            compute_metadata[ worker.id ] = WorkerMetaData( typeof( worker.manager ),
                                                            worker.config.host, cpuspd, ram )
        end
    end
    if length(compute_metadata) == 0
        compute_metadata = missing
    end
    return local_metadata, compute_metadata
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
"""
    CurtailableLoad

Curtailable load
"""
struct CurtailableLoad <: Resource
    index::Int           # Index of the resource
    num_timesteps::Int   # Number of time-steps in the battery's operation

    load::Vector{Float64}     # Load at each time-step, if device is on
    binflag::Bool  # Whether on-off or continuous curtailment

    function CurtailableLoad(;
        index::Integer=0,
        T::Integer=0,
        load::Vector{Float64}=Float64[],
        binflag::Bool=true
    )
        # Dimension checks
        T == size(load, 1) || throw(DimensionMismatch("Invalid load dimension"))
        return new(index, T, load, binflag)
    end

end

function add_resource_to_model!(
    m::MOI.ModelLike,
    h::House,
    l::CurtailableLoad,
    var2idx, con2idx
)
    T = l.num_timesteps

    # ==========================================
    #    I. Add local variables
    # ==========================================
    pnet = MOI.add_variables(m, T)  # Net load
    curt = MOI.add_variables(m, T)  # Curtailment indicator
    # Update list of variable indices
    for t in 1:T
        var2idx[(:curt, l.index, :pnet, t)] = pnet[t]
        var2idx[(:curt, l.index, :uind, t)] = curt[t]
    end

    for t in 1:T
        if l.binflag
            # Add binary constraint
            MOI.add_constraint(m,
                MOI.SingleVariable(curt[t]),
                MOI.ZeroOne()
            )
        else
            # Just add bounds
            MOI.add_constraint(m,
                MOI.SingleVariable(curt[t]),
                MOI.Interval(0.0, 1.0)
            )
        end
    end

    # ==========================================
    #    II. Update linking constraints
    # ==========================================
    for t in 1:T
        cidx = con2idx[(:house, h.index, :link, t)]
        MOI.modify(m, cidx, MOI.ScalarCoefficientChange(pnet[t], -1.0))
    end

    # ==========================================
    #    III. Add local constraints
    # ==========================================
    # Curtailement
    #   p[t] - u[t] * P[t] = 0
    for t in 1:T
        cidx = MOI.add_constraint(m,
            MOI.ScalarAffineFunction{Float64}(
                [
                    MOI.ScalarAffineTerm(1.0, pnet[t]),
                    MOI.ScalarAffineTerm(-l.load[t], curt[t])
                ],
                0.0
            ),
            MOI.EqualTo(0.0)
        )
        con2idx[(:curt, l.index, :curtail, t)] = cidx
    end

    # ==========================================
    #    IV. Add sub-resources to current model
    # ==========================================
    # (None here)

    return nothing
end


text:
###############################################################################
#
# Recursive display function
#
###############################################################################
"""
Display the parameters in a colored way.

$(METHODLIST)

"""
function pretty_display! end

"""
When a dictionary (a pair of key and item) is passed to the `pretty_display!` function, the dictionary will be displayed as `(LEADING_SPACE)key ⇨ item,`.
    However, if the dictionary item is an array of dictionaries, then the array will be displayed recursively.

    pretty_display!(dict::Pair, max_len::Int, spaces = "    ")

Display the key and item in a dictionary (Pair), given
- `dict` Dictionary to display
- `max_len` Maximum length of dictionary key (in an array)
- `spaces` Leading spaces before displaying the dictionary key

---
Examples
```julia
# display the pair in a pretty way (automatic alignments)
pretty_display!("a" => "b", 1);
pretty_display!("a" => "b", 2, "  ");
```
"""
pretty_display!(dict::Pair, max_len::Int, spaces = "    ") =
(
    # print leading spaces
    print(spaces);

    # print the key
    printstyled(string(dict[1]); color=:light_magenta);

    # print spaces after key and arrow
    print( repeat(" ", max_len-length(string(dict[1]))) * " ⇨ " );

    # if the value is a vetor of pairs, recursive display
    if typeof(dict[2]) <: Vector
        # display [ and line break
        print( "[\n" );
        pretty_display!(dict[2], spaces*repeat(" ",max_len+5));

        # display a ] and the next line
        print(spaces * repeat(" ",max_len+4) * "],\n" );

    # if the value is not array of pair, display it
    else
        printstyled(string(dict[2]); color=:cyan);
        print(",\n" );
    end;

    return nothing;
)

"""
When an array of dictionaries is given, the `pretty_display!` function computes the maximum length of the dictionary keys, and then display the
    dictionaries in a colored and pretty manner:

    pretty_display!(
                dicts::Union{Vector{Pair{String,String}},
                             Vector{Pair{String,Any}},
                             Vector{Pair{Any,String}},
                             Vector{Pair{Any,Any}}},
                spaces::String = "    ")

Display array of pairs (dictionary) in recursive manner, given
- `dicts` Parameters to display
- `spaces` Leading spaces

---
Examples
```julia
# display the vector of pairs recursively (automatic alignments)
_dicts = ["A" => "b", "d" => "A", "rr" => ["ra" => "rB", "rD" => "ra"]];
pretty_display!(_dicts);
pretty_display!(_dicts, "  ");
```
"""
pretty_display!(
            dicts::Union{Vector{Pair{String,String}},
                         Vector{Pair{String,Any}},
                         Vector{Pair{Any,String}},
                         Vector{Pair{Any,Any}}},
            spaces::String = "    ") =
(
    # determine the length of the keys
    _max_len = maximum(length.([string(p[1]) for p in dicts]));

    # display the elements
    pretty_display!.(dicts, _max_len, spaces);

    return nothing
)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule tectonic_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("tectonic")
JLLWrappers.@generate_main_file("tectonic", UUID("d7dd28d6-a5e6-559c-9131-7eb760cdacc5"))
end  # module tectonic_jll


text:
# Routines for taking steps using either Euler or SDC

using FastGaussQuadrature

function step_euler(dt, sol, solver, ubc1, ubc2)
    alpha = solver.alpha
    F1, F2 = ugradu(sol)
    RHS1 = alpha^2*(sol.U1 - dt*F1)
    RHS2 = alpha^2*(sol.U2 - dt*F2)
    newsol = solver.solve(RHS1, RHS2, ubc1, ubc2)
    return newsol
end

function step_bdf(dt, prevsol, sol, solver, ubc1, ubc2)
    alpha = solver.alpha
    F1, F2 = ugradu(sol)
    prevF1, prevF2 = ugradu(prevsol)
    RHS1 = 1/3*alpha^2*(4*sol.U1 - prevsol.U1 - 2*dt*(2*F1 - prevF1))
    RHS2 = 1/3*alpha^2*(4*sol.U2 - prevsol.U2 - 2*dt*(2*F2 - prevF2))
    newsol = solver.solve(RHS1, RHS2, ubc1, ubc2)
    return newsol
end

function sdc_substep(dt, sdc_order)
    return sdc_order > 1 ? dt / (sdc_order-1) : dt
end

function sdc_integration(M)
    # Setup substeps
    x = collect(linspace(-1, 1, M))
    # Vandermonde matrix (transposed)
    At = ones(M,M)
    for m=2:M
        @. At[m, :] = At[m-1, :]*x
    end
    # Use Gauss-Legendre for quadrature
    xleg, wleg = gausslegendre(M)
    # Compute substep integration weights
    I = zeros(M,M-1)
    for m=1:M-1
        a, b = x[m], x[m+1]
        wint = wleg*(b-a)/2
        xint = (xleg+1)*(b-a)/2 + a
        xint_i = ones(M)
        p = zeros(M)
        p[1] = sum(wint)
        for i=2:M
            xint_i .*= xint
            p[i] = sum(wint .* xint_i)
        end
        I[:,m] = At\p
    end
    return I/2
end

function step_sdc(dt, order, sol, solver, ubc1, ubc2)
    # TODO: Need to implement BC time dependence
    alpha = solver.alpha
    Re = alpha^2*dt
    P = order # order-1 can also work
    K = order
    if order>1
        W = sdc_integration(P)
    else
        # Fallback to Euler
        P = 2
        K = 1
    end
    # Set initial conditions for all passes
    solutions = Array{Solution}(P, K)
    for k=1:K
        solutions[1, k] = sol
    end
    # First pass (k=0)
    # Eulers steps
    for m=1:P-1
        init = solutions[m, 1]
        F1init, F2init = ugradu(init)
        RHS1 = alpha^2*init.U1 - Re*F1init
        RHS2 = alpha^2*init.U2 - Re*F2init
        solutions[m+1, 1] = solver.solve(RHS1, RHS2, ubc1, ubc2)
    end

    # Subsequent passes
    for k=1:K-1
        for m=1:P-1
            # Load historic stages that we need
            init = solutions[m,   k+1] # Substep initial condition
            pred = solutions[m+1, k  ] # Prediction
            hist = solutions[m  , k  ] # Prediction history
            # Post-process stages
            F1init, F2init = ugradu(init)
            #
            L1pred = alpha^2*pred.U1 - pred.RHS1
            L2pred = alpha^2*pred.U2 - pred.RHS2
            #
            F1hist, F2hist = ugradu(hist)
            # Integrate
            s = size(init.U1)
            I1, I2 = zeros(s), zeros(s)
            L = dt*(order-1)
            w = W[:, m]*L
            for n=1:P
                tmp = solutions[n, k]
                F1tmp, F2tmp = ugradu(tmp)
                L1tmp = alpha^2*tmp.U1 - tmp.RHS1
                L2tmp = alpha^2*tmp.U2 - tmp.RHS2
                I1 += w[n]*( L1tmp - Re*F1tmp )
                I2 += w[n]*( L2tmp - Re*F2tmp )
            end

            # f(t) = 1 + t
            # @show t1 = quadgk(f, dt*(m-1), dt*m)[1]
            # @show t2 = sum(w .* f(0:dt:dt*(P-1)))
            # @show t1-t2
            # @assert abs(t1-t2)/abs(t1) < 1e-14

            # Create new RHS
            RHS1 = alpha^2*init.U1 - Re*(F1init-F1hist) - L1pred + I1/dt
            RHS2 = alpha^2*init.U2 - Re*(F2init-F2hist) - L2pred + I2/dt
            #
            solutions[m+1, k+1] = solver.solve(RHS1, RHS2, ubc1, ubc2)
        end
    end

    if K>1
        # Output error est (UNTESTED)
        maxnorm(v) = maximum(abs.(v))
        l2norm(v) = sqrt(sum(vec(v.^2)))
        corr1 = solutions[P, K].U1 - solutions[P, K-1].U1
        corr2 = solutions[P, K].U2 - solutions[P, K-1].U2
        sol = solutions[P,K]
        info("U1 SDC err est: l2=", l2norm(corr1)/l2norm(sol.U1),", max=", maxnorm(corr1)/maxnorm(sol.U1))
        info("U2 SDC err est: l2=", l2norm(corr2)/l2norm(sol.U2),", max=", maxnorm(corr2)/maxnorm(sol.U2))
    end
    return solutions[P,K]

    # Original order 2 code
    # # First take Euler step
    # F1, F2 = ugradu(sol)
    # RHS1 = alpha^2*U1 - Re*F1
    # RHS2 = alpha^2*U2 - Re*F2
    # prov = solver.solve(RHS1, RHS2, ubc1, ubc2)
    # # Post process Euler step
    # F1prov, F2prov = ugradu(prov)
    # L1prov = alpha^2*prov.U1 - prov.RHS1
    # L2prov = alpha^2*prov.U2 - prov.RHS2
    # # Post process initial conds
    # L1 = alpha^2*sol.U1 - sol.RHS1
    # L2 = alpha^2*sol.U2 - sol.RHS2
    # # Integrate
    # I1 = (L1-Re*F1 + L1prov-Re*F1prov)/2
    # I2 = (L2-Re*F2 + L2prov-Re*F2prov)/2
    # # Take correction step
    # RHS1 = Re/dt*U1 - L1prov + Re*I1
    # RHS2 = Re/dt*U2 - L2prov + Re*I2
    # newsol = solver.solve(RHS1, RHS2, ubc1, ubc2)
    # return newsol
end


text:
import Base: +, -, ==, <, <=, convert, hash, promote_rule, isfinite

"""
Subtype of `Real` defining a lexicographically comparable pair of `Real`.
It is designed to be used by `Scheduler` where standard real numbers run to a
problem of undefined order of undefined order of removal from priority queue.

`PriorityTime` two fields `time` and `priority` may have different types,
but both have to be subtypes of `Real`.
`priority` should be used to determine order of execution of `Action`s that
have the same time. Two actions with identical `time` and `priority`
have undefined oreder of execution so this should be avoided.

`PriorityTime` type has defined lexicographic order and `+`, `-`.
It is immutable, has a custom `hash` function and conversions from `Real` types.
"""
struct PriorityTime{T1<:Real, T2<:Real}<:Real
    time::T1
    priority::T2
end

"""
    PriorityTime(time)

Construct `PriorityTime` with `priority` randomly generated using `rand()`.
"""
PriorityTime(time::T1) where T1<:Real = PriorityTime{T1,Float64}(time, rand())

# need to define conversion from PriorityTime to PriorityTime to avoid recursion
# as PriorityTime is a subtype of Real
function convert(T::Type{PriorityTime{T1,T2}},
                 x::PriorityTime{T3,T4}) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real}
    PriorityTime{T1,T2}(convert(T1, x.time), convert(T2, x.priority))
end

# Reals other than PriorityTime convert with priority equal to zero of their type
function convert(T::Type{PriorityTime{T1,T2}}, x::T3) where {T1<:Real, T2<:Real, T3<:Real}
    PriorityTime{T1,T2}(convert(T1, x), zero(T2))
end

# need to define promotion from PriorityTime to PriorityTime to avoid recursion
# as PriorityTime is a subtype of Real
function promote_rule(::Type{PriorityTime{T1,T2}},
                      ::Type{PriorityTime{T3,T4}}) where {T1<:Real, T2<:Real, T3<:Real, T4<:Real}
    PriorityTime{promote_type(T1,T3),promote_type(T2,T4)}
end

# Reals other than PriorityTime assume that they are converted to time field.
# Result of promotion is PriorityTime
function promote_rule(::Type{PriorityTime{T1,T2}},
                      ::Type{T3}) where {T1<:Real, T2<:Real, T3<:Real}
    PriorityTime{promote_type(T1,T3),T2}
end

# Need to also handle the case when PriorityTime is a second argument in promotion.
function promote_rule(::Type{T3},
                      ::Type{PriorityTime{T1,T2}}) where {T1<:Real, T2<:Real, T3<:Real}
    PriorityTime{promote_type(T1,T3),T2}
end

# We use tuple hashing algorithm.
function hash(x::PriorityTime{T1,T2}, h::UInt) where {T1<:Real, T2<:Real}
    hash((x.time, x.priority), h)
end

# in addition and subtraction both fields are updated
# so that the operations are commutative
function +(x::PriorityTime{T1,T2}, y::PriorityTime{T1,T2}) where {T1<:Real, T2<:Real}
    PriorityTime{T1,T2}(x.time+y.time, x.priority+y.priority)
end

function -(x::PriorityTime{T1,T2}, y::PriorityTime{T1,T2}) where {T1<:Real, T2<:Real}
    PriorityTime{T1,T2}(x.time-y.time, x.priority-y.priority)
end

# we define lexicographic order on PriorityTime
# where both objectives are maximized and time is more important
function ==(x::PriorityTime{T1,T2}, y::PriorityTime{T1,T2}) where {T1<:Real, T2<:Real}
    x.time==y.time && x.priority==y.priority
end

function <(x::PriorityTime{T1,T2}, y::PriorityTime{T1,T2}) where {T1<:Real, T2<:Real}
    x.time < y.time && return true
    x.time == y.time && x.priority < y.priority && return true
    return false
end

function <=(x::PriorityTime{T1,T2}, y::PriorityTime{T1,T2}) where {T1<:Real, T2<:Real}
    x.time < y.time && return true
    x.time == y.time && x.priority <= y.priority && return true
    return false
end


text:
#RSDeltaSigma: exampleHBF() function
#-------------------------------------------------------------------------------

function HBFex1()
	q1 = [
		 0 -4 -5;
		 1  1  1;
		 0 -4 -5;
		-1 -1 -1;
		-1 -3 -5;
		 1  1  1;
		-3 -5  0;
		-1 -1  0;
	]
	q2 = [
		-1 -3 -7;
		 1  1  1;
		-2 -5 -7;
		-1  1  1;
		-3 -8  0;
		 1 -1  0;
		-4 -6 -8;
		-1 -1 -1;
		-4 -8  0;
		 1 -1  0;
		-5 -7 -8;
		-1 -1 -1;
		-5  0  0;
		 1  0  0;
		-6 -7  0;
		-1 -1  0;
		-7 -8  0;
		-1 -1  0;
		-6  0  0;
		 1  0  0;
	]
	return (q1,q2)
end

"""`HBFex2()`

Coefficients from toolbox v2.0 documentation.
"""
function HBFex2()
	q1 = [
		 0 -4 -7;
		 1 -1  1;
		-1 -3 -6;
		-1 -1 -1;
		-2 -4 -7;
		 1 -1  1;
	]
	q2 = [
		-1 -3 -8;
		 1  1 -1;
		-2 -4 -9;
		-1  1 -1;
		-3 -5 -9;
		 1 -1  1;
		-4 -7 -8;
		-1  1  1;
		-5 -8 -11;
		 1 -1 -1;
		-6 -9 -11;
		-1  1 -1;
	]
	return (q1,q2)
end

"""`Saramaki88()`

Coefficients from "Efficient VLSI-Realizable Decimators for Sigma-Delta
Analog-to-Digital Converters," T. Saramaki and H. Tenhunen, ISCAS 1988,
pp 1525-2528.
"""
function Saramaki88()
	q1 = [
		 0 -2;
		 1 -1;
		-2 -9;
		-1  1;
	]
	q2 = [
		-1 -3  0;
		 1  1  0;
		-3 -5 -6;
		-1 -1 -1;
		-4 -6  0;
		 1  1  0;
	]
	return (q1,q2)
end

"""`Saramaki90()`

Coefficients from "Multiplier-Free Decimation Algorithms for Superresolution
Oversampled Converters," T. Saramaki, T. Karema, T. Ritoniemi, and H. Tenhunen
ISCAS 1990, vol. 4. pp 1525-2528.
"""
function Saramaki90()
	q1 = [
		 0 -4 -5;
		 1  1  1;
		 0 -4 -5;
		-1 -1 -1;
		-1 -3 -5;
		 1  1  1;
		-3 -5 -20;
		-1 -1  1;
	]
	q2 = [
		-1 -3 -7;
		 1  1  1;
		-2 -5 -7;
		-1  1  1;
		-3 -8  0;
		 1 -1  0;
		-4 -6 -8;
		-1 -1 -1;
		-4 -8  0;
		 1 -1  0;
		-5 -7 -8;
		-1 -1 -1;
		-5  0  0;
		 1  0  0;
		-6 -7  0;
		-1 -1  0;
		-6  0  0;
		 1  0  0;
		-7 -8  0;
		-1 -1  0;
		-6  0  0;
		 1  0  0;
	]
	return (q1,q2)
end

function exampleHBF(n::Int = 0)
	fnlist = [HBFex1, HBFex2, Saramaki88, Saramaki90]
	if !in(n, keys(fnlist))
		throw("exampleHBF: No code for example $n.")
	end
	fn = fnlist[n]
	(q1, q2) = fn()
	F1 = convertForm(bunquantize(q1),q1)
	F2 = convertForm(bunquantize(q2),q2)
	return (F1, F2)
end

#Last line


text:
add(solver_name::AbstractString, v::Bool=true) = error("""
                                                        `POMDPs.add` has been removed.
                                                        To add POMDPs packages, use `POMDPs.add_registry()` with standard `Pkg.add`
                                                       """)

add_all(;native_only=false, v::Bool=true) = error("""
                                                    `POMDPs.add_all` has been removed.
                                                    To add POMDPs packages, use `POMDPs.add_registry()` with standard `Pkg.add`
                                                  """)
remove_all() = error("""
                     `POMDPs.remove_all` has been removed.
                     To remove POMDPs packages, use `POMDPs.add_registry()` with standard `Pkg.rm`
                     """)

build() = error("""
                `POMDPs.build` has been removed.
                To build POMDPs packages, use `POMDPs.add_registry()` with standard `Pkg.build`
                """)

test_all(v::Bool=false) = error("""
                                `POMDPs.test_all` has been removed.
                                To test POMDPs packages, use `POMDPs.add_registry()` with standard `Pkg.test`
                                """)

"""
    available()

Prints all the available packages in the JuliaPOMDP registry
"""
function available()
    reg_dict = read_registry(joinpath(Pkg.depots1(), "registries", "JuliaPOMDP", "Registry.toml"))
    for (uuid, pkginfo) in reg_dict["packages"]
        println(pkginfo["name"])
    end
end

"""
    add_registry()

Adds the JuliaPOMDP registry
"""
function add_registry(;url=POMDP_REGISTRY)
    depot = Pkg.depots1()
    # clone to temp dir first
    tmp = mktempdir()
    Base.shred!(LibGit2.CachedCredentials()) do creds
        LibGit2.with(Pkg.GitTools.clone(url, tmp; header = "registry from $(repr(url))", credentials = creds)) do repo
        end
    end
    # verify that the clone looks like a registry
    if !isfile(joinpath(tmp, "Registry.toml"))
        Pkg.Types.pkgerror("no `Registry.toml` file in cloned registry")
    end

    registry = read_registry(joinpath(tmp, "Registry.toml"))
    verify_registry(registry)

    # copy to depot
    regpath = joinpath(depot, "registries", registry["name"])
    ispath(dirname(regpath)) || mkpath(dirname(regpath))
    if Pkg.Types.isdir_windows_workaround(regpath)
        existing_registry = read_registry(joinpath(regpath, "Registry.toml"))
        @assert registry["uuid"] == existing_registry["uuid"]
        @info("registry `$(registry["name"])` already exists in `$(Base.contractuser(dirname(regpath)))`")
    else
        cp(tmp, regpath)
        Pkg.Types.printpkgstyle(stdout, :Added, "registry `$(registry["name"])` to `$(Base.contractuser(dirname(regpath)))`")
    end

end

function read_registry(regfile)
    registry = Pkg.TOML.parsefile(regfile)
    return registry
end

const REQUIRED_REGISTRY_ENTRIES = ("name", "uuid", "repo", "packages")

function verify_registry(registry::Dict{String, Any})
    for key in REQUIRED_REGISTRY_ENTRIES
        haskey(registry, key) || Pkg.Types.pkgerror("no `$key` entry in `Registry.toml`.")
    end
end


text:
import LightGraphs; lg = LightGraphs

export LabeledGraph, labels, graph_to_gr, graph_to_cnf
export get_vertex, vertices, nv, add_vertex!, rem_vertex!
export edges, ne, add_edge!, has_edge, rem_edge!
export degree, all_neighbors, eliminate!, cliqueness
export line_graph, tree_from_tree_decompostion, chordal_graph

# **************************************************************************************** #
#                                  Labeled Graph Struct
# **************************************************************************************** #

# TODO: It may be better to implement a subtype of AbstractGraph where the fadjlist
# is replaced with a dictionary of adjacency lists which can be indexed by symbols.
# This would avoid any iteration over vertices required for finding a vertex with a
# particular label.

"""
Struct to represent a labeled graph. Unique symbols are created for each vertex if none are
provided and remain unaltered during the lifetime of the graph.

# Example
```julia-repl
julia> g = LabeledGraph()
LabeledGraph({0, 0} undirected simple Int64 graph, Symbol[])

julia> add_vertex!(g, :a_vertex_label)
1-element Array{Symbol,1}:
 :a_vertex_label

julia> g.labels[1]
:a_vertex_label
```
"""
struct LabeledGraph
    graph::lg.AbstractGraph
    labels::Vector{Symbol}

    LabeledGraph() = new(lg.SimpleGraph(), Symbol[])
    LabeledGraph(G::lg.AbstractGraph) = new(deepcopy(G), Symbol.([1:lg.nv(G)...]))
    LabeledGraph(G::lg.AbstractGraph, labels::Array{Symbol, 1}) = new(deepcopy(G), labels)
    LabeledGraph(N::Int) = LabeledGraph(lg.SimpleGraph(N))
    LabeledGraph(labels::Array{Symbol, 1}) = new(lg.SimpleGraph(length(labels)), labels)
end

"""
    graph_to_gr(G::LabeledGraph, filename::String)

Write the provided graph to a file in gr format.
"""
graph_to_gr(G::LabeledGraph, filename::String) = graph_to_gr(G.graph, filename::String)

function graph_to_gr(G::lg.AbstractGraph, filename::String)
    open(filename, "w") do io
        write(io, "p tw $(lg.nv(G)) $(lg.ne(G))\n")
        for e in lg.edges(G)
            write(io, "$(e.src) $(e.dst)\n")
        end
    end
end

"""
    graph_from_gr(filename::String)

Read a graph from the provided gr file.
"""
function graph_from_gr(filename::String)
    lines = readlines(filename)

    # Create a Graph with the correct number of vertices.
    num_vertices, num_edges = parse.(Int, split(lines[1], ' ')[3:end])
    G = lg.SimpleGraph(num_vertices)

    # Add an edge to the graph for every other line in the file.
    for line in lines[2:end]
        src, dst = parse.(Int, split(line, ' '))
        lg.add_edge!(G, src, dst)
    end
    LabeledGraph(G)
end

"""
    graph_to_cnf(G::LabeledGraph, filename::String)

Write the provided graph to a file in cnf format.
"""
graph_to_cnf(G::LabeledGraph, filename::String) = graph_to_cnf(G.graph, filename::String)

function graph_to_cnf(G::lg.AbstractGraph, filename::String)
    open(filename, "w") do io
        write(io, "p cnf $(lg.nv(G)) $(lg.ne(G))\n")
        for e in lg.edges(G)
            write(io, "$(e.src) $(e.dst) 0\n")
        end
    end
end

"""
    graph_from_cnf(filename::String)

Read a graph from the provided cnf file.
"""
function graph_from_cnf(filename::String)
    lines = readlines(filename)

    # Create a Graph with the correct number of vertices.
    num_vertices, num_edges = parse.(Int, split(lines[1], ' ')[3:end])
    G = lg.SimpleGraph(num_vertices)

    # Add an edge to the graph for every other line in the file.
    for line in lines[2:end]
        src, dst = parse.(Int, split(line, ' ')[1:2])
        lg.add_edge!(G, src, dst)
    end
    LabeledGraph(G)
end

# **************************************************************************************** #
#                                Labeled Graph Interface
# **************************************************************************************** #

"""
    labels(G::LabeledGraph)

Return the labels contained in a LabeledGraph.
"""
labels(G::LabeledGraph) = G.labels

"""
    get_vertex(G::LabledGraph, v_label)

Return the first vertex whose label matches the argument `v_label`. If an array of labels is
provided, an array of the corresponding vertices is return.
"""
function get_vertex(G::LabeledGraph, v_label::Symbol)
    for v in vertices(G)
        if G.labels[v] == v_label
            return v
        end
    end
    nothing
end

function get_vertex(G::LabeledGraph, v_labels::Array{Symbol, 1})
    verts = Array{Int, 1}(undef, length(v_labels))
    for i = 1:length(verts)
        verts[i] = get_vertex(G, v_labels[i])
    end
    verts
end

"""
    vertices(G::LabeledGraph)

Return the vertices of a labeled graph.
"""
function vertices(G::LabeledGraph)
    lg.vertices(G.graph)
end

"""
    nv(G::LabeledGraph)

Return the number of vertices in `G`.
"""
function nv(G::LabeledGraph)
    lg.nv(G.graph)
end

"""
    add_vertex!(G::LabeledGraph, label::Symbol)

Add a new vertex to `G` and assign the given label to it.
"""
function add_vertex!(G::LabeledGraph, label::Symbol)
    lg.add_vertex!(G.graph)
    append!(G.labels, Symbol[label])
end

"""
    rem_vertex!(G::LabeledGraph, v)

Delete the vertex with index or label `v` from `G`.
"""
function rem_vertex!(G::LabeledGraph, v::Int)
    lg.rem_vertex!(G.graph, v)
    G.labels[v] = G.labels[end]
    deleteat!(G.labels, length(G.labels))
end

function rem_vertex!(G::LabeledGraph, v_label::Symbol)
    v = get_vertex(G, v_label)
    if !(v === nothing)
        return rem_vertex!(G, v)
    end
    false
end

"""
    edges(G::LabeledGraph)

Return an iterator of the edges of `G`.
"""
function edges(G::LabeledGraph)
    lg.edges(G.graph)
end

"""
    ne(G::LabeledGraph)

Return the number of edges in `G`.
"""
function ne(G::LabeledGraph)
    lg.ne(G.graph)
end

"""
    add_edge!(G::LabeledGraph, u::Int, v::Int)

Add an edge to `G` connecting vertices `u` and `v`.
"""
function add_edge!(G::LabeledGraph, u::Int, v::Int)
    lg.add_edge!(G.graph, u, v)
end

function add_edge!(G::LabeledGraph, u_label::Symbol, v_label::Symbol)
    u = get_vertex(G, u_label)
    v = get_vertex(G, v_label)
    lg.add_edge!(G.graph, u, v)
end

"""
    has_edge(G::LabeledGraph, u::Int, v::Int)

Return true if `G` has an edge connecting vertices `u` and `v`. Return false otherwise.
"""
function has_edge(G::LabeledGraph, u::Int, v::Int)
    lg.has_edge(G.graph, u, v)
end

"""
    rem_edge!(G::LabeledGraph, u::Int, v::Int)

Remove the edge connecting vertices `u` and `v` if it exists.
"""
function rem_edge!(G::LabeledGraph, u::Int, v::Int)
    lg.rem_edge!(G.graph, u, v)
end

"""
    degree(G::LabeledGraph[, v])

Return an array containing the degree of each vertex of `G`. If `v` is specified, only
return degrees for vertices in `v`.
"""
function degree(G::LabeledGraph)
    lg.degree(G.graph)
end

function degree(G::LabeledGraph, v::Union{Int, Array{Int, 1}})
    lg.degree(G.graph, v)
end

function degree(G::LabeledGraph, v_label::Union{Symbol, Array{Symbol, 1}})
    v = get_vertex(G, v_label)
    (v===nothing) ? false : lg.degree(G.graph, v)
end

"""
    all_neighbors(G::LabeledGraph, v::Int)

Return an array of all neighbors of `v` in `G`.
"""
function all_neighbors(G::LabeledGraph, v::Int)
    lg.all_neighbors(G.graph, v)
end

function all_neighbors(G::LabeledGraph, v_label::Symbol)
    v = get_vertex(G, v_label)
    lg.all_neighbors(G.graph, v)
end

"""
    eliminate!(G::LabledGraph, v)

Connect all the neighbors of `v` together before removing `v` from `G`.
"""
function eliminate!(G::LabeledGraph, v_label::Symbol)
    v = get_vertex(G, v_label)
    eliminate!(G, v)
end

"""
    eliminate!(G::AbstractGraph, v::Integer)

Connect all the neighbours of v together before removing v from G.
"""
function eliminate!(G::LabeledGraph, v::Integer)
    # Loop over all pairs of neighbours (vi, ui) of v and try connect them with an edge
    # before removing v from G.
    Nᵥ = all_neighbors(G, v)::Array{Int64, 1}
    for i in 1:length(Nᵥ)
        vi = Nᵥ[i]
        for j in i+1:length(Nᵥ)
            uj = Nᵥ[j]
            add_edge!(G, vi, uj)
        end
    end
    rem_vertex!(G, v)
end

"""
    eliminate!(G::AbstractGraph, v::Integer, c_map::Dict{Int, Int})

Connect all the neighbours of v together before removing v from G.

The dictionary `c_map` mapping vertices of 'G' to their cliqueness is updated inplace
accordingly.
"""
function eliminate!(G::LabeledGraph, v::Integer, c_map::Dict{Int, Int})
    Nᵥ = all_neighbors(G, v)::Array{Int64, 1}
    # Loop over all pairs of vertices in the neighbourhood of 'v'.
    for i in 1:length(Nᵥ)
        vi = Nᵥ[i]
        for j in i+1:length(Nᵥ)
            ui = Nᵥ[j]

            # Try add an edge connecting vi and ui. If successful, update `c_map`.
            edge_added = add_edge!(G, vi, ui)::Bool
            if edge_added
                # Common neighbours of vi and ui have one less edge to add
                # when being eliminated after vi and ui are connected.
                Nvi = all_neighbors(G, vi)::Array{Int64, 1}
                Nui = all_neighbors(G, ui)::Array{Int64, 1}
                for n in intersect(Nvi, Nui)
                    c_map[n] -= 1
                end

                # ui and vi are now neighbours so their cliqueness may increase.
                for n in Nvi
                    if !(n == ui) && !(has_edge(G, n, ui)::Bool)
                        c_map[vi] += 1
                    end
                end
                for n in Nui
                    if !(n == vi) && !(has_edge(G, n, vi)::Bool)
                        c_map[ui] += 1
                    end
                end
            end
        end
    end

    # Removing v from G means it's also removed from its neighbour's neighbourhood, so their
    # cliqueness may be reduced.
    for n in Nᵥ
        Nₙ = all_neighbors(G, n)::Array{Int64, 1}
        for u in Nₙ
            if !(u == v)
                edge_exists = has_edge(G, v, u)::Bool
                if !edge_exists
                    c_map[n] -= 1
                end
            end
        end
    end
    # Remove v from G and the corresponding entry from c_map.
    N = nv(G)::Int
    c_map[v] = c_map[N]
    delete!(c_map, N)
    rem_vertex!(G, v)
end

"""
    cliqueness(G::LabeledGraph, v::Symbol)

Return the number of edges that need to be added to `G` in order to make the neighborhood of
vertex labeled by the symbol `v` a clique.
"""
cliqueness(G::LabeledGraph, v::Symbol)::Int = cliqueness(G, get_vertex(G, v))

"""
    cliqueness(G::LabeledGraph, v::Integer)

Return the number of edges that need to be added to `G` in order to make the neighborhood of
vertex `v` a clique.
"""
cliqueness(G::LabeledGraph, v::Integer)::Int = cliqueness(G.graph, v)

function cliqueness(G::lg.AbstractGraph, v::Integer)::Int
    neighborhood = lg.all_neighbors(G, v)::Array{Int64, 1}
    count = 0
    for i in 1:length(neighborhood)
        for j in i+1:length(neighborhood)
            vi = neighborhood[i]
            ui = neighborhood[j]
            if !lg.has_edge(G, ui, vi)::Bool
                count += 1
            end
        end
    end
    count
end

# **************************************************************************************** #
#                          Functions for creating various graphs
# **************************************************************************************** #

"""
    line_graph(G::LabeledGraph)

Return a LabeledGraph representing the line graph of the 'G'.

The label for each each vertex of the line graph is created by concatenating the labels of
the corresponding vertices in the LabeledGraph 'G'.
"""
function line_graph(G::LabeledGraph)
    G_edges = collect(edges(G))
    vertex_labels = [combine_labels(G.labels[e.src], G.labels[e.dst]) for e in G_edges]
    line_graph(G.graph; vertex_labels=vertex_labels)
end

"""
    line_graph(G::AbstractGraph;
               vertex_labels::Array{Symbol, 1}=Symbol[])

Return a LabeledGraph representing the line graph of the 'G'.

The label for each each vertex of the line graph is created by
concatenating the labels of the corresponding vertices in 'G'.

The symbols in the array `vertex_labels` are used as labels for the vertices of the returned
line graph. If `vertex_labels` is empty then labels are created by combining the indices of
the corresponding vertices in 'G'.
"""
function line_graph(G::lg.AbstractGraph;
                    vertex_labels::Array{Symbol, 1}=Symbol[])
    # Create a labeled graph LG whose vertices corresponding to the edges of G.
    G_edges = collect(lg.edges(G))
    if isempty(vertex_labels)
        vertex_labels = [combine_labels(e.src, e.dst) for e in G_edges]
    end
    LG = LabeledGraph(lg.SimpleGraph(length(G_edges)), vertex_labels)

    # Connect any two vertices of LG whose corresponding edges in G share a vertex in G.
    for i in 1:length(G_edges)-1
        for j in i+1:length(G_edges)
            u = G_edges[i]
            v = G_edges[j]
            if (u.src == v.src) || (u.src == v.dst) || (u.dst == v.src) || (u.dst == v.dst)
                add_edge!(LG, i, j)
            end
        end
    end
    LG
end

"""
    combine_labels(label_A, label_B)

Concatenate 'label_A' and 'label_B' in lexicographical order.
"""
function combine_labels(label_A, label_B)
    label_A < label_B ? Symbol(label_A, :_, label_B) : Symbol(label_B, :_, label_A)
end

"""
    tree_from_tree_decompostion(td::Dict{Symbol, Any})

Returns a LabeledGraph representing the tree described by the tree decomposition in `td`.
"""
function tree_from_tree_decompostion(td::Dict{Symbol, Any})
    tree = LabeledGraph([Symbol("b_$i") for i in 1:td[:num_bags]])
    for (u, v) in td[:edges]
        add_edge!(tree, u, v)
    end
    tree
end

"""
    chordal_graph(G::LabeledGraph, π̄::Array{Symbol, 1})

Return a chordal graph built from 'G' using the elimination order 'π̄'.

The returned graph is created from 'G' by iterating over the vertices of 'G', according to
the order 'π̄', and for each vertex, connecting all the neighbors that appear later in the
order.
"""
function chordal_graph(G::LabeledGraph, π̄::Array{Symbol, 1})
    if !(Set(π̄) == Set(G.labels)) || !(length(π̄) == nv(G))
        error("π̄ must be an elimination order for G")
    end

    # Use G as a starting point for H and create dictionary to keep track of which vertices
    # in the elimination order haven't been considered yet.
    H = deepcopy(G)
    higher_order = Dict{Symbol, Bool}(G.labels .=> true)

    # for each vertex v in the elimination order, connect the neighbors of v that appear
    # after v in the elimination order.
    for v in π̄
        higher_order[v] = false
        neighbors = all_neighbors(H, v)
        for i = 1:length(neighbors)-1
            for j = i+1:length(neighbors)
                u = H.labels[neighbors[i]]
                w = H.labels[neighbors[j]]
                if higher_order[u] && higher_order[w]
                    add_edge!(H, u, w)
                end
            end
        end
    end
    H
end

chordal_graph(G::LabeledGraph, π̄::Array{<:Integer, 1}) = chordal_graph(G, G.labels[π̄])


text:
"""
    parseRawData(file::String)

Takes as input a csv file whose lines are of the form
```
    p | q | temperature | optimizer | minimum configuration | starting configuration | timestamp  | timespan of optimization | status
```
and parses it into a DataFrame.
"""
function parseRawData(file::String)
    raw = readdlm(file,'|')
    raw[:,4] = map( Symbol , raw[:,4])
    raw[:,5] = map( u -> parse.(Float64,split(u[2:end-1],", ")) , raw[:,5])
    raw[:,6] = map( u -> u=="Any[]" ? Vector{Float64}() : parse.(Float64,split(u[2:end-1],", ")) , raw[:,6])
    raw[:,7] = map( string , raw[:,7])
    raw[:,9] = map( string |> u->u=="FEASIBLE_POINT" ? "good" : u , raw[:,9])
    tup = Vector{Tuple}()
    for row in eachrow(raw)
        if row[9] == "FEASIBLE_POINT"
            row[9] = "good"
        end
        push!(tup,Tuple(row))
    end
    data = DataFrame(tup)
    names = [:p, :q, :temperature, :optimizer, :minimum, :starting, :id, :time, :status]
    names!(data,names)
    data
end

"""
    preprocessRawDataframe(db::DataFrame)

Takes as input a raw-data DataFrame generated by ```parseRawData``` and returns a new DataFrame ovtained by:
1)  selecting only gradient-based simulations;
2)  selecting only the successful simulations based on the exit status;
3)  selecting only zero-temperature simulations;
4)  adding a column with the total number of valid simulations;
5)  adding a column with ```minimum``` rounded to the first decimal digit, and flipped such that the first non-null magnetization is positive;
6)  adding a column with the value of the energy of ```minimum``` rounded to the second decimal digit;
7)  adding a column with the value of the entropy of ```minimum``` rounded to the second decimal digit.
"""
function preprocessRawDataframe(db::DataFrame)
    # select gradient
    db = db[ db.optimizer .|> u->occursin("gradient",u |> string ) , :]
    # select good simulations
    db = db[db.status .== "good", :]
    # select zero temperature
    db = db[db.temperature .== 0., :]
    # totalTrials
    tot = nrow(db)
    db.totalTrials = ones(tot)*tot
    # add info on ground state configuration (rounded to 1 decimal, I'm interested in 0,+-1 basically) normalized to factor out the Z2 symmetry
    db.cleanMinimum = map( u -> roundVector(u[:minimum],digits=1) |> flipZ2 , db |> eachrow)
    # and info on free energy of the exact minimum found, rounded to the second decimal
    db.energy = map( u -> roundVector(modelFreeEnergy(u[:minimum],u[:p],u[:q],u[:temperature]),digits=2) , db |> eachrow)
    # and info on entropy of the clean minimum found, rounded to the second decimal
    db.entropy = map( u -> roundVector(modelEntropy(u[:minimum],u[:q]),digits=2) , db |> eachrow)

    return db
end

"""
    analyzeGS(db::DataFrame, file::String)

Takes a preprocessed DataFrame and extracts informations about the ground state (GS).
It saves the GS energy and the list of distinct GSs into a Dictionary with fields ```energy, states```.
It exports the results using the package ```JLD``` into the file specified by ```file```; ```file``` must have extension ```.jld```.

"""
function analyzeGS(db::DataFrame, file::String)

    groundStateEnergy = min(db.energy...)
    groundStates = db[db.energy .== groundStateEnergy,:].cleanMinimum |> unique

    save(
         file,
         Dict(
            "energy"     => groundStateEnergy,
            "states"     => groundStates
           )
        )
end

"""
    analyzeMinima(db::DataFrame, file::String)

Takes a preprocessed DataFrame and extracts informations about the local minima of the energy landscape.
It saves  the list of distinct minima, the linear size of their basins of attraction and an estimate of the volume of the basins of attraction into a Dictionary with fields ```states, basinsLinear, basinsVolume```.
It exports the results using the package ```JLD``` into the file specified by ```file```; ```file``` must have extension ```.jld```.

The estimate of the volume of a basin of attraction is given by the number of descents that ends in the corresponding minimum.

The estimate of the linear size of a basin of attraction is given by the distance between all initial conditions and the minimum currently under study.
This given, one can compute the maximum distance for example, or any other useful statistic.
"""
function analyzeMinima(db::DataFrame, file::String)

    # add useful columns to db
    db.distance = map( u -> norm(u[1]-u[2]), zip(db.cleanMinimum, db.starting))

    # states
    tmp = db.cleanMinimum |> countmap
    states = tmp |> keys |> collect

    l = length(states)
    basinNumber = Vector{Int64}(undef,l)
    basinDistances = Vector{Vector{Float64}}(undef,l)
    for (i,s) in enumerate(states)
        basinNumber[i] = tmp[s]
        basinDistances[i] = db.distance[db.cleanMinimum .== [s]]
    end
    save(
         file,
         Dict(
            "states"       => states,
            "basinsLinear" => basinDistances,
            "basinsVolume" => basinNumber
           )
        )
end

"""
    checkMinimaFound(db::DataFrame, file::String)

Takes a preprocessed DataFrame and returns the curve of the number of distinct minima found VS the number of descent runs performed as a ```Matrix``` of size ```numberRuns x 2```, and saves it into a dictionary ith fields ```saturationCurve```.
It exports the results using the package ```JLD``` into the file specified by ```file```; ```file``` must have extension ```.jld```.

Every time it is called, it recomputed this saturation curve using a different random permuation of the rows of ```db```.

"""
function checkMinimaFound(db::DataFrame, file::String)
    # saturation of metastable states vs total trials
    vec = db.cleanMinimum
    distinctMin = Set{Vector{Float64}}()
    history = Matrix{Float64}(undef,length(vec),2)
    perm = randperm(length(vec))
    for i in 1:length(vec)
        if !(vec[perm[i]] in distinctMin)
            push!(distinctMin, vec[perm[i]])
        end
        history[i,:] = [ i,length(distinctMin)]
    end
    save(
         file,
         Dict(
            "saturationCurve"       => history
           )
        )
end

"""
    whichQP(directory::String)

Assuming that ```directory``` contains only the raw data simulated, with one file per pair (p,q) named, for example, q17_p4.csv, it returns a list of tuples with the simulated values of (q,p).

"""
function whichQP(directory::String)
    res = []
    for file in readdir(directory)
        push!(res, Tuple(split(file[2:end-4],"_p") .|> u->parse(Int64,u) ))
    end
    return res
end

"""
    getTotalRuns(file::String)

Counts the numer of lines of file.

NOTICE: one could also use the built-in function ```countlines(file)```. It seems that on small files ```countlines``` is faster than the shell command ```wc```, while on a file of ~300000 lines ```wc``` is faster.
In all cases, the difference is not very important, and ```countlines``` allocates lesser memory than ```wc```.
"""
function getTotalRuns(file::String)
    read(pipeline(`wc -l $file`,`cut -d ' ' -f 1 `), String) |> u -> parse(Int,u)
end


text:
"""
$(TYPEDEF)

Abstract supertype used for information storage object the directed acyclic graph.
"""
abstract type AbstractCache end

"""
$(TYPEDSIGNATURES)

Function used to initialize the storage cache `d::AbstractCache`
for a given type of directed acycle graph `g::AbstractDirectedGraph`.
"""
function initialize!(::AbstractCache, ::AbstractDirectedGraph)
    error("Initialization method not defined for this combination
           of AbstractCache and AbstractDirectedGraph.")
end

"""
$(TYPEDEF)

Abstract supertype used for attributes stored in a cache.
"""
abstract type AbstractCacheAttribute end

"""
$(TYPEDSIGNATURES)

Initializes information in cache `c` for each node in `g` that may be used in a forward-pass of attribute `t`.
"""
f_init!(t::AbstractCacheAttribute, g::AbstractDirectedGraph, c::AbstractCache)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a variable `v` at index `k` in
cache `c` associated with graph `g` using information at index `k`.
"""
fprop!(t::AbstractCacheAttribute, v::Variable, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a subexpression `v` at index `k` in
cache `c` associated with graph `g` using information taken from the children of `k`.
"""
fprop!(t::AbstractCacheAttribute, v::Subexpression, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a expression `v` at index `k` in
cache `c` associated with graph `g` using information taken from the children of `k`.
"""
fprop!(t::AbstractCacheAttribute, v::Expression, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a parameter `v` at index `k` in
cache `c` associated with graph `g` using information at index `k`.
"""
fprop!(t::AbstractCacheAttribute, v::Parameter, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a constant `v` at index `k` in
cache `c` associated with graph `g` using information at index `k`.
"""
fprop!(t::AbstractCacheAttribute, v::Constant, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Initializes information in cache `c` for each node in `g` that may be used in a reverse-pass of attribute `t`.
"""
r_init!(t::AbstractCacheAttribute, g::AbstractDirectedGraph, c::AbstractCache)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a variable `v` at index `k` in  cache `c`
associated with graph `g` using information at index `k` taken from the parents of `k`.
"""
rprop!(t::AbstractCacheAttribute, v::Variable, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a subexpressions `v` at index `k` in  cache `c`
associated with graph `g` using information at index `k` taken from the parents of `k`.
"""
rprop!(t::AbstractCacheAttribute, v::Subexpression, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a expressions `v` at index `k` in  cache `c`
associated with graph `g` using information at index `k` taken from the parents of `k`.
"""
rprop!(t::AbstractCacheAttribute, v::Expression, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a parameters `v` at index `k` in  cache `c`
associated with graph `g` using information at index `k` taken from the parents of `k`.
"""
rprop!(t::AbstractCacheAttribute, v::Parameter, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

"""
$(TYPEDSIGNATURES)

Populates information associated with attribute `t` for a constant `v` at index `k` in  cache `c`
associated with graph `g` using information at index `k` taken from the parents of `k`.
"""
rprop!(t::AbstractCacheAttribute, v::Constant, g::AbstractDirectedGraph, c::AbstractCache, k::Int)

Base.@kwdef mutable struct VariableValues{T<:Real}
    x0::Vector{T}                          = T[]
    x::Vector{T}                           = T[]
    lower_variable_bounds::Vector{T}       = T[]
    upper_variable_bounds::Vector{T}       = T[]
    node_to_variable_map::Vector{Int}      = Int[]
    variable_to_node_map::Vector{Int}      = Int[]
    variable_types::Vector{VariableType}   = VariableType[]
end

@inline val(b::VariableValues{T}, i::Int) where T = @inbounds b.x[i]
@inline lbd(b::VariableValues{T}, i::Int) where T = @inbounds b.lower_variable_bounds[i]
@inline ubd(b::VariableValues{T}, i::Int) where T = @inbounds b.upper_variable_bounds[i]
val(b::VariableValues{T}) where T = b.x
lbd(b::VariableValues{T}) where T = b.lower_variable_bounds
ubd(b::VariableValues{T}) where T = b.upper_variable_bounds

function _get_x!(::Type{BranchVar}, out::Vector{T}, v::VariableValues{T}) where T<:Real
    @inbounds for i = 1:length(v.node_to_variable_map)
        out[i] = v.x[v.node_to_variable_map[i]]
    end
    return nothing
end

function _initialize_or_copy!(y::VariableValues{T}, x::VariableValues{T}, s::Symbol) where {T<:Real}
    isempty(getfield(y, s)) ? setfield!(y, s, copy(getfield(x, s))) : copy!(getfield(y, s), getfield(x, s))
end
function update_box_and_pnt!(y::VariableValues{T}, x::VariableValues{T}, update_box::Bool) where {T<:Real}
    if update_box
        _initialize_or_copy!(y, x, :x0)
        _initialize_or_copy!(y, x, :lower_variable_bounds)
        _initialize_or_copy!(y, x, :upper_variable_bounds)
    end
    _initialize_or_copy!(y, x, :x)
end

"""
    DirectedTree

A tree graph with a single sink node.
"""
Base.@kwdef mutable struct DirectedTree <: AbstractDirectedAcyclicGraph
    "List of nodes"
    nodes::Vector{Node}                         = Node[]
    "List of index of variables in this tree"
    variables::Dict{Int,Int}                    = Dict{Int,Int}()
    "Information on all variables..."
    v::VariableValues{Float64}                  = VariableValues{Float64}()
    "List of constant values"
    constant_values::Vector{Float64}            = Float64[]
    "List of constant values"
    parameter_values::Vector{Float64}           = Float64[]
    "Number of nodes"
    node_count::Int                             = 0
    "Number of variables"
    variable_count::Int                         = 0
    "Number of constants"
    constant_count::Int                         = 0
    sink_bnd::Interval{Float64}                 = Interval{Float64}(-Inf,Inf)
    ""
    sparsity::Vector{Int}                       = Int[]
    ""
    rev_sparsity::Dict{Int,Int}                 = Dict{Int,Int}()
    dependent_variable_count::Int               = 0
    dep_subexpr_count::Int                      = 0
    dependent_subexpressions::Vector{Int}       = Int[]
    dependent_subexpression_dict::Dict{Int,Int} = Dict{Int,Int}()
    linearity::Linearity                        = LIN_CONSTANT
    user_operators::OperatorRegistry            = OperatorRegistry()
    #children::SpraseMatrixCSC{Bool,Int}        = spzeros(Bool,1,1)
    #parents::SparseMatrixCSC{Bool,Int}         = spzeros(Bool,1,1)
end
#DirectedTree(n::Int) = DirectedTree(children = spzeros(Bool,n,n), parents=spzeros(Bool,n,n))
const DAT = DirectedTree

# node property access functions that can be defined at abstract type
node(g::DAT, i)            = g.nodes[i]
nodes(g::DAT)              = g.nodes

variable(g::DAT, i)        = g.variables[i]
variables(g::DAT)          = g.variables

constant_value(g::DAT, i)  = g.constant_values[i]
constant_values(g::DAT)    = g.constant_values
parameter_value(g::DAT, i) = g.parameter_values[i]
parameter_values(g::DAT)   = g.parameter_values

node_class(g::DAT, i)      = node_class(node(g, i))
ex_type(g::DAT, i)         = ex_type(node(g, i))
first_index(g::DAT, i)     = first_index(node(g, i))
secondary_index(g::DAT, i) = secondary_index(node(g, i))
arity(g::DAT, i)           = arity(node(g, i))
children(g::DAT, i)        = children(node(g, i))
child(g::DAT, i, j)        = child(node(g, j), i)

is_binary(g::DAT, i) = arity(g, i) == 2

node_count(g::DAT)     = g.node_count
variable_count(g::DAT) = g.variable_count
constant_count(g::DAT) = g.constant_count

dependent_subexpression_index(g::DAT, i) = g.dependent_subexpression_dict[i]
dep_subexpr_count(g::DAT)                = length(g.dependent_subexpressions)
sparsity(g::DAT, i)                      = g.sparsity
rev_sparsity(g::DAT, i::Int, k::Int)     = g.rev_sparsity[i]

user_univariate_operator(g::DAT, i) = g.user_operators.univariate_operator_f[i]
user_multivariate_operator(g::DAT, i) = g.user_operators.multivariate_operator_evaluator[i]

function DirectedTree(aux_info, d, op::OperatorRegistry, sub_sparsity::Dict{Int,Vector{Int}}, subexpr_linearity, parameter_values, is_sub, subexpr_indx)

    nd = copy(d.nd)
    adj = copy(d.adj)
    const_values = copy(d.const_values)

    sparsity, dependent_subexpressions = _compute_sparsity(d, sub_sparsity, is_sub, subexpr_indx)
    dependent_subexpression_dict = Dict{Int,Int}()
    for (i, v) in enumerate(dependent_subexpressions)
        dependent_subexpression_dict[v] = i
    end
    rev_sparsity = Dict{Int,Int}()
    for (i,s) in enumerate(sparsity)
        rev_sparsity[s] = i
    end

    nodes = _convert_node_list(aux_info, d.nd, op)
    lin = linearity(nd, adj, subexpr_linearity)
    DirectedTree(nodes = nodes,
                    variables = rev_sparsity,
                    constant_values = const_values,
                    parameter_values = parameter_values,
                    node_count = length(nodes),
                    variable_count = length(sparsity),
                    constant_count = length(const_values),
                    sparsity = sparsity,
                    rev_sparsity = rev_sparsity,
                    dependent_variable_count = length(sparsity),
                    dep_subexpr_count = length(dependent_subexpressions),
                    dependent_subexpressions = copy(dependent_subexpressions),
                    dependent_subexpression_dict = dependent_subexpression_dict,
                    linearity = lin[1],
                    user_operators = op
                    )
end

forward_uni = [i for i in instances(AtomType)]
setdiff!(forward_uni, [VAR_ATOM; PARAM_ATOM; CONST_ATOM; SELECT_ATOM; SUBEXPR])
f_switch = binary_switch(forward_uni, is_forward = true)
@eval function fprop!(t::T, ex::Expression, g::DAT, c::AbstractCache , k::Int) where T<:AbstractCacheAttribute
    id = ex_type(g, k)
    $f_switch
    error("fprop! for ex_type = $id not defined.")
    return
end

reverse_uni = [i for i in instances(AtomType)]
setdiff!(reverse_uni, [VAR_ATOM; PARAM_ATOM; CONST_ATOM; SELECT_ATOM; SUBEXPR])
r_switch = binary_switch(reverse_uni, is_forward = false)
@eval function rprop!(t::T, ex::Expression, g::DAT, c::AbstractCache, k::Int) where T<:AbstractCacheAttribute
    id = ex_type(g, k)
    $r_switch
    error("rprop! for ex_type = $id not defined.")
    return
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
#=
downsample.jl
Copyright 2019-03-05, Jeff Fessler, University of Michigan

see also MIRT/time/downsample.jl
=#

export downsample_dim1
export downsample1
export downsample2
export downsample3
export downsample

using Test: @test, @inferred

"""
    y = downsample_dim1(x, down ; warn::Bool)
downsample by factor `down` along first dimension by averaging

in
- `x [n1 (Nd)]`
- `down::Int` downsampling factor

option
- `warn::Bool` warn if noninteger multiple; default `isinteractive()`

out
- `y [n1/down (Nd)]`

Copyright 2019-03-05, Jeff Fessler, University of Michigan
"""
function downsample_dim1(x::AbstractArray{<:Number}, down::Int
		; warn::Bool = isinteractive())

	dim = size(x)
	dim1 = dim[1]
	x = reshape(x, dim1, :) # [n1 *Nd]
	m1 = Int(floor(dim1 / down))
	if m1 * down < dim1
		warn && @warn("truncating input size $dim1 to $(m1 * down)")
		x = x[1:(m1*down),:]
	end
	y = reshape(x, down, :)
	y = sum(y, dims=1) / down # mean
	y = reshape(y, m1, dim[2:end]...)
#	out = similar(x, (m1, dim[2:end]...))
#	vec(out) .= vec(y) # failed attempt to help @inferred for 2D input arrays
	return y
end

"""
    downsample_dim1_test()
"""
function downsample_dim1_test()
	@inferred downsample_dim1(1:4, 2)
	x = reshape(2:2:48, 4, 6)
#	y = @inferred downsample_dim1(x, 2) # todo: fails
	y = downsample_dim1(x, 2)
	@test y == reshape(3:4:47, 2, 6)
	true
end

"""
    y = downsample1(x, down ; warn=true)
downsample 1D vector by factor `down`

in
- `x [n1]`
- `down::Int` downsampling factor

option
- `warn::Bool` warn if noninteger multiple; default `isinteractive()`

out
- `y [n1/down]`

Copyright 2019-03-05, Jeff Fessler, University of Michigan
"""
function downsample1(x::AbstractVector{<:Number}, down::Int
		; warn::Bool = isinteractive())

	dim = size(x)
	dim1 = dim[1]
	m1 = floor(Int, dim1 / down)
	if m1 * down < dim1
		warn && @warn("truncating input size $dim1 to $(m1 * down)")
		y = reshape((@view x[1:(m1*down)]), down, :)
	else
		y = reshape(x, down, :)
	end
	y = sum(y, dims=1) / down # mean
	y = reshape(y, m1, dim[2:end]...)
	return vec(y)
end

"""
downsample1_test()
"""
function downsample1_test()
	x = 2:2:48
	y = @inferred downsample1(x, 2)
	@test y == 3:4:47
	@inferred downsample1(1:5, 2, warn=false) # warn
	true
end

"""
    y = downsample2(x, down ; warn=true, T)

downsample by averaging by integer factors
in
- `x [nx ny]`
- `down` can be a scalar (same factor for both dimensions) or a 2-vector

option
- `warn::Bool` warn if noninteger multiple; default `isinteractive()`
- `T::DataType` specify output eltype; default `eltype(x[1] / down[1])`

out
- `y [nx/down ny/down]`
"""
function downsample2(x::AbstractMatrix{<:Number},
		down::AbstractVector{Int},
		; warn::Bool = isinteractive(),
		T::DataType = eltype(x[1] / down[1])
	)

	length(down) != 2 && throw("bad down $down")

	idim = size(x)
	odim = floor.(Int, idim ./ down)

	warn && any(odim .* down .!= idim) && @warn("truncating to $odim")

#	y = similar(x, odim) # fails!?
	y = Array{T}(undef, odim[1], odim[2])
	d1 = down[1]
	d2 = down[2]
	for i2=1:odim[2]
		for i1=1:odim[1]
			y[i1,i2] =
				sum(@view x[(i1-1)*d1 .+ (1:d1), (i2-1)*d2 .+ (1:d2)]) / d1 / d2
		end
	end

#=
	fun = (x, d) -> downsample_dim1(x, d, warn=warn)

#=
#	this old way returned an Adjoint type:
	y = fun(x, down[1])
	y = fun(y', down[end])'
=#

#	this way avoids the Adjoint:
	y = fun(x', down[end])' # doing adjoint first
	y = fun(y, down[1])

#=	failed attempt at helping @inferred:
	f1 = (x, d) -> downsample1(x, d, warn=warn)
	tmp = x'
	tmp = hcat([f1(tmp[:,n], down[end]) for n=1:size(tmp,2)]...)'
	y = hcat([f1(tmp[:,n], down[1]) for n=1:size(tmp,2)]...)
@shows y
=#
=#

	return y
end

downsample2(x::AbstractArray{<:Number,2}, down::Int ; args...) =
	downsample2(x, [down, down] ; args...)

"""
downsample2_test()
"""
function downsample2_test()
	x = reshape(0:2:46, 4, 6)
	y = @inferred downsample2(x, 2)
	@test y == [5 21 37; 9 25 41]
	true
end

"""
    y = downsample3(x, down ; warn=true, T)

downsample by averaging by integer factors
in
- `x [nx ny nz]`
- `down` can be a scalar (same factor for all dimensions) or a 3-vector

option
- `warn::Bool` warn if noninteger multiple; default true
- `T::DataType` specify output eltype; default `eltype(x[1] / down[1])`

out
- `y [nx/down ny/down nz/down]`
"""
function downsample3(x::AbstractArray{<:Number,3},
		down::AbstractVector{Int}
		; warn::Bool = isinteractive(),
		T::DataType = eltype(x[1] / down[1])
	)

	length(down) != 3 && throw(DimensionMismatch("down $down"))

	idim = size(x)
	odim = floor.(Int, idim ./ down)

	warn && any(odim .* down .!= idim) && @warn("truncating to $odim")

	return downsample3_perm(x, Tuple(down)) # because it is faster
end

downsample3(x::AbstractArray{<:Number,3}, down::Int ; args...) =
	downsample3(x, [down, down, down] ; args...)

# this method is good for @inferred but is slower!
function downsample3_loop(x::AbstractArray{<:Number,3},
		down::AbstractVector{Int} ;
		T::DataType = eltype(x[1] / down[1]),
	)

	odim = floor.(Int, size(x) ./ down)

	y = Array{T}(undef, odim[1], odim[2], odim[3])
	d1 = down[1]
	d2 = down[2]
	d3 = down[3]
	d123 = d1 * d2 * d3

	for i3=1:odim[3]
		for i2=1:odim[2]
			for i1=1:odim[1]
				y[i1,i2,i3] = sum(@view x[(i1-1)*d1 .+ (1:d1),
					(i2-1)*d2 .+ (1:d2), (i3-1)*d3 .+ (1:d3)]) / d123
			end
		end
	end
	return y
end

# this method fails @inferred but is faster!
function downsample3_perm(x::AbstractArray{<:Number,3}, down::Dims{3})

	# down sample along each dimension
	y = downsample_dim1(x, down[1])
	y = downsample_dim1(permutedims(y, [2, 1, 3]), down[2])
	y = downsample_dim1(permutedims(y, [3, 2, 1]), down[3]) # [3 1 2] order
	y = permutedims(y, [2, 3, 1])

	return y
end

"""
downsample3_test()
"""
function downsample3_test()
	x = [6, 5, 2]
	x = reshape(2*(1:prod(x)), x...)
	for down = 1:2
		y = downsample3(x, down)
	#	y = @inferred downsample3(x, down) # todo fails
		if down == 1
			@test y == x
		elseif down == 2
			@test y == reshape([39 63; 43 67; 47 71], 3,2,1) # squeeze
		end
	end

	x = ones(4, 6, 8)
	@test downsample3_loop(x, [2, 2, 2]) == downsample3_perm(x, (2, 2, 2))

	true
end

#=
if has_aspire # todo
	filex = [test_dir filesep 'testx.fld']
	filey = [test_dir filesep 'testy.fld']
	fld_write(filex, x)
	delete(filey)
	com = ['op sample3 mean ' filey ' ' filex ' %d %d %d']
	com = sprintf(com, down, down, down)
	os_run(com)
	z = fld_read(filey)
	!isequal(y, z) && throw("aspire/matlab mismatch")
end
=#

"""
    downsample(:test)
self test
"""
function downsample(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	@test downsample_dim1_test()
	@test downsample1_test()
	@test downsample2_test()
	@test downsample3_test()
	true
end


text:
#####
##### High level training procedure
#####

"""
    Env{Game, Network, Board}

Type for an AlphZero environment.

The environment features the current neural network, the best neural network
seen so far that is used for data generation, a memory buffer
and an iteration counter.

# Constructor

    Env{Game}(params, curnn, bestnn=copy(curnn), experience=[], itc=0)

Construct a new AlphaZero environment.
- `Game` is the type of the game being played
- `params` has type [`Params`](@ref)
- `curnn` is the current neural network and has type [`AbstractNetwork`](@ref)
- `bestnn` is the best neural network so far, which is used for data generation
- `experience` is the initial content of the memory buffer
   as a vector of [`TrainingSample`](@ref)
- `itc` is the value of the iteration counter (0 at the start of training)
"""
mutable struct Env{Game, Network, Board}
  params :: Params
  curnn  :: Network
  bestnn :: Network
  memory :: MemoryBuffer{Board}
  itc    :: Int
  function Env{Game}(
      params, curnn, bestnn=copy(curnn), experience=[], itc=0) where Game
    Board = GI.Board(Game)
    msize = max(params.mem_buffer_size[itc], length(experience))
    memory = MemoryBuffer{Board}(msize, experience)
    return new{Game, typeof(curnn), Board}(
      params, curnn, bestnn, memory, itc)
  end
end

GameType(env::Env{Game}) where Game = Game

#####
##### Training handlers
#####

"""
    Handlers

Namespace for the callback functions that are used during training.
This enables logging, saving and plotting to be implemented separately.
An example handler object is `Session`.

All callback functions take a handler object `h` as their first argument
and sometimes a second argment `r` that consists in a report.

| Callback                    | Comment                                        |
|:----------------------------|:-----------------------------------------------|
| `iteration_started(h)`      | called at the beggining of an iteration        |
| `self_play_started(h)`      | called once per iter before self play starts   |
| `game_played(h)`            | called after each game of self play            |
| `self_play_finished(h, r)`  | sends report: [`Report.SelfPlay`](@ref)        |
| `memory_analyzed(h, r)`     | sends report: [`Report.Memory`](@ref)          |
| `learning_started(h, r)`    | sends report: [`Report.LearningStatus`](@ref)  |
| `updates_started(h)`        | called before each series of batch updates     |
| `updates_finished(h, r)`    | sends report: [`Report.LearningStatus`](@ref)  |
| `checkpoint_started(h)`     | called before a checkpoint evaluation starts   |
| `checkpoint_game_played(h)` | called after each arena game                   |
| `checkpoint_finished(h, r)` | sends report: [`Report.Checkpoint`](@ref)      |
| `learning_finished(h, r)`   | sends report: [`Report.Learning`](@ref)        |
| `iteration_finished(h, r)`  | sends report: [`Report.Iteration`](@ref)       |
| `training_finished(h)`      | called once at the end of training             |
"""
module Handlers

  import ..Report

  function iteration_started(h)      return end
  function self_play_started(h)      return end
  function game_played(h)            return end
  function self_play_finished(h, r)  return end
  function memory_analyzed(h, r)     return end
  function learning_started(h, r)    return end
  function updates_started(h)        return end
  function updates_finished(h, r)    return end
  function checkpoint_started(h)     return end
  function checkpoint_game_played(h) return end
  function checkpoint_finished(h, r) return end
  function learning_finished(h, r)   return end
  function iteration_finished(h, r)  return end
  function training_finished(h)      return end

end

import .Handlers

#####
##### Public utilities
#####

"""
    get_experience(env::Env)

Return the content of the agent's memory as a
vector of [`TrainingSample`](@ref).
"""
get_experience(env::Env) = get(env.memory)

"""
    initial_report(env::Env)

Return a report summarizing the configuration of agent before training starts,
as an object of type [`Report.Initial`](@ref).
"""
function initial_report(env::Env)
  num_network_parameters = Network.num_parameters(env.curnn)
  num_reg_params = Network.num_regularized_parameters(env.curnn)
  player = MctsPlayer(env.curnn, env.params.self_play.mcts)
  mcts_footprint_per_node = MCTS.memory_footprint_per_node(player.mcts)
  return Report.Initial(
    num_network_parameters, num_reg_params, mcts_footprint_per_node)
end

#####
##### Training loop
#####

function resize_memory!(env::Env{G,N,B}, n) where {G,N,B}
  exp = get(env.memory)
  env.memory = MemoryBuffer{B}(n, exp)
  return
end

function evaluate_network(contender, baseline, params, handler)
  contender = MctsPlayer(contender, params.mcts)
  baseline = MctsPlayer(baseline, params.mcts)
  ngames = params.num_games
  rec = Recorder{GameType(baseline)}()
  avgz = pit(contender, baseline, ngames; memory=rec,
      reset_every=params.reset_mcts_every,
      flip_probability=params.flip_probability) do i, z
    Handlers.checkpoint_game_played(handler)
  end
  redundancy = compute_redundancy(rec)
  return avgz, redundancy
end

function learning!(env::Env, handler)
  ap = env.params.arena
  lp = env.params.learning
  checkpoints = Report.Checkpoint[]
  losses = Float32[]
  tloss, teval, ttrain = 0., 0., 0.
  experience = get(env.memory)
  if env.params.use_symmetries
    experience = augment_with_symmetries(GameType(env), experience)
  end
  trainer, tconvert = @timed Trainer(env.curnn, experience, lp)
  init_status = learning_status(trainer)
  Handlers.learning_started(handler, init_status)
  # Compute the number of batches between each checkpoint
  nbatches = lp.max_batches_per_checkpoint
  if !iszero(lp.min_checkpoints_per_epoch)
    ntotal = num_batches_total(trainer)
    nbatches = min(nbatches, ntotal ÷ lp.min_checkpoints_per_epoch)
  end
  # Loop state variables
  best_evalz = ap.update_threshold
  nn_replaced = false

  for k in 1:lp.num_checkpoints
    # Execute a series of batch updates
    Handlers.updates_started(handler)
    dlosses, dttrain = @timed batch_updates!(trainer, nbatches)
    status, dtloss = @timed learning_status(trainer)
    Handlers.updates_finished(handler, status)
    tloss += dtloss
    ttrain += dttrain
    append!(losses, dlosses)
    # Run a checkpoint evaluation
    Handlers.checkpoint_started(handler)
    env.curnn = get_trained_network(trainer)
    (evalz, redundancy), dteval =
      @timed evaluate_network(env.curnn, env.bestnn, ap, handler)
    teval += dteval
    # If eval is good enough, replace network
    success = (evalz >= best_evalz)
    if success
      nn_replaced = true
      env.bestnn = copy(env.curnn)
      best_evalz = evalz
    end
    checkpoint_report = Report.Checkpoint(
      k * nbatches, status, evalz, redundancy, success)
    push!(checkpoints, checkpoint_report)
    Handlers.checkpoint_finished(handler, checkpoint_report)
  end

  report = Report.Learning(
    tconvert, tloss, ttrain, teval,
    init_status, losses, checkpoints, nn_replaced)
  Handlers.learning_finished(handler, report)
  return report
end

function simple_memory_stats(env)
  mem = get_experience(env)
  nsamples = length(mem)
  ndistinct = length(merge_by_board(mem))
  return nsamples, ndistinct
end

function self_play!(env::Env{G}, handler) where G
  params = env.params.self_play
  player = MctsPlayer(env.bestnn, params.mcts)
  new_batch!(env.memory)
  Handlers.self_play_started(handler)
  mem_footprint = 0
  elapsed = @elapsed begin
    for i in 1:params.num_games
      self_play!(player, env.memory)
      Handlers.game_played(handler)
      reset_every = params.reset_mcts_every
      if (!isnothing(reset_every) && i % reset_every == 0) ||
          i == params.num_games
        mem_footprint = max(mem_footprint,
          MCTS.approximate_memory_footprint(player.mcts))
        MCTS.reset!(player.mcts)
      end
      if !isnothing(params.gc_every) && i % params.gc_every == 0
        Network.gc(env.bestnn)
      end
    end
  end
  MCTS.memory_footprint(player.mcts)
  inference_tr = MCTS.inference_time_ratio(player.mcts)
  speed = last_batch_size(env.memory) / elapsed
  expdepth = MCTS.average_exploration_depth(player.mcts)
  memsize, memdistinct = simple_memory_stats(env)
  report = Report.SelfPlay(
    inference_tr, speed, expdepth, mem_footprint, memsize, memdistinct)
  Handlers.self_play_finished(handler, report)
  return report
end

function memory_report(env::Env, handler)
  if isnothing(env.params.memory_analysis)
    return nothing
  else
    report = memory_report(
      env.memory, env.curnn, env.params.learning, env.params.memory_analysis)
    Handlers.memory_analyzed(handler, report)
    return report
  end
end

"""
    train!(env::Env, handler=nothing)

Start or resume the training of an AlphaZero agent.

A `handler` object can be passed that implements a subset of the callback
functions defined in [`Handlers`](@ref).
"""
function train!(env::Env, handler=nothing)
  while env.itc < env.params.num_iters
    Handlers.iteration_started(handler)
    resize_memory!(env, env.params.mem_buffer_size[env.itc])
    sprep, spperfs = Report.@timed self_play!(env, handler)
    mrep, mperfs = Report.@timed memory_report(env, handler)
    lrep, lperfs = Report.@timed learning!(env, handler)
    rep = Report.Iteration(spperfs, mperfs, lperfs, sprep, mrep, lrep)
    env.itc += 1
    Handlers.iteration_finished(handler, rep)
  end
  Handlers.training_finished(handler)
end


text:
using Genie.Router
using Genie.Renderer
using Genie.Renderer.Html
using VentisController

route("/", VentisController.billgatesbooks, method = POST,named = :index)

# route("/bgbooks", VentisController.billgatesbooks,method = POST,named = :index1)

# route("/bgbooks", BooksController.billgatesbooks)

# route("/") do
#   serve_static_file("welcome.html")
# end

# routes.jl
# route("/hello") do
#   "Welcome to Genie!"
# end

# route("/upload_img", VentisController.upload_img)

# route("/threshold", VentisController.threshold, method = POST, named = :threshold_img)