text:
####################
# Default API URIs #
####################

const API_ENDPOINT = "API_ENDPOINT" in keys(ENV) ?
    HttpCommon.URI(ENV["API_ENDPOINT"]) : HttpCommon.URI("http://104.197.141.88/")

api_uri(path) = HttpCommon.URI(API_ENDPOINT, path = path)

#######################
# GitLab REST Methods #
#######################

function gitlab_request(request_method, endpoint;
                        auth = AnonymousAuth(), handle_error = true,
                        headers = Dict(), params = Dict())
    authenticate_headers!(headers, auth)
    params = gitlab2json(params)
    api_endpoint = api_uri(endpoint)
    @show api_endpoint
    if request_method == Requests.get
        r = request_method(api_endpoint; headers = headers, query = params)
    else
        r = request_method(api_endpoint; headers = headers, json = params)
    end
    ## @show UTF8String(r.data)
    handle_error && handle_response_error(r)
    return r
end

gh_get(endpoint = ""; options...) = gitlab_request(Requests.get, endpoint; options...)
gh_get(endpoint = ""; options...) = gitlab_request(Requests.get, endpoint; options...)
gh_post(endpoint = ""; options...) = gitlab_request(Requests.post, endpoint; options...)
gh_put(endpoint = ""; options...) = gitlab_request(Requests.put, endpoint; options...)
gh_delete(endpoint = ""; options...) = gitlab_request(Requests.delete, endpoint; options...)
gh_patch(endpoint = ""; options...) = gitlab_request(Requests.patch, endpoint; options...)

gh_get_json(endpoint = ""; options...) = Requests.json(gh_get(endpoint; options...))
gh_post_json(endpoint = ""; options...) = Requests.json(gh_post(endpoint; options...))
gh_put_json(endpoint = ""; options...) = Requests.json(gh_put(endpoint; options...))
gh_delete_json(endpoint = ""; options...) = Requests.json(gh_delete(endpoint; options...))
gh_patch_json(endpoint = ""; options...) = Requests.json(gh_patch(endpoint; options...))

#################
# Rate Limiting #
#################

## There are no rate_limit APIs in GitLab
## rate_limit(; options...) = gh_get_json("/rate_limit"; options...)

##############
# Pagination #
##############

has_page_links(r) = haskey(r.headers, "Link")
get_page_links(r) = split(r.headers["Link"], ',')

function find_page_link(links, rel)
    relstr = "rel=\"$(rel)\""
    for i in 1:length(links)
        if contains(links[i], relstr)
            return i
        end
    end
    return 0
end

extract_page_url(link) = match(r"<.*?>", link).match[2:end-1]

function gitlab_paged_get(endpoint; page_limit = Inf, start_page = "", handle_error = true,
                          headers = Dict(), params = Dict(), options...)
    if isempty(start_page)
        r = gh_get(endpoint; handle_error = handle_error, headers = headers, params = params, options...)
    else
        @assert isempty(params) "`start_page` kwarg is incompatible with `params` kwarg"
        r = Requests.get(start_page, headers = headers)
    end
    results = HttpCommon.Response[r]
    page_data = Dict{GitLabString, GitLabString}()
    if has_page_links(r)
        page_count = 1
        while page_count < page_limit
            links = get_page_links(r)
            next_index = find_page_link(links, "next")
            next_index == 0 && break
            r = Requests.get(extract_page_url(links[next_index]), headers = headers)
            handle_error && handle_response_error(r)
            push!(results, r)
            page_count += 1
        end
        links = get_page_links(r)
        for page in ("next", "last", "first", "prev")
            page_index = find_page_link(links, page)
            if page_index != 0
                page_data[page] = extract_page_url(links[page_index])
            end
        end
    end
    return results, page_data
end

function gh_get_paged_json(endpoint = ""; options...)
    results, page_data = gitlab_paged_get(endpoint; options...)
    return mapreduce(Requests.json, vcat, results), page_data
end

##################
# Error Handling #
##################

function handle_response_error(r::HttpCommon.Response)
    if r.status >= 400
        message, docs_url, errors = "", "", ""
        try
            data = Requests.json(r)
            message = get(data, "message", "")
            docs_url = get(data, "documentation_url", "")
            errors = get(data, "errors", "")
        end
        error("Error found in GitLab response:\n",
              "\tStatus Code: $(r.status)\n",
              "\tMessage: $message\n",
              "\tDocs URL: $docs_url\n",
              "\tErrors: $errors")
    end
end


text:
# ---
# title: 84. Largest Rectangle in Histogram
# id: problem84
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Stack
# link: <https://leetcode.com/problems/largest-rectangle-in-histogram/description/>
# hidden: true
# ---
#
# Given _n_ non-negative integers representing the histogram's bar height where
# the width of each bar is 1, find the area of largest rectangle in the
# histogram.
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/histogram.png)
# Above is a histogram where width of each bar is 1, given height =
# `[2,1,5,6,2,3]`.
#
#
#
# ![](https://assets.leetcode.com/uploads/2018/10/12/histogram_area.png)
# The largest rectangle is shown in the shaded area, which has area = `10` unit.
#
#
#
# **Example:**
#
#
#
#     Input: [2,1,5,6,2,3]
#     Output: 10
#
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
struct DijkstraResult{T, U<:Integer}  <: ShortestPathResult
    parents::Vector{U}
    dists::Vector{T}
    predecessors::Vector{Vector{U}}
    pathcounts::Vector{UInt64}
    closest_vertices::Vector{U}
end

"""
    struct Dijkstra <: ShortestPathAlgorithm

The structure used to configure and specify that [`shortest_paths`](@ref)
should use [Dijkstra's algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)
to compute shortest paths. Optional fields for this structure include
- all_paths::Bool - set to `true` to calculate all (redundant, equivalent) paths to a given destination
- track_vertices::Bool - set to `true` to keep a running list of visited vertices (used for specific
  centrality calculations; generally not needed).

`Dijkstra` is the default algorithm used when a distance matrix is specified.

### Implementation Notes
`Dijkstra` supports the following shortest-path functionality:
- non-negative distance matrices / weights
- (optional) multiple sources
- all destinations
- redundant equivalent path tracking
- vertex tracking

### Performance
If using a sparse matrix for `distmx` in [`shortest_paths`](@ref), you *may* achieve better performance
by passing in a transpose of its sparse transpose. That is, assuming `D` is the sparse distance matrix:
```
D = transpose(sparse(transpose(D)))
```
Be aware that realizing the sparse transpose of `D` incurs a heavy one-time penalty, so this strategy
should only be used when multiple calls to [`shortest_paths`](@ref) with the distance matrix are planned.
"""
struct Dijkstra <: ShortestPathAlgorithm
    all_paths::Bool
    track_vertices::Bool
end

Dijkstra(;all_paths=false, track_vertices=false) = Dijkstra(all_paths, track_vertices)

function shortest_paths(g::AbstractGraph, srcs::Vector{U}, distmx::AbstractMatrix{T}, alg::Dijkstra) where {T, U<:Integer}
    nvg = nv(g)
    dists = fill(typemax(T), nvg)
    parents = zeros(U, nvg)
    visited = zeros(Bool, nvg)

    pathcounts = zeros(UInt64, nvg)
    preds = fill(Vector{U}(), nvg)
    H = PriorityQueue{U,T}()
    # fill creates only one array.

    for src in srcs
        dists[src] = zero(T)
        visited[src] = true
        pathcounts[src] = 1
        H[src] = zero(T)
    end

    closest_vertices = Vector{U}()  # Maintains vertices in order of distances from source
    sizehint!(closest_vertices, nvg)

    while !isempty(H)
        u = dequeue!(H)

        if alg.track_vertices
            push!(closest_vertices, u)
        end

        d = dists[u] # Cannot be typemax if `u` is in the queue
        for v in outneighbors(g, u)
            alt = d + distmx[u, v]

            if !visited[v]
                visited[v] = true
                dists[v] = alt
                parents[v] = u

                pathcounts[v] += pathcounts[u]
                if alg.all_paths
                    preds[v] = [u;]
                end
                H[v] = alt
            elseif alt < dists[v]
                dists[v] = alt
                parents[v] = u
                #615
                pathcounts[v] = pathcounts[u]
                if alg.all_paths
                    resize!(preds[v], 1)
                    preds[v][1] = u
                end
                H[v] = alt
            elseif alt == dists[v]
                pathcounts[v] += pathcounts[u]
                if alg.all_paths
                    push!(preds[v], u)
                end
            end
        end
    end

    if alg.track_vertices
        for s in vertices(g)
            if !visited[s]
                push!(closest_vertices, s)
            end
        end
    end

    for src in srcs
        pathcounts[src] = 1
        parents[src] = 0
        empty!(preds[src])
    end

    return DijkstraResult{T, U}(parents, dists, preds, pathcounts, closest_vertices)
end

shortest_paths(g::AbstractGraph, s::Integer, distmx::AbstractMatrix, alg::Dijkstra) = shortest_paths(g, [s], distmx, alg)
# If we don't specify an algorithm, use dijkstra.
shortest_paths(g::AbstractGraph, s, distmx::AbstractMatrix) = shortest_paths(g, s, distmx, Dijkstra())


text:
using Test
import RegisterMismatchCommon

# Tests are in RegisterMismatch and RegisterMismatchCuda


text:
""" Contains default TreeEnsemble parameters, and modifies them with kwargs. """
function TreeEnsemble_defaults(kwargs...)
    d = Dict(:max_depth => 2,
        :cp => 1e-4,
        :hypertol => 0.005, # hyperplane separation tolerance
        :minbucket => 0.02,
        :regression => false)
    set_params!(d; kwargs...)
    return d
end

mutable struct TreeEnsemble
    trees::Array
    weights::Array
    params::Dict
    classes::Union{Nothing, Array{Any}}
    solver

    function TreeEnsemble(solver; kwargs...)
        te = new([],
                 [],
                 TreeEnsemble_defaults(),
                 nothing,
                 solver)
        set_params!(te.params; kwargs...)
        return te
    end
end

set_param!(te::TreeEnsemble, s::Symbol, v::Any) = set_param!(te.params, s, v)
set_params!(te::TreeEnsemble; kwargs...) = set_params!(te.params; kwargs...)
get_param(te::TreeEnsemble, s::Symbol) = get_param(te.params, s)

""" Plants a set number of MIOTrees in a TreeEnsemble. """
function plant_trees(te::TreeEnsemble, n_trees::Int)
    solver = te.solver
    for i = 1:n_trees
        push!(te.trees, MIOTree(solver; te.params...))
    end
end

function fit!(te::TreeEnsemble, method::String, X, Y)
    data = split_data(X, Y, bins = length(te.trees))
    if !get_param(te, :regression)
        te.classes = sort(unique(Y)) # Make sure that all trees have the same classes.
        for mt in te.trees
            mt.classes = te.classes
        end
    end
    @showprogress 1 "Training ensemble of $(length(te.trees)) trees. " for i = 1:length(te.trees)
        fit!(te.trees[i], method, data[i][1], data[i][2])
    end
    return
end

""" Computes the optimal weights for a regressing TreeEnsemble. """
function weigh_trees(te, X, Y)
    get_param(te, :regression) || throw(ErrorException("Can only weight TreeEnsembles for regression."))
    m = JuMP.Model(te.solver)
    @variable(m, w[1:length(te.trees)])
    @constraint(m, sum(w) == 1)
    @variable(m, preds[1:length(Y), 1:length(te.trees)])
    evals = hcat([predict(mt, X) for mt in te.trees]...)
    @constraint(m, preds .== evals * w)
    @objective(m, Min, 1/length(Y)*sum((Y .- preds).^2)) # Minimize squared error
    optimize!(m)
    te.weights = getvalue.(w)
    return
end

function predict(te::TreeEnsemble, X)
    evals = hcat([predict(mt, X) for mt in te.trees]...)
    if get_param(te, :regression)
        if isnothing(te.weights)
            throw(ErrorException("TreeEnsemble must be weighted before prediction. Please use the weigh_trees function."))
        else
            return evals * te.weights
        end
    else
        return [mode(evals[i,:])[1] for i = 1:size(X, 1)] # TODO: find a better way to tie-break?
    end
end

function score(te::TreeEnsemble, X, Y)
    preds = predict(te, X)
    if get_param(te, :regression)
        return 1 - sum((preds .- Y).^2) / sum((preds .- sum(Y)/length(Y)*ones(length(Y))).^2)
    else
        return sum(preds .== Y)/length(Y)
    end
end


text:
#############################################################################
# proximal_level_bundle.jl
# Implements the proximal level bundle method.
# For reference look at:
#############################################################################

export proximal_level_bundle

TOL = 1e-3

# maxiters = 20

# λ = 0.5 * 2^0.5
λ = 1

function proximal_level_bundle(p::Problem, f::LovaszExtAtom)
  if length(get_v(p.objective)) > 1
    error("Cannot optimize functions with multiple variables.")
  end

  # the approximated problem
  q = Problem(p.head, p.objective, p.constraints)
  obj = p.objective + f                          # the objective function

  # warm start
  solve!(q, SCSSolver(verbose=false))
  t = Variable(1)
  q.objective += t
  push!(q.solution.primal, q.solution.primal[end])
  q.solution.primal[end - 1] = 0
  q.solution.status = "not yet solved"
  var = get_v(q.objective.children[1])[1]
  a = greedy_rand(f.func, var.value[:])
  push!(q.constraints, dot(a, var) <= t)
  approx = []                                    # the affine functions
  push!(approx, q.constraints[end].lhs)

  # Step 0: Initialization
  solve!(q, SCSSolver(verbose=false), warmstart = true)
  t.value = evaluate(approx[1])[1]
  upper = evaluate(obj)[1]
  lower = evaluate(q.objective)[1]
  gap = upper - lower
  a = greedy_rand(f.func, var.value[:])
  push!(q.constraints, dot(a, var) <= t)
  push!(approx, q.constraints[end].lhs)

  # level, set as a variable so its value updates automatically
  level = Variable(1)
  fix!(level, λ * lower + (1 - λ) * upper)       # the level set

  # projection center
  size = var.size
  x₀ = Variable(size)                            # the prox center
  fix!(x₀, var.value[:])

  # projection
  proj = Problem(:minimize, norm(var - x₀))
  dotprod = dot(a, var)
  diff = evaluate(obj)[1] - evaluate(dotprod)[1]
  push!(proj.constraints, dotprod + diff - level <= 0)
  pos_test = []
  push!(pos_test, proj.constraints[end].lhs)
  println("diff = $(evaluate(pos_test[1])[1])")

  iters = 0                                      # TODO: update TODO: k(l)
  approx_iters = 0
  optimal = var.value[:]                         # the current optimal solution
  cur_val = 0
  pos(x) = evaluate(x)[1] > -1e-1                # find positive elements

  for i = 1:maxiters
    println("gap = $gap")
    println()
    # Step 1: optimality check
    if gap < TOL
      println(iters)
      println(approx_iters)
      return optimal
      break                                      # optimality reached
    else
      # Step 2: level feasibility check
      solve!(proj, SCSSolver(verbose=false, suppress_warnings=true), warmstart = iters > 0)
      iters += 1
      println("proj.optval = $(proj.optval)")
      if proj.status == :Infeasible
        # Step 3: update lower bound and level
        println("lower bound update")
        solve!(q, SCSSolver(verbose=false), warmstart = true)
        approx_iters += 1
        t.value = minimum(evaluate.(approx))[1]
        lower = evaluate(q.objective)[1]
        level.value = λ * lower + (1 - λ) * upper
      elseif proj.status == :Unbounded
        break
        error("The convex function is not strongly convex")
      else
        # Step 4: update linear constraints
        # start_index = length(proj.solution.dual) - length(approx)
        # pos_index = find(pos_test, proj.solution.dual[(start_index + 1): end]) # update the trial set
        maxx = -10000
        for i = 1:length(pos_test)
          maxx = max(maxx, evaluate(pos_test[i])[1])
        end
        println("maxx = $maxx")
        pos_index = find(pos, pos_test)        # update the trial set
        println("pox_index = $pos_index")
        q.constraints = q.constraints[pos_index]
        approx = approx[pos_index]
        proj.constraints = proj.constraints[pos_index]
        pos_test = pos_test[pos_index]
        a = greedy_rand(f.func, var.value[:])
        push!(q.constraints, dot(a, var) <= t)
        push!(approx, q.constraints[end].lhs)
        dotprod = dot(a, var)
        diff = evaluate(obj)[1] - evaluate(dotprod)[1]
        push!(proj.constraints, dotprod + diff - level <= 0)
        push!(pos_test, proj.constraints[end].lhs)
        maxx1 = -10000
        for i = 1:length(pos_test)
          maxx1 = max(maxx1, evaluate(pos_test[i])[1])
        end
        println("maxx1 = $maxx1")
        # Step 5: add new candidate and update upper bound and level
        x₀.value = var.value[:]                # update the prox centers
        cur_val = evaluate(obj)[1]
        if upper > cur_val
          optimal = var.value[:]
          upper = cur_val
          gap = upper - lower
          level.value = λ * lower + (1 - λ) * upper
        # else
        #   var.value[:] = optimal
        end
      end
    end
  end
end


text:
using Test
using DataFrames

using JMcDM

@testset "Euclidean distance" begin
    @test euclidean([0.0, 1.0, 2.0], [0.0, 1.0, 2.0]) == 0.0
    @test euclidean([0.0, 0.0, 0.0]) == 0.0
    @test euclidean([0.0, 0.0, 1.0]) == 1.0
    @test euclidean([0.0, 0.0, 1.0], [0.0, 0.0, 2.0]) == 1.0
end

@testset "Normalization" begin
    tol = 0.00001
    nz = normalize([1.0, 2.0, 3.0, -1.0, 0.0])
    @test isapprox(nz[1], 0.2581989, atol=tol)
    @test isapprox(nz[2], 0.5163978, atol=tol)
    @test isapprox(nz[3], 0.7745967, atol=tol)
    @test isapprox(nz[4], -0.2581989, atol=tol)
    @test isapprox(nz[5], 0.0000000, atol=tol)
end

@testset "Column min and max vectors" begin
    df = DataFrame()
    df[:,:x] = [0.0, 1.0, 10.0]
    df[:,:y] = [0.0, -1.0, -10.0]
    @test colmins(df) == [0.0, -10.0]
    @test colmaxs(df) == [10.0, 0.0]
end

@testset "Unitize vector" begin
    x = [1.0, 2.0, 3.0, 4.0, 5.0]
    result = x |> unitize |> sum
    @test result == 1.0
end

@testset "Product weights with DataFrame" begin
    df = DataFrame()
    df[:, :x] = [1.0, 2.0, 4.0, 8.0]
    df[:, :y] = [10.0, 20.0, 30.0, 40.0]
    w = [0.60, 0.40]
    result = w * df

    @test result[:, :x] == [0.6, 1.2, 2.4, 4.8]
    @test result[:, :y] == [4.0, 8.0, 12.0, 16.0]
end

@testset "Make Decision Matrix" begin
    m = rand(5, 10)
    df = makeDecisionMatrix(m)

    @test isa(df, DataFrame)
    @test size(df) == (5, 10)
end

@testset "TOPSIS" begin
    tol = 0.00001
    df = DataFrame()
    df[:, :x] = Float64[9, 8, 7]
    df[:, :y] = Float64[7, 7, 8]
    df[:, :z] = Float64[6, 9, 6]
    df[:, :q] = Float64[7, 6, 6]
    w = Float64[4, 2, 6, 8]

    fns = makeminmax([maximum, maximum, maximum, maximum])

    result = topsis(df, w, fns)

    @test isa(result, TopsisResult)
    @test result.bestIndex == 2
    @test isapprox(result.scores, [0.3876870, 0.6503238, 0.0834767], atol=tol)
end

@testset "VIKOR" begin
    tol = 0.00001
    w =  [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
    Amat = [
      100 92 10 2 80 70 95 80 ;
      80  70 8  4 100 80 80 90 ;
      90 85 5 0 75 95 70 70 ;
      70 88 20 18 60 90 95 85
    ]
    dmat = makeDecisionMatrix(Amat)
    result = vikor(dmat, w)

    @test isa(result, VikorResult)
    @test result.bestIndex == 4

    @test isapprox(result.scores[1], 0.1975012087551764, atol=tol)
    @test isapprox(result.scores[2], 0.2194064473270817, atol=tol)
    @test isapprox(result.scores[3], 0.3507643203516215, atol=tol)
    @test isapprox(result.scores[4], -0.16727341435277993, atol=tol)
end

@testset "ELECTRE" begin
    tol = 0.00001
    w =  [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
    Amat = [
      100 92 10 2 80 70 95 80 ;
      80  70 8  4 100 80 80 90 ;
      90 85 5 0 75 95 70 70 ;
      70 88 20 18 60 90 95 85
    ]
    dmat = makeDecisionMatrix(Amat)
    result = electre(dmat, w)

    @test isa(result, ElectreResult)
    @test isa(result.bestIndex, Tuple)
    @test result.bestIndex[1] == 4

    @test isapprox(result.C, [0.36936937,  0.01501502, -2.47347347,  2.08908909], atol=tol)
    @test isapprox(result.D, [0.1914244, -0.1903929,  2.8843076, -2.8853391], atol=tol)
end

@testset "MOORA" begin
    tol = 0.00001
    w =  [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
    Amat = [
      100 92 10 2 80 70 95 80 ;
      80  70 8  4 100 80 80 90 ;
      90 85 5 0 75 95 70 70 ;
      70 88 20 18 60 90 95 85
    ]
    dmat = makeDecisionMatrix(Amat)
    result = moora(dmat, w)

    @test isa(result, MooraResult)
    @test isa(result.bestIndex, Int64)
    @test result.bestIndex == 4

    @test isapprox(result.scores, [0.33159387, 0.29014464, 0.37304311, 0.01926526], atol=tol)
end

@testset "DEMATEL" begin
    tol = 0.00001
    K = [
        0 3 0 2 0 0 0 0 3 0;
        3 0 0 0 0 0 0 0 0 2;
        4 1 0 2 1 3 1 2 3 2;
        4 1 4 0 1 2 0 1 0 0;
        3 2 3 1 0 3 0 2 0 0;
        4 1 4 4 0 0 0 1 1 3;
        3 0 0 0 0 2 0 0 0 0;
        3 0 4 3 2 3 1 0 0 0;
        4 3 2 0 0 1 0 0 0 2;
        2 1 0 0 0 0 0 0 3 0
    ]

    dmat = makeDecisionMatrix(K)

    result = dematel(dmat)

    @test isapprox(result.threshold, 0.062945, atol=tol)

    @test isapprox(result.c, [0.3991458, 0.2261648, 1.0204318, 0.7538625,
                           0.8096760, 0.9780926, 0.2717874,
                           0.9455390, 0.5960514, 0.2937537], atol=tol)

    @test isapprox(result.r, [1.5527024, 0.7251791, 0.8551461, 0.6895615,
                           0.2059141, 0.6790404, 0.1057168,
                           0.3163574, 0.6484014, 0.5164858], atol=tol)

    @test isapprox(result.influenceMatrix,
            [ 0.0  1.0  0.0  1.0  0.0  0.0  0.0  0.0  1.0  0.0;
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0;
            1.0  1.0  0.0  1.0  0.0  1.0  0.0  1.0  1.0  1.0;
            1.0  1.0  1.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0;
            1.0  1.0  1.0  1.0  0.0  1.0  0.0  1.0  0.0  0.0;
            1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0  1.0  1.0;
            1.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0;
            1.0  0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0  0.0;
            1.0  1.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0;
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0], atol=tol)

    @test isapprox(result.weights,
        [0.1686568559124561,
        0.07991375718719543,
        0.14006200243438863,
        0.10748052790517183,
        0.08789022388276985,
        0.12526272598854982,
        0.03067915023486491,
        0.10489168834828348,
        0.092654758940811,
        0.06250830916550884
    ], atol=tol)

end

@testset "AHP - Consistency" begin

    tol = 0.00001

    K = [
        1 7 (1 / 5) (1 / 8) (1 / 2) (1 / 3) (1 / 5) 1;
        (1 / 7) 1 (1 / 8) (1 / 9) (1 / 4) (1 / 5) (1 / 9) (1 / 8);
        5 8 1 (1 / 3) 4 2 1 1;
        8 9 3 1 7 5 3 3;
        2 4 (1 / 4) (1 / 7) 1 (1 / 2) (1 / 5) (1 / 5);
        3 5 (1 / 2) (1 / 5) 2 1 (1 / 3) (1 / 3);
        5 9 1 (1 / 3) 5 3 1 1;
        1 8 1 (1 / 3) 5 3 1 1
    ]

    dmat = makeDecisionMatrix(K)
    result::AHPConsistencyResult = ahp_consistency(dmat)

    @test isa(result, AHPConsistencyResult)

    @test result.isConsistent

    @test isapprox(result.CR, 0.07359957154133831, atol=tol)

    @test isapprox(result.CI, 0.10377539587328702, atol=tol)

    @test isapprox(result.lambda_max, 8.72642777111301, atol=tol)

    @test isapprox(result.pc, [
        8.40982182534543
        8.227923797501001
        8.95201017080425
        8.848075127995905
        8.860427038870013
        8.941497805932498
        8.946071053879708
        8.62559534857526
    ], atol=tol)
end

@testset "AHP" begin

    tol = 0.00001

    K = [
    1 7 1 / 5 1 / 8 1 / 2 1 / 3 1 / 5 1;
    1 / 7 1 1 / 8 1 / 9 1 / 4 1 / 5 1 / 9 1 / 8;
    5 8 1 1 / 3 4 2 1 1;
    8 9 3 1 7 5 3 3;
    2 4 1 / 4 1 / 7 1 1 / 2 1 / 5 1 / 5;
    3 5 1 / 2 1 / 5 2 1 1 / 3 1 / 3;
    5 9 1 1 / 3 5 3 1 1;
    1 8 1 1 / 3 5 3 1 1
]

    A1 = [1 3 1 / 5 2;
                  1 / 3 1 1 / 7 1 / 3;
                  5 7 1 4;
                  1 / 2 3 1 / 4 1]
    A2 = [1 1 / 2 4 5;
                  2 1 6 7;
                  1 / 4 1 / 6 1 3;
                  1 / 5 1 / 7 1 / 3 1]
    A3 = [1 1 / 2 1 / 6 3;
                  2 1 1 / 4 5;
                  6 4 1 9;
                  1 / 3 1 / 5 1 / 9 1]
    A4 = [1 7 1 / 4 2;
                  1 / 7 1 1 / 9 1 / 5;
                  4 9 1 5;
                  1 / 2 5 1 / 5 1]
    A5 = [1 6 2 3;
                  1 / 6 1 1 / 4 1 / 3;
                  1 / 2 4 1 2;
                  1 / 3 3 1 / 2 1]
    A6 = [1 1 / 4 1 / 2 1 / 7;
                  4 1 2 1 / 3;
                  2 1 / 2 1 1 / 5;
                  7 3 5 1]
    A7 = [1 3 7 1;
                  1 / 3 1 4 1 / 3;
                  1 / 7 1 / 4 1 1 / 7;
                  1 3 7 1]
    A8 = [1 2 5 8;
                  1 / 2 1 3 6;
                  1 / 5 1 / 3 1 3;
                  1 / 8 1 / 6 1 / 3 1];

    km = makeDecisionMatrix(K)
    as = map(makeDecisionMatrix, [A1, A2, A3, A4, A5, A6, A7, A8])

    result = ahp(as, km)

    @test isa(result, AHPResult)

    @test result.bestIndex == 3

    @test isapprox(result.scores, [0.2801050, 0.1482273, 0.3813036, 0.1903641], atol=tol)
end

@testset "NDS" begin
    cases = [
        1.0 2.0 3.0;
        2.0 1.0 3.0;
        1.0 3.0 2.0;
        4.0 5.0 6.0
    ]

    nd = makeDecisionMatrix(cases)

    result = nds(nd)

    @test isa(result, NDSResult)

    @test result.bestIndex == 4

    @test result.ranks == [0, 0, 0, 3]

end

@testset "Laplace" begin

    tol = 0.00001

    mat = [
        3000 2750 2500 2250;
        1500 4750 8000 7750;
        2000 5250 8500 11750
    ]

    dm = makeDecisionMatrix(mat)

    result = laplace(dm)

    @test isa(result, LaplaceResult)

    @test isapprox(result.expected_values, [2625.0, 5500, 6875], atol=tol)

    @test result.bestIndex == 3
end

@testset "Maximin" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    dm = makeDecisionMatrix(mat)

    result = maximin(dm)

    @test isa(result, MaximinResult)

    @test isapprox(result.rowmins, [18, 18, 24, 20], atol=tol)

    @test result.bestIndex == 3
end

@testset "Maximax" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    dm = makeDecisionMatrix(mat)

    result = maximax(dm)

    @test isa(result, MaximaxResult)

    @test isapprox(result.rowmaxs, [26, 34, 34, 30], atol=tol)

    @test result.bestIndex in [2, 3]
end

@testset "Minimax" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    dm = makeDecisionMatrix(mat)

    result = minimax(dm)

    @test isa(result, MinimaxResult)

    @test isapprox(result.rowmaxs, [26, 34, 34, 30], atol=tol)

    @test result.bestIndex == 1
end

@testset "Minimin" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    dm = makeDecisionMatrix(mat)

    result = minimin(dm)

    @test isa(result, MiniminResult)

    @test isapprox(result.rowmins, [18, 18,24,20], atol=tol)

    @test result.bestIndex in [1, 2]
end

@testset "Savage" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    dm = makeDecisionMatrix(mat)

    result = savage(dm)

    @test isa(result, SavageResult)

    @test result.bestIndex == 4

end

@testset "Hurwicz" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    dm = makeDecisionMatrix(mat)

    result = hurwicz(dm)

    @test isa(result, HurwiczResult)

    @test result.bestIndex == 3

end

@testset "Maximum likelihood for single criterion" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    weights = [0.2, 0.5, 0.2, 0.1]

    dm = makeDecisionMatrix(mat)

    result = mle(dm, weights)

    @test isa(result, MLEResult)

    @test result.bestIndex == 2

    @test isapprox(result.scores, [24, 29.2, 27, 27], atol=tol)
end

@testset "Expected Regret" begin

    tol = 0.00001

    mat = [
        26 26 18 22;
        22 34 30 18;
        28 24 34 26;
        22 30 28 20
    ]

    weights = [0.2, 0.5, 0.2, 0.1]

    dm = makeDecisionMatrix(mat)

    result = expectedregret(dm, weights)

    @test isa(result, ExpectedRegretResult)

    @test result.bestIndex == 2

    @test isapprox(result.scores, [8, 2.8, 5, 5], atol=tol)
end

@testset "Game" begin

    tol = 0.00001

    mat = [-2 6 3;
        3 -4 7;
        -1 2 4]

    dm = makeDecisionMatrix(mat)
    result = game(dm)

    @test isa(result, GameResult)

    @test result.value == 0.6666666666666661

    @test isapprox(result.row_player_probabilities, [0.4666667, 0.5333333, 0.0000000], atol=tol)

end

@testset "Game - Rock & Paper & Scissors" begin

    tol = 0.00001

    mat = [0 -1 1;
      1 0 -1;
       -1 1 0]

    dm = makeDecisionMatrix(mat)
    result = game(dm)

    @test isa(result, GameResult)

    @test result.value == 0.0

    @test isapprox(result.row_player_probabilities, [0.333333, 0.333333, 0.33333], atol=tol)

end

@testset "Data Envelopment" begin

    tol = 0.00001

    x1 = [96.0, 84, 90, 81, 102, 83, 108, 99, 95]
    x2 = [300.0, 282, 273, 270, 309, 285, 294, 288, 306]

    out = [166.0, 150, 140, 136, 171, 144, 172, 170, 165]
    inp = hcat(x1, x2)

    result::DataEnvelopResult = dataenvelop(inp, out)

    @test isa(result, DataEnvelopResult)

    @test result.orderedcases == [:Case8, :Case2, :Case7, :Case1, :Case9, :Case6, :Case5, :Case4, :Case3]

    @test isapprox(result.efficiencies,
        [0.9879815986198964, 0.9999999999999999,
        0.8959653733189055, 0.9421686746987951,
        0.9659435120753173, 0.9715662650602411,
        0.9911164465786314, 1.0, 0.9841048789857857], atol=tol)

   @test isapprox(result.references[:, :Case1],
        [0.0, 0.544106, 0.0, 0.0, 0.0,
        0.0, 0.0, 0.496377, 0.0], atol=tol)

   @test isapprox(result.references[:, :Case2],
        [0.0, 1, 0, 0, 0, 0, 0, 0, 0], atol=tol)

    @test isapprox(result.references[:, :Case3],
        [0.0, 0.266193, 0, 0, 0, 0, 0, 0.588654, 0], atol=tol)

    @test isapprox(result.references[:, :Case4],
        [0.0, 0.860241, 0, 0, 0, 0, 0, 0.0409639, 0], atol=tol)

    @test isapprox(result.references[:, :Case5],
        [0.0, 0.314982, 0, 0, 0, 0, 0, 0.727957, 0], atol=tol)

    @test isapprox(result.references[:, :Case6],
        [0.0, 0.96 , 0, 0, 0, 0, 0, 0, 0], atol=tol)

    @test isapprox(result.references[:, :Case7],
        [0.0, 0 , 0, 0, 0, 0, 0, 1.01176, 0], atol=tol)

    @test isapprox(result.references[:, :Case8],
        [0.0, 0 , 0, 0, 0, 0, 0, 1.0, 0], atol=tol)

    @test isapprox(result.references[:, :Case9],
        [0.0, 0.774923 , 0, 0, 0, 0, 0, 0.286833, 0], atol=tol)

end

@testset "Grey Relational Analysis" begin

    tol = 0.0001

    df = DataFrame(
        :K1 => [105000.0, 120000, 150000, 115000, 135000],
        :K2 => [105.0, 110, 120, 105, 115],
        :K3 => [10.0, 15, 12, 20, 15],
        :K4 => [4.0, 4, 3, 4, 5],
        :K5 => [300.0, 500, 550, 600, 400],
        :K6 => [10.0, 8, 12, 9, 9]
    )
    functionlist = [minimum, maximum, minimum, maximum, maximum, minimum]

    w = [0.05, 0.20, 0.10, 0.15, 0.10, 0.40]

    result = grey(df, w, functionlist)

    @test isa(result, GreyResult)

    @test isapprox(result.scores, [0.525, 0.7007142857142857,
    0.5464285714285715, 0.5762820512820512, 0.650952380952381], atol=tol)

    @test result.bestIndex == 2
end

@testset "SAW" begin

    @testset "Example 1: 4 criteria × 4 alternatives" begin
        tol = 0.0001
        df = DataFrame(
            :c1 => [25.0, 21, 19, 22],
            :c2 => [65.0, 78, 53, 25],
            :c3 => [7.0, 6, 5, 2],
            :c4 => [20.0, 24, 33, 31]
        )
        weights = [0.25, 0.25, 0.25, 0.25]
        fns = [maximum, maximum, minimum, maximum]
        result = saw(df, weights, fns)

        @test result isa SawResult

        @test isapprox(result.scores, [0.681277, 0.725151, 0.709871, 0.784976], atol=tol)
    end

    @testset "Example 2: 7 criteria × 5 alternatives " begin
        tol = 0.0001
        decmat = [4.0  7  3  2  2  2  2;
                  4.0  4  6  4  4  3  7;
                  7.0  6  4  2  5  5  3;
                  3.0  2  5  3  3  2  5;
                  4.0  2  2  5  5  3  6]

        df = makeDecisionMatrix(decmat)
        weights = [0.283, 0.162, 0.162, 0.07, 0.085, 0.162, 0.076]
        fns = makeminmax([maximum for i in 1:7])
        result = saw(df, weights, fns)

        @test result isa SawResult
        @test isapprox(result.scores, [0.553228, 0.713485, 0.837428, 0.514657, 0.579342], atol=tol)
        @test result.bestIndex == 3
        @test result.ranking == [3, 2, 5, 1, 4]
    end
end

@testset "ARAS Additive Ratio Assessment" begin

    tol = 0.0001

    df = DataFrame(
        :K1 => [105000.0, 120000, 150000, 115000, 135000],
        :K2 => [105.0, 110, 120, 105, 115],
        :K3 => [10.0, 15, 12, 20, 15],
        :K4 => [4.0, 4, 3, 4, 5],
        :K5 => [300.0, 500, 550, 600, 400],
        :K6 => [10.0, 8, 12, 9, 9]
    )
    functionlist = [minimum, maximum, minimum, maximum, maximum, minimum]

    w = [0.05, 0.20, 0.10, 0.15, 0.10, 0.40]

    result = aras(df, w, functionlist)

    @test isa(result, ARASResult)

    @test isapprox(result.scores, [0.81424068, 0.89288620, 0.76415790, 0.84225462, 0.86540635], atol=tol)

    @test result.bestIndex == 2
end

@testset "WPM" begin
    tol = 0.0001
    decmat = [3	12.5	2	120	14	3;
                5	15	3	110	38	4;
                3	13	2	120	19	3;
                4	14	2	100	31	4;
                3	15	1.5	125	40	4]

    df = makeDecisionMatrix(decmat)

    weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201]

    fns = [maximum, minimum, minimum, maximum, minimum, maximum]

    result = wpm(df, weights, fns)

    @test result isa WPMResult

    @test isapprox(result.scores, [0.7975224331331, 0.7532541470585, 0.7647463553356, 0.7873956894791, 0.7674278741782], atol=tol)

end

@testset "WASPAS" begin
    tol = 0.0001
    decmat = [3	12.5	2	120	14	3;
                5	15	3	110	38	4;
                3	13	2	120	19	3;
                4	14	2	100	31	4;
                3	15	1.5	125	40	4]

    df = makeDecisionMatrix(decmat)

    weights = [0.221, 0.159, 0.175, 0.127, 0.117, 0.201]

    fns = [maximum, minimum, minimum, maximum, minimum, maximum]

    result = waspas(df, weights, fns)

    @test result isa WASPASResult

    @test isapprox(result.scores, [0.805021, 0.775060, 0.770181, 0.796424, 0.788239], atol=tol)

end

@testset "EDAS" begin
    tol = 0.0001
    decmat = [5000 5 5300 450;
                4500 5 5000 400;
                4500 4 4700 400;
                4000 4 4200 400;
                5000 4 7100 500;
                5000 5 5400 450;
                5500 5 6200 500;
                5000 4 5800 450]

    df = makeDecisionMatrix(decmat)

    weights = [0.25, 0.25, 0.25, 0.25];

    fns = [maximum, maximum, minimum, minimum];

    result = edas(df, weights, fns)

    @test result isa EDASResult

    @test isapprox(result.scores, [0.759594, 0.886016, 0.697472, 0.739658, 0.059083, 0.731833, 0.641691, 0.385194], atol=tol)

end

@testset "MARCOS" begin

    tol = 0.0001
    decmat = [8.675 8.433 8.000 7.800 8.025 8.043;
    8.825 8.600 7.420 7.463 7.825 8.229;
    8.325 7.600 8.040 7.700 7.925 7.600;
    8.525 8.667 7.180 7.375 7.750 8.071]

    df = makeDecisionMatrix(decmat)

    weights = [0.19019, 0.15915, 0.19819, 0.19019, 0.15115, 0.11111];

    fns = [maximum, maximum, maximum, maximum, maximum, maximum];

    Fns = convert(Array{Function,1}, fns)

    result = marcos(df, weights, Fns)

    @test result isa MARCOSResult

    @test isapprox(result.scores, [0.684865943528, 0.672767106696, 0.662596906139, 0.661103207660], atol=tol)

end

@testset "MABAC" begin

    tol = 0.0001
    decmat = [2 1 4 7 6 6 7 3000;
    4 1 5 6 7 7 6 3500;
    3 2 6 6 5 6 8 4000;
    5 1 5 7 6 7 7 3000;
    4 2 5 6 7 7 6 3000;
    3 2 6 6 6 6 6 3500]

    df = makeDecisionMatrix(decmat)

    weights = [0.293, 0.427, 0.067, 0.027, 0.053, 0.027, 0.053, 0.053]

    fns = [maximum, maximum, maximum, maximum, maximum, maximum, maximum, minimum];

    result = mabac(df, weights, fns)

    @test result isa MABACResult

    @test isapprox(result.scores, [-0.31132, -0.10898, 0.20035, 0.04218, 0.34452, 0.20035], atol=tol)

end

@testset "Reverse minimum & maximum array" begin

    fns = [minimum, maximum, maximum, minimum, maximum]
    revfns = [maximum, minimum, minimum, maximum, minimum]

    @test reverseminmax(fns) == revfns

    @test reverseminmax(revfns) == fns
end

@testset "Make Array of minimum and maximum" begin
    result1 = makeminmax([maximum, maximum, maximum, maximum])

    @test typeof(result1) == Array{Function,1}

    @test typeof(result1[1]([1.0, 2.0, 3.0])) == Float64

    @test result1[1]([1.0, 2.0, 3.0]) == 3.0

end

@testset "MAIRCA" begin

    tol = 0.0001
    decmat = [6.952 8.000 6.649 7.268 8.000 7.652 6.316;
    7.319 7.319 6.604 7.319 8.000 7.652 5.313;
    7.000 7.319 7.652 6.952 7.652 6.952 4.642;
    7.319 6.952 6.649 7.319 7.652 6.649 5.000]

    df = makeDecisionMatrix(decmat)

    weights = [0.172, 0.165, 0.159, 0.129, 0.112, 0.122, 0.140];

    fns = [maximum, maximum, maximum, maximum, maximum, maximum, minimum];

    result = mairca(df, weights, fns)

    @test result isa MAIRCAResult

    @test isapprox(result.scores, [0.1206454, 0.0806646, 0.1458627, 0.1454237], atol=tol)

end

@testset "COPRAS" begin

    tol = 0.0001
    decmat = [2.50 240 57 45 1.10 0.333333;
    2.50 285 60 75 4.00 0.428571;
    4.50 320 100 65 7.50 1.111111;
    4.50 365 100 90 7.50 1.111111;
    5.00 400 100 90 11.00 1.111111;
    2.50 225 60 45 1.10 0.333333;
    2.50 270 57 60 4.00 0.428571;
    4.50 330 100 70 7.50 1.111111;
    4.50 365 100 80 7.50 1.111111;
    5.00 380 110 65 8.00 1.111111;
    2.50 285 65 80 4.00 0.400000;
    4.00 280 75 65 4.00 0.400000;
    4.50 365 102 95 7.50 1.111111;
    4.50 400 102 95 7.50 1.111111;
    6.00 450 110 95 11.00 1.176471;
    6.00 510 110 105 11.00 1.176471;
    6.00 330 140 110 18.50 1.395349;
    2.50 240 65 80 4.00 0.400000;
    4.00 280 75 75 4.00 0.400000;
    4.50 355 102 95 7.50 1.111111;
    4.50 385 102 90 7.50 1.111111;
    5.00 385 114 95 7.50 1.000000;
    6.00 400 110 90 11.00 1.000000;
    6.00 480 110 95 15.00 1.000000;
    6.00 440 140 100 18.50 1.200000;
    6.00 500 140 100 18.50 1.200000;
    5.00 450 125 100 15.00 1.714286;
    6.00 500 150 125 18.50 1.714286;
    6.00 515 180 140 22.00 2.307692;
    7.00 550 200 150 30.00 2.307692;
    6.00 500 180 140 15.00 2.307692;
    6.00 500 180 140 18.50 2.307692;
    6.00 500 180 140 22.00 2.307692;
    7.00 500 180 140 30.00 2.307692;
    7.00 500 200 140 37.00 2.307692;
    7.00 500 200 140 45.00 2.307692;
    7.00 500 200 140 55.00 2.307692;
    7.00 500 200 140 75.00 2.307692]

    df = makeDecisionMatrix(decmat)

    weights = [0.1667, 0.1667, 0.1667, 0.1667, 0.1667, 0.1667];

    fns = [maximum, maximum, maximum, maximum, maximum, minimum];

    result = copras(df, weights, fns)

    @test result isa COPRASResult

    @test isapprox(result.scores, [0.44194, 0.44395, 0.41042, 0.44403, 0.48177, 0.44074, 0.42430, 0.41737, 0.43474, 0.44382, 0.46625, 0.48602, 0.45019, 0.45825, 0.51953, 0.54265, 0.56134, 0.45588, 0.49532, 0.44788, 0.45014, 0.48126, 0.51586, 0.56243, 0.58709, 0.60091, 0.51850, 0.61085, 0.65888, 0.75650, 0.61430, 0.63486, 0.65542, 0.72065, 0.77680, 0.82379, 0.88253, 1.00000] , atol=tol)

end

@testset "PROMETHEE" begin
    tol = 0.005
    decmat = [42.0 35 43 51;
              89 72 92 85;
              14 85 17 40;
              57 60 45 80;
              48 32 43 40;
              71 45 60 85;
              69 40 72 55;
              64 35 70 60]
    df = makeDecisionMatrix(decmat)
    qs = [49, nothing, 45, 30]
    ps = [100, 98, 95, 80]
    weights = [0.25, 0.35, 0.22, 0.18]
    fns = makeminmax([maximum, maximum, maximum, maximum])
    prefs = convert(Array{Function,1}, [prometLinear, prometVShape, prometLinear, prometLinear])

    result = promethee(df, prefs, weights, fns, qs, ps)

    @test result isa PrometheeResult

    @test isapprox(result.scores, [0.07, -0.15, -0.06, -0.05, 0.10, 0.0, 0.03, 0.06], atol=tol)

    @test result.bestIndex == 2
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
# do_print_tests.jl

# Optional: put using statements at the top.
# My personal preference is to have every test file be a runnable script
# This means that if someone was to run just this file, it should work
# This makes it easy to grab a part of your testset and dig into what the error is.
using ExamplePackage, Base.Test

# Write a few tests
@test output_string(:MySymbol) == "MySymbol"
@test output_string(:x) == "x"
@test_throws MethodError output_string("MySymbol")


text:
using Clang
import Xorg_libxcb_jll
import Xorg_xcb_util_jll
import xkbcommon_jll

xcb_include_dir = joinpath(Xorg_libxcb_jll.artifact_dir, "include", "xcb")
xcb_util_dir = joinpath(dirname(dirname(Xorg_xcb_util_jll.libxcb_util_path)), "include", "xcb")
xkb_include_dir = joinpath(xkbcommon_jll.artifact_dir, "include", "xkbcommon")
xkb_headers = joinpath.(Ref(xkb_include_dir), ["xkbcommon.h", "xkbcommon-x11.h"])
xcb_headers = [joinpath.(Ref(xcb_include_dir), ["xkb.h"])..., joinpath(xcb_util_dir, "xcb_event.h")]

# Set up include paths
clang_xcb_includes = [xcb_include_dir, xcb_util_dir]
clang_xkb_includes = [xkb_include_dir]

# Clang arguments
clang_extraargs = ["-v"]

# Test if a header should be wrapped
function wrap_header(top_hdr, cursor_header)
    any(startswith.(Ref(dirname(cursor_header)), vcat(clang_xcb_includes, clang_xkb_includes)))
end

wc_xcb = init(;
                        headers=xcb_headers,
                        output_file=joinpath(@__DIR__, "..", "gen", "xcb_api.jl"),
                        common_file=joinpath(@__DIR__, "..", "gen", "xcb_common.jl"),
                        clang_includes=clang_xcb_includes,
                        clang_args=clang_extraargs,
                        header_wrapped=wrap_header,
                        header_library=x -> basename(x) == "xkb.h" ? "libxcb_xkb" : basename(x) == "xcb_event.h" ? "libxcb_util" : "libxcb",
                        clang_diagnostics=false,
                        )
run(wc_xcb)

wc_xkb = init(;
                        headers=xkb_headers,
                        output_file=joinpath(@__DIR__, "..", "gen", "xkb_api.jl"),
                        common_file=joinpath(@__DIR__, "..", "gen", "xkb_common.jl"),
                        clang_includes=clang_xkb_includes,
                        clang_args=clang_extraargs,
                        header_wrapped=wrap_header,
                        header_library=x -> endswith(x, "xkbcommon-x11.h") ? "libxkbcommon_x11" : "libxkbcommon",
                        clang_diagnostics=true,
                        )
run(wc_xkb)


text:
using PushGP: PushInstruction, inputtypes, outputtypes, neededstacks, hasstacks, split_type_name_in_parts
using PushGP: numoutputs, arity, name, opname, execute!, oldname, oldopname, desc
using PushGP: PushInstructionA1tA1, Lit, Add, names, Sub, Mul, Div, AfromB
using PushGP: Dup, Swap, Flush, StackDepth, Eq, Lt, Gt, Mod, Max, Min
using PushGP: PushInterpreter, stack

struct EmptyInstr <: PushInstruction end
struct EmptyA1tA1{T} <: PushInstructionA1tA1{T} end

# To make the tests architecture-independent
const IntStringL = lowercase(string(Int)) # Since it can be Int64 or Int32
const IntStringU = uppercase(string(Int)) # Since it can be Int64 or Int32

@testset "Instructions: basic" begin
    ei = EmptyInstr()
    @test inputtypes(ei) == DataType[]
    @test outputtypes(ei) == DataType[]
    @test neededstacks(ei) == DataType[]
    @test arity(ei) == 0
    @test numoutputs(ei) == 0
    @test name(ei) == "emptyinstr"
    @test desc(ei) == "emptyinstr"
    @test opname(ei) == "emptyinstr"
    @test oldname(ei) == "EMPTYINSTR"
    @test oldopname(ei) == "EMPTYINSTR"

    interp = PushInterpreter()
    @test hasstacks(interp, ei) == true
    @test execute!(ei, interp) === nothing # Default behavior is to be a No-Op

    eia1a1 = EmptyA1tA1{Int}()
    @test inputtypes(eia1a1) == DataType[Int]
    @test outputtypes(eia1a1) == DataType[Int]
    @test neededstacks(eia1a1) == DataType[Int]
    @test arity(eia1a1) == 1
    @test numoutputs(eia1a1) == 1
    @test name(eia1a1) == (IntStringL * "_emptya1ta1")
    @test desc(eia1a1) == (IntStringL * "_emptya1ta1")
    @test opname(eia1a1) == "emptya1ta1"
    @test oldname(eia1a1) == (IntStringU * ".EMPTYA1TA1")
    @test oldopname(eia1a1) == "EMPTYA1TA1"

    eia1a1_b = EmptyA1tA1{String}()
    @test name(eia1a1_b) == "string_emptya1ta1"
    @test desc(eia1a1_b) == "string_emptya1ta1"
    @test opname(eia1a1_b) == "emptya1ta1"
    @test oldname(eia1a1_b) == "STRING.EMPTYA1TA1"
    @test oldopname(eia1a1_b) == "EMPTYA1TA1"
    @test sort(names(eia1a1_b)) == sort(["STRING.EMPTYA1TA1", "string_emptya1ta1"])
end # @testset "Instructions" begin

@testset "split_type_name_in_parts" begin
    tn, st1, rest = split_type_name_in_parts(PushInstruction)
    @test tn    == "PushInstruction"
    @test st1  === nothing
    @test rest === nothing

    tn, st1, rest = split_type_name_in_parts(Add{Int})
    @test tn   == "Add"
    @test st1  == string(Int)
    @test rest === nothing

    tn, st1, rest = split_type_name_in_parts(AfromB{Int,Float64})
    @test tn   == "AfromB"
    @test st1  == string(Int)
    @test rest == ["Float64"]
end # @testset "Literal instruction" begin

@testset "Literal instruction" begin
    li = Lit(23)
    @test inputtypes(li) == DataType[]
    @test outputtypes(li) == DataType[Int]
    @test neededstacks(li) == DataType[Int]
    @test arity(li) == 0
    @test numoutputs(li) == 1
    @test name(li) == IntStringL * "_lit"
    @test oldname(li) == IntStringU * ".LIT"
    @test desc(li) == "23"

    interp = PushInterpreter()
    execute!(li, interp)
    @test last(stack(interp, Int)) == 23
    @test length(stack(interp, Int)) == 1
end # @testset "Literal instruction" begin

@testset "Add instruction" begin
    i = Add{Int}()
    @test isa(i, Add{Int})
    @test name(i) == IntStringL * "_add"
    @test opname(i) == "add"
    @test oldopname(i) == "+"
    @test oldname(i) == IntStringU * ".+"

    interp = PushInterpreter()
    push!(interp, 1)
    push!(interp, 2)
    execute!(i, interp)
    @test length(stack(interp, Int)) == 1
    @test last(stack(interp, Int)) == 3
end # @testset "Add instructions" begin

@testset "Sub instruction" begin
    i = Sub{Float64}()
    @test isa(i, Sub{Float64})
    @test name(i) == "float64_sub"
    @test opname(i) == "sub"
    @test oldopname(i) == "-"
    @test oldname(i) == "FLOAT64.-"

    interp = PushInterpreter()
    push!(interp, 1.0)
    push!(interp, 2.1)
    execute!(i, interp)
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == (1.0 - 2.1)
end # @testset "Sub instructions" begin

@testset "Mul instruction" begin
    i = Mul{Float64}()
    @test isa(i, Mul{Float64})
    @test name(i) == "float64_mul"
    @test opname(i) == "mul"
    @test oldopname(i) == "*"
    @test oldname(i) == "FLOAT64.*"

    interp = PushInterpreter()
    push!(interp, 1.7)
    push!(interp, 2.1)
    execute!(i, interp)
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == (1.7 * 2.1)

    # But we can also multiply strings since Julia uses this to concat...
    # This instruction not available in normal Push languages so is unique to PushGP.jl (to the best of my knowledge).
    i = Mul{String}()
    @test isa(i, Mul{String})
    @test name(i) == "string_mul"
    @test opname(i) == "mul"
    @test oldname(i) == "STRING.*"

    interp = PushInterpreter()
    push!(interp, "a")
    push!(interp, "bc34")
    execute!(i, interp)
    @test length(stack(interp, String)) == 1
    @test last(stack(interp, String)) == "abc34" # Note order
end # @testset "Mul instructions" begin

@testset "Div instruction" begin
    i = Div{Float64}()
    @test isa(i, Div{Float64})
    @test name(i) == "float64_div"
    @test opname(i) == "div"
    @test oldopname(i) == "/"
    @test oldname(i) == "FLOAT64./"

    interp = PushInterpreter()
    push!(interp, 1.5)
    push!(interp, 2.0)
    execute!(i, interp)
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == (1.5 / 2.0)

    # Note that if denominator is 0.0 the Div instruction is a No-Op
    interp = PushInterpreter()
    push!(interp, 3.1)
    push!(interp, 0.0)
    @test execute!(i, interp) === nothing
    @test length(stack(interp, Float64)) == 2
    @test last(stack(interp, Float64)) == 0.0
end # @testset "Div instructions" begin

@testset "Mod instruction" begin
    i = Mod{Float64}()
    @test name(i) == "float64_mod"
    @test opname(i) == "mod"
    @test oldopname(i) == "%"
    @test oldname(i) == "FLOAT64.%"

    interp = PushInterpreter()
    push!(interp, 4.5)
    push!(interp, 2.0)
    execute!(i, interp)
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == 2.0

    # Note that if denominator is 0.0 the Mod instruction is a No-Op
    interp = PushInterpreter()
    push!(interp, 3.1)
    push!(interp, 0.0)
    @test execute!(i, interp) === nothing
    @test length(stack(interp, Float64)) == 2
    @test last(stack(interp, Float64)) == 0.0
end # @testset "Mod instructions" begin

@testset "Max instructions" begin
    i = Max{Float64}()
    @test name(i) == "float64_max"
    @test opname(i) == "max"
    @test oldopname(i) == "MAX"
    @test oldname(i) == "FLOAT64.MAX"

    interp = PushInterpreter()
    push!(interp, 3.1)
    push!(interp, 17.54)
    execute!(i, interp)
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == 17.54
end # @testset "Max instruction" begin

@testset "Max instructions" begin
    i = Min{Float64}()
    @test name(i) == "float64_min"
    @test opname(i) == "min"
    @test oldopname(i) == "MIN"
    @test oldname(i) == "FLOAT64.MIN"

    interp = PushInterpreter()
    push!(interp, 3.1)
    push!(interp, 17.54)
    execute!(i, interp)
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == 3.1
end # @testset "Max instruction" begin

@testset "AfromB instructions" begin
    i = AfromB{String,Float64}()
    @test name(i) == "string_from_float64"
    @test oldname(i) == "STRING.FROMFLOAT64"
    @test opname(i) == "from"

    interp = PushInterpreter()
    push!(interp, 3.0)
    execute!(i, interp)
    @test length(stack(interp, String)) == 1
    @test last(stack(interp, String)) == "3.0"

    i = AfromB{Float64,Int}()
    interp = PushInterpreter()
    push!(interp, 5)
    execute!(i, interp)
    @test length(stack(interp, Int)) == 0
    @test length(stack(interp, Float64)) == 1
    @test last(stack(interp, Float64)) == 5.0
end

@testset "Eq instructions" begin
    i = Eq{Float64}()
    @test name(i)         == "float64_eq"
    @test oldname(i)      == "FLOAT64.="
    @test opname(i)       == "eq"
    @test arity(i)        == 2
    @test inputtypes(i)   == DataType[Float64, Float64]
    @test outputtypes(i)  == DataType[Bool]
    @test neededstacks(i) == DataType[Float64, Bool]

    interp = PushInterpreter()
    push!(interp, 5.0)
    push!(interp, 5.0)
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 1
    @test last(stack(interp, Bool)) == true
    @test length(stack(interp, Float64)) == 0

    push!(interp, 1.1)
    push!(interp, 1.2)
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 2
    @test last(stack(interp, Bool)) == false
    @test length(stack(interp, Float64)) == 0
end # @testset "Eq instructions" begin

@testset "Lt instructions" begin
    i = Lt{Int}()
    @test name(i)         == IntStringL * "_lt"
    @test oldname(i)      == IntStringU * ".<"
    @test opname(i)       == "lt"
    @test oldopname(i)    == "<"
    @test arity(i)        == 2
    @test inputtypes(i)   == DataType[Int, Int]
    @test outputtypes(i)  == DataType[Bool]
    @test neededstacks(i) == DataType[Int, Bool]

    interp = PushInterpreter()
    push!(interp, 2)
    push!(interp, 3)
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 1
    @test last(stack(interp, Bool)) == true
    @test length(stack(interp, Int)) == 0

    push!(interp, 3)
    push!(interp, 2)
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 2
    @test last(stack(interp, Bool)) == false
    @test length(stack(interp, Int)) == 0
end # @testset "Lt instructions" begin

@testset "Gt instructions" begin
    i = Gt{String}()
    @test name(i)         == "string_gt"
    @test oldname(i)      == "STRING.>"
    @test opname(i)       == "gt"
    @test oldopname(i)    == ">"
    @test arity(i)        == 2
    @test inputtypes(i)   == DataType[String, String]
    @test outputtypes(i)  == DataType[Bool]
    @test neededstacks(i) == DataType[String, Bool]

    interp = PushInterpreter()
    push!(interp, "c")
    push!(interp, "a")
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 1
    @test last(stack(interp, Bool)) == true
    @test length(stack(interp, String)) == 0

    push!(interp, "a")
    push!(interp, "d")
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 2
    @test last(stack(interp, Bool)) == false
    @test length(stack(interp, String)) == 0
end # @testset "Lt instructions" begin

@testset "Dup stack instruction" begin
    i = Dup{String}()
    @test name(i) == "string_dup"
    @test oldname(i) == "STRING.DUP"
    @test opname(i) == "dup"

    interp = PushInterpreter()
    push!(interp, "42")
    execute!(i, interp)
    @test length(stack(interp, String)) == 2
    @test last(stack(interp, String)) == "42"
end

@testset "Swap stack instruction" begin
    i = Swap{Bool}()
    @test name(i) == "bool_swap"
    @test oldname(i) == "BOOL.SWAP"
    @test opname(i) == "swap"

    interp = PushInterpreter()
    push!(interp, true)
    push!(interp, false)
    execute!(i, interp)
    @test length(stack(interp, Bool)) == 2
    @test last(stack(interp, Bool)) == true # Before swap true was next to last now it is last
    @test stack(interp, Bool)[end-1] == false
end

@testset "Flush stack instruction" begin
    i = Flush{Int}()
    @test name(i) == IntStringL * "_flush"
    @test oldname(i) == IntStringU * ".FLUSH"
    @test opname(i) == "flush"

    interp = PushInterpreter()
    push!(interp, 1)
    push!(interp, 2)
    push!(interp, 3)
    execute!(i, interp)
    @test length(stack(interp, Int)) == 0
end

@testset "StackDepth stack instruction" begin
    i = StackDepth{Int}()
    @test name(i) == IntStringL * "_stackdepth"
    @test oldname(i) == IntStringU * ".STACKDEPTH"
    @test opname(i) == "stackdepth"

    interp = PushInterpreter()
    push!(interp, 1)
    push!(interp, 2)
    push!(interp, 3)
    execute!(i, interp)
    @test length(stack(interp, Int)) == 4
    @test last(stack(interp, Int)) == 3 # Since length was 3 when we executed instruction
end


text:
################## random samples of code using package multiobjectiveGraphs #######################

multObj = multiobj(2,1,2, typeofvalues=Float32)

println(multObj)
println(descr(multObj))

#=
for i in Squares(7)
           println(i)
       end
=#

#=
for obj in multObj
	println(descr(obj))
end
=#

# println(descr(multObj[2]))

max_stars_repo_name: DRobbes/multiobjectiveGraphs

max_issues_repo_name: DRobbes/multiobjectiveGraphs

max_forks_repo_name: DRobbes/multiobjectiveGraphs


text:
##
using Revise
using GLMakie
using Pkg
Pkg.activate(".")
using ProteinStructureReconstruction
using HDF5
using Images
##
filepath = "./Data/Fed_X63_Z3_SIM.h5"
h5file = h5open(filepath, "r")

x_range, y_range = 2000:2300, 1600:1900;
img_data = h5file["Data"][x_range, y_range, :, :] |> x -> reinterpret(N0f16, x);
channels = read(h5file["Data"]["channels"])
##
import ImageBinarization: binarize, MinimumError
img_eq = equalize_zstack(img_data);
##
using Interpolations
function interpolate_zstack(zstack::AbstractArray{T, 3}, pixel_sizes) where {T}
    scale_x, scale_y, scale_z = pixel_sizes/minimum(pixel_sizes)
    nx, ny, nz = size(zstack)
    grids = ((0:nx-1)*scale_x, (0:ny-1)*scale_y, (0:nz-1)*scale_z)
    algorithm = BSpline(Linear())
    itp = interpolate(zstack, algorithm)
    sitp = scale(itp, grids...)
    return sitp(0:(nx-1)*scale_x, 0:(ny-1)*scale_y, 0:(nz-1)*scale_z)
end

pixel_sizes = read(h5file["Data"]["pixel_sizes"])[end:-1:1]
img_interp = mapslices(x -> interpolate_zstack(x, pixel_sizes), img_eq, dims = [1, 2, 3]);
size(img_interp)
##
img_mask = mapslices(x -> binarize(x, MinimumError()), img_interp, dims = [1, 2, 3]) |> x -> float32.(x);
##
using StaticArrays
function pointcloud_from_image(image::AbstractArray{T, N}) where {T, N}
    coords = map(findall(x -> x > zero(T), image)) do ind
        ind |> Tuple |> SVector{N, T}
    end
    return coords
end

img_coords = [pointcloud_from_image(vol) for vol in eachslice(img_mask, dims = 4)];
size(img_coords[1])
##
using NearestNeighbors
trees = map(KDTree, img_coords)

##
# for each point in desmin, find closest distance to actin set
_, desmin_point_actin_set_distance = nn(trees[1], img_coords[2])
##
maximum(desmin_point_actin_set_distance)
##
let
    fig = Figure()
    ax = GLMakie.Axis(fig[1, 1])
    hist!(ax, desmin_point_actin_set_distance*pixel_sizes[1], bins = 30)
    fig
end
##
import StatsBase: fit, Histogram
bin_size = 0.8
distances = collect(0.0:bin_size:maximum(desmin_point_actin_set_distance)+bin_size)
h = fit(Histogram, desmin_point_actin_set_distance, distances, closed=:left)
##
let
    fig = Figure(resolution = (1200, 600))
    ax1 = GLMakie.Axis(fig[1, 1], title = "Point(desmin)-to-Set(actin) correlation function (Δr = 0.8)", xlabel = "Distance (pixels)", ylabel = "Correlation (a.u.)")
    ax2 = GLMakie.Axis(fig[1, 2], yscale = log10, title = "Point(desmin)-to-Set(actin) correlation function (Δr = 0.8)", xlabel = "Distance (pixels)", ylabel = "Correlation (a.u.)")
    scatterlines!(ax1, distances[1:end-1], h.weights)
    scatterlines!(ax2, distances[1:end-1], h.weights)
    fig
end
##
pixel_sizes


text:
module FormatStreams

using Formats
using Formats:
		FormatHandler, Formatted, FormattedIO, FormattedFilename,
		resolveformat, resolvecoding

export FormattedStream, streamf, eachval, eachval!

"""
	FormattedStream{T}

Base type for formatted streams.

Parameter `T` is the type of the objects contained in the stream.
"""
abstract type FormattedStream{T} <: IO end

include("registry.jl")
include("io.jl")

"""
Read, write and detect streams of formatted data

```julia
using Formats        # Framework for data formats/codings (separate package)
using FormatCodecs   # Codecs for common codings (separate package)
using FormatStreams  # Streams of formatted data (this package)
```

To open a file as a stream of formatted objects, automatically inferring its
format/coding if unspecified:

```julia
s = streamf("trajectory.xtc")
```

An existing IO stream can be wrapped as a stream of formatted objects:

```julia
io = open("trajectory.xtc")
s = streamf(io)
```

Format/coding can be specified using functions from the `Formats` package:

```julia
s = streamf(specify("lysozyme.dat", "trajectory/x-trr"))
```

The `FormattedStream` returned by `streamf` can be used with standard IO
functions. Basic functions are always supported:

```julia
s = streamf("trajectory.xtc")
position(s)
eof(s)
seekstart(s)
close(s)
```

The `read` function is typically supported, allowing sequential reading:

```julia
s = streamf("trajectory.xtc")
while !eof(s)
	frame = read(s)
	[...]
end
close(s)
```

Function `read!` may be supported to read into a pre-allocated output buffer:

```julia
s = streamf("trajectory.xtc")
frame = MolecularModel()
while !eof(s)
	read!(s, frame)
	[...]
end
close(s)
```

Some formatted streams are seekable:

```julia
s = streamf("trajectory.trr")
seek(s, 500)
seekend(s)
```

Some have a known length:

```julia
n = length(s1)
```

Some allow replacing and/or appending values with `write`:

```julia
s1 = streamf("input.xtc")
s2 = streamf(openf("output.trr", "w"))
frame = MolecularModel()
while !eof(s1)
	read!(s1, frame)
	write(s2, frame)
end
close(s1)
close(s2)
```

And some can be truncated:

```julia
truncate(s, 5)
```

The `streamf` function supports `do` block syntax, automatically closing the
formatted stream upon completion:

```julia
streamf("trajectory.xtc") do s
	frame = MolecularModel()
	while !eof(s)
		read!(s, frame)
		[...]
	end
end
```

The convenience function `eachval` iterates over all values in a formatted
streamf:

```julia
s = streamf("trajectory.xtc")
for frame in eachval(s)
	[...]
end
```

The `eachval!` function iterates over a stream by reading into a pre-allocated
output buffer. Combined with a `do` block, it provides a compact and efficient
way to iterate over formatted objects:

```julia
streamf("trajectory.xtc") do s
	for frame in eachval!(s, MolecularModel())
		[...]
	end
end
```

Complete documentation is provided in `README.md`.
"""
FormatStreams

end # module


text:
const known_vargroups = Dict(
"Atmosphere"=>[
            "aerosol_optical_thickness_1610",
            "aerosol_optical_thickness_550",
            "aerosol_optical_thickness_555",
            "aerosol_optical_thickness_659",
            "aerosol_optical_thickness_865",
            "air_temperature_2m",
            "ozone",
            "potential_evaporation",
            "precipitation",
            "water_vapour"],
"_Atmosphere"=>["Rg",
            "air_temperature_2m",
            "air_temperature_max",
            "air_temperature_min",
            "cloud_cover",
            "relative_humidity",
            "potential_evaporation",
            "precipitation",
            "vapour_pressure_deficit",
            "wind"],
"Biosphere"=>["bare_soil_evaporation",
            "black_sky_albedo",
            "burnt_area",
            "c_emissions",
            "evaporation",
            "evaporative_stress",
            "gross_primary_productivity",
            "interception_loss",
            "land_surface_temperature",
            "latent_energy",
            "net_ecosystem_exchange",
            "open_water_evaporation",
            "root_moisture",
            "sensible_heat",
            "soil_moisture",
            "surface_moisture",
            "terrestrial_ecosystem_respiration",
            "transpiration",
            "white_sky_albedo"],
"_Biosphere"=>["bare_soil_evaporation",
            "black_sky_albedo",
            "burnt_area",
            "c_emissions",
            "evaporation",
            "evaporative_stress",
            "fpar_fluxcom",
            "gross_primary_productivity",
            "interception_loss",
            "latent_energy",
            "net_ecosystem_exchange",
            "open_water_evaporation",
            "root_moisture",
            "sensible_heat",
            "snow_sublimation",
            "surface_moisture",
            "terrestrial_ecosystem_respiration",
            "transpiration",
            "white_sky_albedo"],
"Fluxcom"=>["gross_primary_productivity",
            "net_ecosystem_exchange",
            "terrestrial_ecosystem_respiration"],
"GLEAM"=>  ["evaporation",
            "evaporative_stress",
            "potential_evaporation",
            "interception_loss",
            "root_moisture",
            "surface_moisture",
            "bare_soil_evaporation",
            "snow_sublimation",
            "transpiration",
            "open_water_evaporation"]
)


text:
## internals.jl : stores estimator internals
#
# A type that stores estimator internals for specific Estimators.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

abstract type AbstractEstimatorInternals end

# estimator internals
struct EstimatorInternals{T1, T2, T3} <: AbstractEstimatorInternals
    default_internals::T1
    index_set_internals::T2
    sample_method_internals::T3
end

EstimatorInternals(index_set::AbstractIndexSet, sample_method::AbstractSampleMethod, options) =
EstimatorInternals(
                   DefaultInternals(index_set, sample_method, options),
                   IndexSetInternals(index_set, sample_method, options),
                   SampleMethodInternals(index_set, sample_method, options)
)

# default internals are shared by all Estimators
struct DefaultInternals{T1, T2, T3, T4, T5} <: AbstractEstimatorInternals
    samples::T1
    samples_diff::T1
    nb_of_samples::T2
    total_work::T3
    total_time::T3
    current_index_set::T4
    index_set_size::T5
end

function DefaultInternals(index_set, sample_method, options)
    max_search_space = index_set isa Union{AD, U} ? options[:max_search_space] : index_set
    indices = get_index_set(max_search_space, options[:max_index_set_param])

    sz = maximum(indices)
    m = options[:nb_of_qoi]
    n = sample_method isa MC ? 1 : maximum(options[:nb_of_shifts].(collect(indices)))

    T = Float64
    samples = [[Vector{T}(undef, 0) for index in CartesianIndices(sz + one(sz))] for i in 1:m, j in 1:n]
    samples_diff = [[Vector{T}(undef, 0) for index in CartesianIndices(sz + one(sz))] for i in 1:m, j in 1:n]

    nb_of_samples = [0 for index in CartesianIndices(sz + one(sz))]
    total_work = [0. for index in CartesianIndices(sz + one(sz))]
    total_time = [0. for index in CartesianIndices(sz + one(sz))]
    current_index_set = Set{Index{ndims(index_set)}}()
    index_set_size = IndexSetSize(0, 0)

    DefaultInternals(samples, samples_diff, nb_of_samples, total_work, total_time, current_index_set, index_set_size)
end

samples(estimator::Estimator) = estimator.internals.default_internals.samples

samples(estimator::Estimator, n_qoi, n_shift) = estimator.internals.default_internals.samples[n_qoi, n_shift]

samples(estimator::Estimator, n_qoi) = samples(estimator, n_qoi, 1)

samples(estimator::Estimator, n_qoi, n_shift, index::Index) = estimator.internals.default_internals.samples[n_qoi, n_shift][index + one(index)]

samples(estimator::Estimator, n_qoi, index::Index) = samples(estimator, n_qoi, 1, index)

append_samples!(estimator::Estimator, n_qoi, n_shift, index::Index, samples_to_append) = append!(estimator.internals.default_internals.samples[n_qoi, n_shift][index + one(index)], samples_to_append)

append_samples!(estimator::Estimator, n_qoi, index::Index, samples_to_append) = append_samples!(estimator, n_qoi, 1, index, samples_to_append)

samples_diff(estimator::Estimator) = estimator.internals.default_internals.samples_diff

samples_diff(estimator::Estimator, n_qoi, n_shift) = estimator.internals.default_internals.samples_diff[n_qoi, n_shift]

samples_diff(estimator::Estimator, n_qoi) = samples_diff(estimator, n_qoi, 1)

samples_diff(estimator::Estimator, n_qoi, n_shift, index::Index) = estimator.internals.default_internals.samples_diff[n_qoi, n_shift][index + one(index)]

samples_diff(estimator::Estimator, n_qoi, index::Index) = samples_diff(estimator, n_qoi, 1, index)

append_samples_diff!(estimator::Estimator, n_qoi, n_shift, index::Index, samples_to_append) = append!(estimator.internals.default_internals.samples_diff[n_qoi, n_shift][index + one(index)], samples_to_append)

append_samples_diff!(estimator::Estimator, n_qoi, index::Index, samples_to_append) = append_samples_diff!(estimator, n_qoi, 1, index, samples_to_append)

has_samples_at_index(estimator::Estimator, index::Index) = !isempty(estimator.internals.default_internals.samples_diff[1][index + one(index)])

nb_of_samples(estimator::Estimator) = estimator.internals.default_internals.nb_of_samples

nb_of_samples(estimator::Estimator, index::Index) = estimator.internals.default_internals.nb_of_samples[index + one(index)]

add_to_nb_of_samples(estimator::Estimator, index::Index, amount::Integer) = estimator.internals.default_internals.nb_of_samples[index + one(index)] += amount

work(estimator::Estimator, index::Index) = total_work(estimator, index)/nb_of_samples(estimator, index)

total_work(estimator::Estimator) = estimator.internals.default_internals.total_work

total_work(estimator::Estimator, index::Index) = estimator.internals.default_internals.total_work[index + one(index)]

add_to_total_work(estimator::Estimator, index::Index, amount::Real) = estimator.internals.default_internals.total_work[index + one(index)] += estimator[:cost_model] isa EmptyFunction ? 0. : amount*estimator[:cost_model](index)

time(estimator::Estimator, index::Index) = total_time(estimator, index)/nb_of_samples(estimator, index)

total_time(estimator::Estimator) = estimator.internals.default_internals.total_time

total_time(estimator::Estimator, index::Index) = estimator.internals.default_internals.total_time[index + one(index)]

add_to_total_time(estimator::Estimator, index::Index, amount::Real) = estimator.internals.default_internals.total_time[index + one(index)] += amount

current_index_set(estimator::Estimator) = estimator.internals.default_internals.current_index_set

keys(estimator::Estimator) = sort(collect(current_index_set(estimator)))

all_keys(estimator::Estimator) = broadcast(i -> i - one(i), filter(i -> has_samples_at_index(estimator, i -one(i)), CartesianIndices(samples(estimator, 1))))

push!(estimator::Estimator, index::Index) = push!(estimator.internals.default_internals.current_index_set, index)

clear(estimator::Estimator) = empty!(current_index_set(estimator))

# specific internals related to the sample method
abstract type SampleMethodInternals <: AbstractEstimatorInternals end

struct MCInternals <: SampleMethodInternals end

SampleMethodInternals(index_set, sample_method, options) = MCInternals()

struct QMCInternals{T1} <: SampleMethodInternals
    generators::T1
end

generator(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index) = estimator.internals.sample_method_internals.generators[index + one(index)]
generator(estimator::Estimator{<:AbstractIndexSet, <:QMC}, index, shift) = estimator.internals.sample_method_internals.generators[index + one(index)][shift]

function SampleMethodInternals(index_set, sample_method::QMC, options)
    max_search_space = index_set isa Union{AD, U} ? options[:max_search_space] : index_set
    indices = get_index_set(max_search_space, options[:max_index_set_param])

    sz = maximum(indices)
    R = CartesianIndices(sz + one(sz))
    generators = [[ShiftedLatticeRule(options[:point_generator]) for i in 1:options[:nb_of_shifts](I)] for I in R]

    QMCInternals(generators)
end

# specific internals related to the index set
abstract type IndexSetInternals <: AbstractEstimatorInternals end

struct GenericIndexSetInternals <: SampleMethodInternals end

IndexSetInternals(index_set, sample_method, options) = GenericIndexSetInternals()

struct ADInternals{T1, T2, T3, T4} <: SampleMethodInternals
    old_set::T1
    active_set::T1
    max_index_set::T1
    boundary::T2
    max_search_space::T3
	logbook::T4
end

function IndexSetInternals(index_set::AD, sample_method, options)
	d = ndims(index_set)
	s_type = Set{Index{d}}
    old_set = s_type()
    active_set = s_type()
    max_index_set = s_type()
    boundary = copy(active_set)
	push!(boundary, Index(ntuple(i -> 0, d)))
    max_search_space = options[:max_search_space]
	logbook = Vector{Tuple{s_type, s_type, Index{d}}}(undef, 0)

    ADInternals(old_set, active_set, max_index_set, boundary, max_search_space, logbook)
end

active_set(estimator::Estimator{<:AD}) = estimator.internals.index_set_internals.active_set

add_to_active_set(estimator::Estimator{<:AD}, index::Index) = push!(estimator.internals.index_set_internals.active_set, index)

remove_from_active_set(estimator::Estimator{<:AD}, index::Index) = delete!(estimator.internals.index_set_internals.active_set, index)

old_set(estimator::Estimator{<:AD}) = estimator.internals.index_set_internals.old_set

add_to_old_set(estimator::Estimator{<:AD}, index::Index) = push!(estimator.internals.index_set_internals.old_set, index)

is_admissable(estimator::Estimator{<:AD}, index::Index) = is_admissable(estimator.internals.index_set_internals.old_set, index)

max_index_set(estimator::Estimator{<:AD}) = estimator.internals.index_set_internals.max_index_set

add_to_max_index_set(estimator::Estimator{<:AD}, index::Index) = push!(estimator.internals.index_set_internals.max_index_set, index)

boundary(estimator::Estimator{<:AD}) = collect(estimator.internals.index_set_internals.boundary)

update_boundary(estimator::Estimator{<:AD}) = begin
    empty!(estimator.internals.index_set_internals.boundary)
	union!(estimator.internals.index_set_internals.boundary, copy(estimator.internals.index_set_internals.active_set))
end

logbook(estimator::Estimator{<:AD}) = estimator.internals.index_set_internals.logbook

log_adaptive_index_set(estimator::Estimator{<:AD}, max_index) = push!(estimator.internals.index_set_internals.logbook, tuple(copy(old_set(estimator)), copy(active_set(estimator)), max_index))

clear(estimator::Estimator{<:AD}) = begin
    empty!(current_index_set(estimator))
    empty!(active_set(estimator))
    empty!(old_set(estimator))
    empty!(max_index_set(estimator))
	empty!(logbook(estimator))
end

struct UInternals{T1, T2} <: SampleMethodInternals
    accumulator::T1
    pmf::T2
end

function IndexSetInternals(index_set::U, sample_method, options)
    indices = get_index_set(options[:max_search_space], options[:max_index_set_param])

    m = options[:nb_of_qoi]
    n = sample_method isa MC ? 1 : maximum(options[:nb_of_shifts].(collect(indices)))
    d = ndims(index_set)

    accumulator = [Vector{Float64}(undef, 0) for i in 1:m, j in 1:n]
    pmf = Dict(index => prod(broadcast(i -> Geometric(1 - exp(-1.5), index[i]), 1:d)) for index in indices)
    normalize!(pmf)

    UInternals(accumulator, pmf)
end

accumulator(estimator::Estimator{<:U}) = estimator.internals.index_set_internals.accumulator

accumulator(estimator::Estimator{<:U}, n_qoi, n_shift) = estimator.internals.index_set_internals.accumulator[n_qoi, n_shift]

accumulator(estimator::Estimator{<:U}, n_qoi) = estimator.internals.index_set_internals.accumulator[n_qoi, 1]

append_to_accumulator!(estimator::Estimator{<:U}, n_qoi, n_shift, samples_to_append) = append!(estimator.internals.index_set_internals.accumulator[n_qoi, n_shift], samples_to_append)

append_to_accumulator!(estimator::Estimator{<:U}, n_qoi, samples_to_append) = append_to_accumulator!(estimator, n_qoi, 1, samples_to_append)

add_to_accumulator!(estimator::Estimator{<:U}, n_qoi, n_shift, samples_to_add) = estimator.internals.index_set_internals.accumulator[n_qoi, n_shift][end-length(samples_to_add)+1:end] .+= samples_to_add

add_to_accumulator!(estimator::Estimator{<:U}, n_qoi, samples_to_add) = add_to_accumulator!(estimator, n_qoi, 1, samples_to_add)

pmf(estimator::Estimator{<:U}) = estimator.internals.index_set_internals.pmf

pmf(estimator::Estimator{<:U}, index::Index) = estimator.internals.index_set_internals.pmf[index]

set_pmf_key(estimator::Estimator{<:U}, index, val) = estimator.internals.index_set_internals.pmf[index] = val

# keeps track of max index size parameter sz
mutable struct IndexSetSize{N}
    sz::N
    max_sz::N
end

sz(estimator::Estimator) = estimator.internals.default_internals.index_set_size.sz
max_sz(estimator::Estimator) = estimator.internals.default_internals.index_set_size.max_sz

set_sz(estimator::Estimator, n) = begin
    if estimator isa Estimator{<:AD} && n > max_sz(estimator)
        update_boundary(estimator)
    end
    estimator.internals.default_internals.index_set_size.sz = n
    estimator.internals.default_internals.index_set_size.max_sz = max(sz(estimator), max_sz(estimator))
end