using Random, LinearAlgebra, Pickle;
using Statistics, StatsPlots, ColorSchemes;
using LaTeXStrings;
import Distributions;
include("../binary_search.jl");
include("../expfam.jl");

## Getting identification strategy

function everybody(expe, wstar)
    if expe == "fast"
        iss = [(AdaTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCI", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIs", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIsp", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "RS", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "UCBI", "TC", "B"), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (LUCBhalf(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "random"
        iss = [(AdaTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCI", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIs+", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIs", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIsp+", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIsp", "B"), GLRT("GK16")),
               (AdaTopTwo("ada", "EB", "TCI", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "TC", "B"), GLRT("GK16")),
               (AdaTopTwo("ada", "TS", "TC", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "RS", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "RS", "B"), GLRT("GK16")),
               (AdaTopTwo("ada", "TS", "RS", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "UCBI", "TC", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "UCBI", "TC", ""), GLRT("GK16")),
               (AdaTopTwo("ada", "UCBI", "TC", "B"), GLRT("GK16")),
               (FWSampling(CTracking), GLRT("GK16")),
               (DKM(CTracking), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (LUCBhalf(), GLRT("GK16")),
               (TaS(CTracking), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "penaTCI"
        iss = [(AdaTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCI", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIs+", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIs", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIsp+", "B"), GLRT("GK16")),
               (AdaTopTwo("cst", "EB", "TCIsp", "B"), GLRT("GK16"))];
    elseif expe == "fastAdapt"
        iss = [(AdaTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaTopTwo("ada", "EB", "TCI", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaTopTwo("ada", "TS", "TC", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "TS", "RS", ""), GLRT("GK16")),
               (AdaTopTwo("ada", "TS", "RS", ""), GLRT("GK16")),
               (AdaTopTwo("cst", "UCBI", "TC", "B"), GLRT("GK16")),
               (AdaTopTwo("ada", "UCBI", "TC", "B"), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (LUCBhalf(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "testAdapt"
        iss = [(RoundRobin(), GLRT("GK16"))];
    elseif expe == "epsFull"
        iss = [(AnytimeTopTwo(true, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.05), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.5, "log"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.1, "log"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.05, "log"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.5, "pol"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.1, "pol"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.05, "pol"), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "UCB", "TC", "B"), GLRT("GK16")),
               (BAITopTwo(true, true, "TS", "TC"), GLRT("GK16")),
               (BAITopTwo(true, true, "EB", "TCI"), GLRT("GK16")),
               (BAITopTwo(true, false, "UCB", "TC"), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.05), GLRT("GK16")),
               (AdaEpsTopTwo("ada", "TS", "TC", ""), GLRT("GK16")),
               (AdaEpsTopTwo("ada", "EB", "TCI", ""), GLRT("GK16")),
               (AdaEpsTopTwo("ada", "UCB", "TC", "B"), GLRT("GK16")),
               (TaS(CTracking), GLRT("GK16")),
               (EpsFWSampling(CTracking), GLRT("GK16")),
               (DKM(CTracking), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "eps2G"
        iss = [(AnytimeTopTwo(true, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.05), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "UCB", "TC", "B"), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.05), GLRT("GK16")),
               (TaS(CTracking), GLRT("GK16")),
               (EpsFWSampling(CTracking), GLRT("GK16")),
               (DKM(CTracking), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "epsFast"
        iss = [(AnytimeTopTwo(true, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.05), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "UCB", "TC", "B"), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.05), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "epsPoE"
        iss = [(AnytimeTopTwo(true, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.05), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "TS", "TC", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "EB", "TCI", ""), GLRT("GK16")),
               (AdaEpsTopTwo("cst", "UCB", "TC", "B"), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(false, 0.05), GLRT("GK16")),
               (TaS(CTracking), GLRT("GK16")),
               (EpsFWSampling(CTracking), GLRT("GK16")),
               (DKM(CTracking), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (DSR(), GLRT("GK16")),
               (DSH(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "epsVSBench"
        iss = [(AnytimeTopTwo(true, 0.15), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.1), GLRT("GK16")),
               (AnytimeTopTwo(true, 0.05), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.5, "log"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.1, "log"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.05, "log"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.5, "pol"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.1, "pol"), GLRT("GK16")),
               (VSAnytimeTopTwo(true, 0.05, "pol"), GLRT("GK16")),
               (BAITopTwo(true, true, "TS", "TC"), GLRT("GK16")),
               (BAITopTwo(true, true, "EB", "TCI"), GLRT("GK16")),
               (BAITopTwo(true, false, "UCB", "TC"), GLRT("GK16")),
               (TaS(CTracking), GLRT("GK16")),
               (EpsFWSampling(CTracking), GLRT("GK16")),
               (DKM(CTracking), GLRT("GK16")),
               (LUCB(), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "epsTest"
        iss = [(RoundRobin(), GLRT("GK16"))];
    else
        @error "Not implemented";
    end
    return iss;
end


## Instances

function get_instance_experiment(param_inst)
    K = param_inst["nK"];
    if param_inst["inst"] == "sparse"
        μs = [(k == 1) ? 0.25 : 0. for k in 1:K];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "alpha3"
        μs = [1. - ((k - 1) / K)^(0.3) for k in 1:K];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "alpha6"
        μs = [1. - ((k - 1) / K)^(0.6) for k in 1:K];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "eqgap"
        μs = [(k == 1) ? 0 : -0.1 for k in 1:K];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "eqgapE"
        μs = [(k == 1) ? 0 : -0.5 for k in 1:K];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "trivial"
        μs = [(k == 1) ? 0 : -1 - k for k in 1:K];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "GK21_1"
        μs = [0.7, 0.55, 0.5, 0.4, 0.2];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "GK21_2"
        μs = [0.8, 0.75, 0.7, 0.6, 0.5, 0.4];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "GK21_3"
        μs = [0.6, 0.6, 0.55, 0.45, 0.3, 0.2];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "epsInst1"
        μs = [1, 0.97, 0.96, 0.9, 0.8];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "epsInst2"
        μs = [1, 0.97, 0.92, 0.9, 0.8];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "epsInst3"
        μs = [1, 0.96, 0.92, 0.9, 0.8];
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    elseif param_inst["inst"] == "2G"
        Ie = param_inst["Ie"];
        μs = 0.4 * ones(K);
        μs[1:Ie] .= 0.6;
        dists = [Gaussian() for μ in μs];
        return μs, dists;
    else
        @error "Not Implemented";
    end
end

## Summary

function print_summary(pep, dists, μs, Tstar, wstar, δs, iss, datas, Nruns, file)
    nK = nanswers(pep, μs);
    astar = istar(pep, μs);

    ios = (stdout, open(file, "a"));
    for i in 1:length(δs)
        δ = δs[i];
        data = getindex.(datas, i);

        # lower bound
        kl = (1 - 2 * δ) * log((1 - δ) / δ);
        lbd = kl * Tstar;

        rule = repeat("-", 90);
        for io in ios
            println(io, "");
            println(io, rule);
            println(io, long(pep));
            println(io, rule);
            println(io, @sprintf("%10s", "a*"), "  ",
                    @sprintf("%10s", "δ"), "  ",
                    @sprintf("%10s", "K"), "  ",
                    @sprintf("%10s", "LBD"));
            println(io, @sprintf("%10.3f", astar), "  ",
                    @sprintf("%10.3f", δ), "  ",
                    @sprintf("%10.0f", nK), "  ",
                    @sprintf("%10.0f", lbd));
            println(io, rule);
            println(io, @sprintf("%-30s", "Arm"),
                    join(map(k -> @sprintf("%6s", k), 1:nK)), " ",
                    @sprintf("%7s", "total"), "  ",
                    @sprintf("%7s", "std"), "  ",
                    @sprintf("%7s", "err"), "  ",
                    @sprintf("%7s", "time"));
            println(io, @sprintf("%-30s", "μs"),
                    join(map(x -> @sprintf("%6.2f", x), μs)));
            println(io, @sprintf("%-30s", "w*"),
                    join(map(x -> @sprintf("%6.2f", x), wstar)));
            println(io, rule);
            println(io, @sprintf("%-30s", "Sampling rule"));
        end

        # iss
        for r in eachindex(iss)
            τs = [sum(x[2]) for x in data[r,:]];
            Eτ = mean(τs);
            στ = std(τs);
            err = sum(x->x[1] .!= astar, data[r,:])/Nruns;
            tim = sum(x->x[3], data[r,:])/Nruns;

            for io in ios
                println(io, @sprintf("%-30s", abbrev(iss[r][1]) * " & " * abbrev(iss[r][2])),
                        join(map(k -> @sprintf("%6.0f", sum(x->x[2][k], data[r,:])/Nruns), 1:nK)), " ",
                        @sprintf("%7.0f", Eτ), "  ",
                        @sprintf("%7.0f", στ), "  ",
                        @sprintf("%7.5f", err), "  ",
                        @sprintf("%7.5f", tim/1e6)
                        );
            end
            if err > δ
                @warn "too many errors for $(iss[r])";
            end
        end

        for io in ios
            println(io, rule);
        end
    end
end

function print_eps_summary(pep, dists, μs, Tstar, wstar, δs, iss, datas, Nruns, file)
    nK = nanswers(pep, μs);
    astar = istar(pep, μs);
    ieps = iepss(pep, μs);

    ios = (stdout, open(file, "a"));
    for i in 1:length(δs)
        δ = δs[i];
        data = getindex.(datas, i);

        # lower bound
        kl = (1 - 2 * δ) * log((1 - δ) / δ);
        lbd = kl * Tstar;

        rule = repeat("-", 90);
        for io in ios
            println(io, "");
            println(io, rule);
            println(io, long(pep));
            println(io, rule);
            println(io, @sprintf("%10s", "a*"), "  ",
                     @sprintf("%10s", "ϵ"), "  ",
                    @sprintf("%10s", "δ"), "  ",
                    @sprintf("%10s", "K"), "  ",
                    @sprintf("%10s", "|I_ϵ|"), "  ",
                    @sprintf("%10s", "LBD"));
            println(io, @sprintf("%10.3f", astar), "  ",
                    @sprintf("%10.3f", pep.ϵ), "  ",
                    @sprintf("%10.3f", δ), "  ",
                    @sprintf("%10.0f", nK), "  ",
                    @sprintf("%10.0f", length(ieps)), "  ",
                    @sprintf("%10.0f", lbd));
            println(io, rule);
            println(io, @sprintf("%-30s", "Arm"),
                    join(map(k -> @sprintf("%6s", k), 1:nK)), " ",
                    @sprintf("%7s", "total"), "  ",
                    @sprintf("%7s", "std"), "  ",
                    @sprintf("%7s", "err"), "  ",
                    @sprintf("%7s", "err_istar"), "  ",
                    @sprintf("%7s", "time"));
            println(io, @sprintf("%-30s", "μs"),
                    join(map(x -> @sprintf("%6.2f", x), μs)));
            println(io, @sprintf("%-30s", "w*"),
                    join(map(x -> @sprintf("%6.2f", x), wstar)));
            println(io, rule);
            println(io, @sprintf("%-30s", "Sampling rule"));
        end

        # iss
        for r in eachindex(iss)
            τs = [sum(x[2]) for x in data[r,:]];
            Eτ = mean(τs);
            στ = std(τs);
            err = sum(x->!(x[1] in ieps), data[r,:])/Nruns;
            err_star = sum(x->x[1] .!= astar, data[r,:])/Nruns;
            tim = sum(x->x[3], data[r,:])/Nruns;

            for io in ios
                println(io, @sprintf("%-30s", abbrev(iss[r][1]) * " & " * abbrev(iss[r][2])),
                        join(map(k -> @sprintf("%6.0f", sum(x->x[2][k], data[r,:])/Nruns), 1:nK)), " ",
                        @sprintf("%7.0f", Eτ), "  ",
                        @sprintf("%7.0f", στ), "  ",
                        @sprintf("%7.5f", err), "  ",
                        @sprintf("%7.5f", err_star), "  ",
                        @sprintf("%7.5f", tim/1e6)
                        );
            end
            if err > δ
                @warn "too many errors for $(iss[r])";
            end
        end

        for io in ios
            println(io, rule);
        end
    end
end

function print_rand_summary(δs, iss, data, iss_index, param_inst, Nruns, file)
    ios = (stdout, open(file, "a"));
    δ = δs[1];
    K = param_inst["nK"];

    rule = repeat("-", 90);
    for io in ios
        println(io, "");
        println(io, rule);
        println(io, @sprintf("%10s", "δ"), "  ",
                    @sprintf("%10s", "K"), "  ",
                    @sprintf("%10s", "Δmin"), "  ",
                    @sprintf("%10s", "Δmax"), "  ",
                    @sprintf("%10s", "μ1"));
        println(io, @sprintf("%10.3f", δ), "  ",
                    @sprintf("%10.0f", K), "  ",
                    @sprintf("%10.0f", param_inst["gapmin"]), "  ",
                    @sprintf("%10.0f", param_inst["gapmax"]), "  ",
                    @sprintf("%10.0f", param_inst["mu1"]));
        println(io, rule);
        println(io, @sprintf("%-30s", ""),
                    @sprintf("%7s", "total"), "  ",
                    @sprintf("%7s", "std"), "  ",
                    @sprintf("%7s", "err"), "  ",
                    @sprintf("%7s", "time"));
        println(io, rule);
        println(io, @sprintf("%-30s", "Sampling rule"));
    end

    for is in iss
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        instances = getindex.(data[r, :], 2);
        τs = [sum(x[2]) for x in results];
        Eτ = mean(τs);
        στ = std(τs);
        err = sum(x->x[1] .!= 1, results)/Nruns;
        tim = sum(x->x[3], results)/Nruns;

        for io in ios
            println(io, @sprintf("%-30s", abbrev_is(is)),
                        @sprintf("%7.0f", Eτ), "  ",
                        @sprintf("%7.0f", στ), "  ",
                        @sprintf("%7.5f", err), "  ",
                        @sprintf("%7.5f", tim/1e6));
        end

        if err > δ
            @warn "too many errors for " * abbrev_is(is);
        end
    end
    for io in ios
        println(io, rule);
    end
end


function print_eps_rand_summary(δs, iss, data, iss_index, param_inst, Nruns, file)
    ios = (stdout, open(file, "a"));
    δ = δs[1];
    K = param_inst["nK"];
    ϵ = param_inst["eps"];
    rϵ = param_inst["reps"];
    opt = param_inst["opt"];
    Kϵ = Int64(ceil(K * rϵ));

    rule = repeat("-", 90);
    for io in ios
        println(io, "");
        println(io, rule);
        println(io, (opt == "mul" ? "Multiplicative" : "Additive") * " ϵ-BAI for Gaussian bandits");
        println(io, rule);
        println(io, @sprintf("%10s", "δ"), "  ",
                    @sprintf("%10s", "K"), "  ",
                    @sprintf("%10s", "ϵ"), "  ",
                    @sprintf("%10s", "|I_ϵ|"), "  ",
                    @sprintf("%10s", "μ1"));
        println(io, @sprintf("%10.3f", δ), "  ",
                    @sprintf("%10.0f", K), "  ",
                    @sprintf("%10.3f", ϵ), "  ",
                    @sprintf("%10.0f", Kϵ), "  ",
                    @sprintf("%10.0f", param_inst["mu1"]));
        println(io, rule);
        println(io, @sprintf("%-30s", ""),
                    @sprintf("%7s", "total"), "  ",
                    @sprintf("%7s", "std"), "  ",
                    @sprintf("%7s", "err"), "  ",
                    @sprintf("%7s", "err_istar"), "  ",
                    @sprintf("%7s", "time"));
        println(io, rule);
        println(io, @sprintf("%-30s", "Sampling rule"));
    end

    for is in iss
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        instances = getindex.(data[r, :], 2);
        τs = [sum(x[2]) for x in results];
        Eτ = mean(τs);
        στ = std(τs);
        err = sum(x->x[1] .> Kϵ, results)/Nruns;
        err_istar = sum(x->x[1] .!= 1, results)/Nruns;
        tim = sum(x->x[3], results)/Nruns;

        for io in ios
            println(io, @sprintf("%-30s", abbrev_is(is)),
                        @sprintf("%7.0f", Eτ), "  ",
                        @sprintf("%7.0f", στ), "  ",
                        @sprintf("%7.5f", err), "  ",
                        @sprintf("%7.5f", err_istar), "  ",
                        @sprintf("%7.5f", tim/1e6));
        end

        if err > δ
            @warn "too many errors for " * abbrev_is(is);
        end
    end
    for io in ios
        println(io, rule);
    end
end

abbrev_is(is) = abbrev(is[1]) * "-" * abbrev(is[2]);

## Plots

name_algos_plots = Dict("opt-G-GK16" => "Fixed",
                        "DKM-C-G-GK16" => "DKM",
                        "EB-TC-e15-G-GK16" => "EB-TC-e15",
                        "EB-TC-e1-G-GK16" => "EB-TC-e1",
                        "EB-TC-e05-G-GK16" => "EB-TC-e05",
                        "cst-EB-TC-e15-G-GK16" => "cst-EB-TC-e15",
                        "cst-EB-TC-e1-G-GK16" => "cst-EB-TC-e1",
                        "cst-EB-TC-e05-G-GK16" => "cst-EB-TC-e05",
                        "cst-UCB-TC-TB-G-GK16" => "TTUCB",
                        "cst-EB-TCI-G-GK16" => "EB-TCI",
                        "cst-UCB-TC-G-GK16" => "S-TTUCB",
                        "cst-UCBI-TC-G-GK16" => "S-TTUCB",
                        "cst-TS-TC-G-GK16" => "T3C",
                        "cst-TS-RS-G-GK16" => "TTTS",
                        "cst-IMED-TC-G-GK16" => "S-IMED-TC",
                        "cst-EB-TCI-T-G-GK16" => "T-EB-TCI",
                        "cst-UCB-TC-T-G-GK16" => "TTUCB",
                        "cst-UCBI-TC-T-G-GK16" => "TTUCB",
                        "cst-TS-TC-T-G-GK16" => "T-T3C",
                        "cst-TS-RS-T-G-GK16" => "T-TTTS",
                        "cst-IMED-TC-T-G-GK16" => "IMED-TC",
                        "ada-EB-TCI-G-GK16" => "A-EB-TCI",
                        "ada-UCB-TC-G-GK16" => "AS-TTUCB",
                        "ada-UCBI-TC-G-GK16" => "AS-TTUCB",
                        "ada-TS-TC-G-GK16" => "A-T3C",
                        "ada-TS-RS-G-GK16" => "A-TTTS",
                        "ada-IMED-TC-G-GK16" => "AS-IMED-TC",
                        "ada-EB-TCI-T-G-GK16" => "AT-EB-TCI",
                        "ada-UCB-TC-T-G-GK16" => "A-TTUCB",
                        "ada-UCBI-TC-T-G-GK16" => "A-TTUCB",
                        "ada-TS-TC-T-G-GK16" => "AT-T3C",
                        "ada-TS-RS-T-G-GK16" => "AT-TTTS",
                        "ada-IMED-TC-T-G-GK16" => "A-IMED-TC",
                        "FWS-C-G-GK16" => "FWS",
                        "EpsFWS-C-G-GK16" => "FWS",
                        "TaS-C-G-GK16" => "TaS",
                        "LUCB-G-GK16" => "LUCB",
                        "LUCBhalf-G-GK16" => "1/2-LUCB",
                        "RR-G-GK16" => "Uniform");

name_algos_plots_bis = Dict("opt-G-GK16" => "Fixed",
                            "DKM-C-G-GK16" => L"\varepsilon" * "-DKM",
                            #"DKM-C-G-GK16" => "DKM",
                            "cst-EB-TCI-G-GK16" => L"\varepsilon" * "-EB-TCI",
                            "cst-EB-BTCI-G-GK16" => "EB-TCI",
                            "cst-UCBI-TC-T-G-GK16" => L"\varepsilon" * "-TTUCB",
                            "cst-UCB-TC-TB-G-GK16" => L"\varepsilon" * "-TTUCB",
                            "cst-UCB-BTC-T-G-GK16" => "TTUCB",
                            "cst-TS-TC-G-GK16" => L"\varepsilon" * "-T3C",
                            "cst-TS-BTC-G-GK16" => "T3C",
                            "cst-TS-RS-G-GK16" => "TTTS",
                            "EB-TC-e15-G-GK16" => "EB-TC" * L"_{0.15}",
                            #"EB-TC-e1-G-GK16" => "EB-TC" * L"_{\varepsilon}",
                            "EB-TC-e1-G-GK16" => "EB-TC" * L"_{0.1}",
                            "EB-TC-e05-G-GK16" => "EB-TC" * L"_{0.05}",
                            "cst-EB-TC-e15-G-GK16" => "c-EB-TC" * L"_{0.15}",
                            #"cst-EB-TC-e1-G-GK16" => "c-EB-TC" * L"_{0.1}",
                            "cst-EB-TC-e1-G-GK16" => "EB-TC" * L"_{\varepsilon}",
                            "cst-EB-TC-e05-G-GK16" => "c-EB-TC" * L"_{0.05}",
                            "cst-EB-TC-La5-G-GK16" => "la5",
                            "cst-EB-TC-La1-G-GK16" => "la1",
                            "cst-EB-TC-La05-G-GK16" => "la05",
                            "cst-EB-TC-Pa5-G-GK16" => "pa5",
                            "cst-EB-TC-Pa1-G-GK16" => "pa1",
                            "cst-EB-TC-Pa05-G-GK16" => "pa05",
                            "FWS-C-G-GK16" => "FWS",
                            "EpsFWS-C-G-GK16" => L"\varepsilon" * "-FWS",
                            "TaS-C-G-GK16" => L"\varepsilon" * "-TaS",
                            #"TaS-C-G-GK16" => "TaS",
                            "DSH-G-GK16" => "DSH",
                            "DSR-G-GK16" => "DSR",
                            "LUCB-G-GK16" => "LUCB",
                            "LUCBhalf-G-GK16" => "hLUCB",
                            "RR-G-GK16" => "Uniform");


pal_colors = palette(:tab20);
blues_colors = palette(:Blues_9);
greens_colors = palette(:Greens_9);
set19_colors = palette(:Set1_9);
set28_colors = palette(:Set2_8);
cblind_colors = palette(:seaborn_colorblind6);
dict_colors = Dict("EB-TC-e15" => blues_colors[9],
                   "cst-EB-TC-e15" => blues_colors[8],
                   "EB-TC-e1" => blues_colors[7],
                   "cst-EB-TC-e1" => blues_colors[6],
                   "EB-TC-e05" => blues_colors[5],
                   "cst-EB-TC-e05" => blues_colors[4],
                   "EB-TCI" => set19_colors[3],
                   "T-EB-TCI" => pal_colors[2],
                   "A-EB-TCI" => pal_colors[2],
                   "AT-EB-TCI" => pal_colors[2],
                   "Original" => pal_colors[5],
                   "κ = 1.2" => pal_colors[1],
                   "κ = 2" => pal_colors[2],
                   "α = 2" => pal_colors[3],
                   "α = 1.2" => pal_colors[4],
                   "T3C" => set19_colors[1],
                   "T-T3C" => pal_colors[20],
                   "A-T3C" => pal_colors[20],
                   "AT-T3C" => pal_colors[20],
                   "TTTS" => pal_colors[17],
                   "T-TTTS" => pal_colors[18],
                   "A-TTTS" => pal_colors[18],
                   "AT-TTTS" => pal_colors[18],
                   "TTUCB" => set19_colors[4],
                   "S-TTUCB" => pal_colors[8],
                   "A-TTUCB" => pal_colors[8],
                   "AS-TTUCB" => pal_colors[8],
                   "IMED-TC" => pal_colors[10],
                   "IMED-TC-T" => pal_colors[11],
                   "ada-IMED-TC" => pal_colors[10],
                   "ada-IMED-TC-T" => pal_colors[11],
                   "TaS" => set19_colors[5],
                   "FWS" => set28_colors[2],
                   "EpsFWS" => pal_colors[15],
                   "DKM" => set19_colors[7],
                   "DSR" => set19_colors[3],
                   "DSH" => set19_colors[4],
                   "LUCB" => set19_colors[8],
                   "hLUCB" => pal_colors[17],
                   "Fixed" => set28_colors[8],
                   "Uniform" => set19_colors[9]);


dict_colors_Bench_All = Dict("EB-TC" * L"_{0.1}" => cblind_colors[3],
                             "EB-TC" * L"_{0.15}" => cblind_colors[1],
                             "c-EB-TC" * L"_{0.15}" => cblind_colors[2],
                             "EB-TC" * L"_{0.1}" => cblind_colors[3],
                             "c-EB-TC" * L"_{0.1}" => cblind_colors[4],
                             "EB-TC" * L"_{0.05}" => cblind_colors[5],
                             "c-EB-TC" * L"_{0.05}" => cblind_colors[6],
                             "EB-TC" * L"_{\varepsilon}" => cblind_colors[3],
                             L"\varepsilon" * "-T3C" => set19_colors[1],
                             L"\varepsilon" * "-EB-TCI" => set19_colors[3],
                             L"\varepsilon" * "-TTUCB" => set19_colors[4],
                             L"\varepsilon" * "-TaS" => cblind_colors[1],
                             L"\varepsilon" * "-FWS" => cblind_colors[5],
                             L"\varepsilon" * "-DKM" => cblind_colors[6],
                             "T3C" => set19_colors[1],
                             "EB-TCI" => set19_colors[3],
                             "TTUCB" => set19_colors[4],
                             "TaS" => cblind_colors[1],
                             "FWS" => cblind_colors[5],
                             "DKM" => cblind_colors[6],
                             "DSR" => set19_colors[3],
                             "DSH" => set19_colors[4],
                             "LUCB" => set19_colors[8],
                             "la5" => cblind_colors[3],
                             "la1" => cblind_colors[5],
                             "la05" => cblind_colors[6],
                             "pa5" => cblind_colors[4],
                             "pa1" => set19_colors[2],
                             "pa05" => cblind_colors[2],
                             "Uniform" => set19_colors[9]);

function plot_samp(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_samp(iss[i], iss_to_keep), 1:length(iss));
    _plot_samp(pep, μs, Tstar, wstar, δ, iss[keep], data[keep,:], N, file, iss_to_keep, name_iss);
end

function filter_plot_samp(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_samp(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    # lower bound
    kl = (1 - 2 * δ) * log((1 - δ) / δ);
    lbd = kl * Tstar;

    for r in eachindex(iss)
        _τs = [sum(x[2]) for x in data[r, :]];
        τs = [(τ <= 15 * lbd) ? τ : 15 * lbd for τ in _τs];
        xname = name_algos_plots_bis[abbrev(iss[r][1]) * "-" * abbrev(iss[r][2])];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];
        if r == 1
            boxplot(
                [xname],
                [τs],
                label="",
                #ylabel=L"E_{\nu}[\tau_{0,\delta}]",
                #ylabel=L"E_{\nu}[\tau_{0.15,0.01}]",
                ylabel=L"E_{\nu}[\tau_{\varepsilon,\delta}]",
                #ylim=(0, 20000),
                notch=true,
                xtickfontsize=8,
                ytickfontsize=12,
                color=xcolor,
                outliers=true);
        else
            boxplot!(
                [xname],
                [τs],
                label="",
                notch=true,
                xtickfontsize=8,
                ytickfontsize=12,
                color=xcolor,
                outliers=true);
        end
        plot!([xname], [mean(τs)], marker=(:star4,10,:black), label="");
    end

    #plot!([lbd], seriestype=:hline, legend=:topright, label=L"T_{0}(\mu) \log(1/\delta)");
    plot!([lbd], seriestype=:hline, legend=:topright, label=L"T_{\varepsilon}(\mu) \log(1/\delta)");
    #plot!([lbd], seriestype=:hline, legend=:topright, label=L"4.5 T_{0.15}(\mu)");

    savefig(file);
end


function plot_perror(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_perror(iss[i], iss_to_keep), 1:length(iss));
    _plot_perror(pep, μs, Tstar, wstar, δ, iss[keep], data[keep,:], N, file, iss_to_keep, name_iss);
end

function filter_plot_perror(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_perror(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    ieps = iepss(pep, μs);

    # Hard-coded default parameter
    freqHist = 50;
    history = "partial";

    med_list = Float64[];
    is_dict = Dict();
    for r in eachindex(iss)
        results = data[r, :];
        push!(med_list, median([sum(x[2]) for x in results]));
    end
    #budget = minimum(med_list);
    #budget = 2000;
    #budget = 1000;
    budget = 500;

    α = 0.05;
    k = Distributions.quantile(Distributions.Normal(), 1 - α / 2);
    start = 1;

    for r in eachindex(iss)
        xname = name_algos_plots_bis[abbrev(iss[r][1]) * "-" * abbrev(iss[r][2])];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];

        # Compute probability of errors: Mean and Wilson CI
        perrors = zeros(1, Int64(ceil(budget / freqHist)));
        real_Nruns = zeros(1, Int64(ceil(budget / freqHist)));
        times = [i*freqHist for i in eachindex(perrors)];
        times[end] = Int64(ceil(budget));

        results = data[r, :];
        for x in results
            if length(x[4]) > length(perrors)
                _recos = x[4][1:length(perrors)];
            else
                _recos = x[4];
            end
            if pep.opt == "add"
                perrors[1:length(_recos)] .+= maximum(μs) .- μs[_recos] .> pep.ϵ;
            else
                perrors[1:length(_recos)] .+= maximum(μs) .- μs[_recos] .> pep.ϵ * maximum(μs);
            end
            real_Nruns[1:length(_recos)] .+= 1;
        end

        p0_errors = (perrors .+ k^2) ./ (real_Nruns .+ k^2);
        perrors ./= real_Nruns;
        CIs_errors = k * sqrt.(real_Nruns .* perrors .* (1 .- perrors) .+ k^2 / 4) ./ (real_Nruns .+ k^2);

        if r == 1
            #plot(times[start:end], p0_errors[start:end], ribbon=CIs_errors[start:end], label="", linecolor = invisible(), fillalpha=.1, color=dict_colors[xname]);
            plot(times[start:end], perrors[start:end], label=xname, color=xcolor, xlabel=L"n", ylabel=L"P_{\nu}[\mu_{\hat \imath_{n}} < \mu^\star -" * string(pep.ϵ) * "]", linestyle=:auto, linewidth=3);
        else
            #plot!(times[start:end], p0_errors[start:end], ribbon=CIs_errors[start:end], label="", linecolor = invisible(), fillalpha=.1, color=dict_colors[xname]);
            plot!(times[start:end], perrors[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        end
    end

    savefig(file);
end


function plot_sregret(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_sregret(iss[i], iss_to_keep), 1:length(iss));
    _plot_sregret(pep, μs, Tstar, wstar, δ, iss[keep], data[keep,:], N, file, iss_to_keep, name_iss);
end

function filter_plot_sregret(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_sregret(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    # Hard-coded default parameter
    #freqHist = 50;
    #freqHist = 20;
    freqHist = 4;
    history = "partial";

    med_list = Float64[];
    is_dict = Dict();
    for r in eachindex(iss)
        results = data[r, :];
        push!(med_list, median([sum(x[2]) for x in results]));
    end
    #budget = minimum(med_list);
    #budget = 2000;
    budget = 1000;
    #budget = 500;
    start = 1;

    for r in eachindex(iss)
        xname = name_algos_plots_bis[abbrev(iss[r][1]) * "-" * abbrev(iss[r][2])];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];

        counts_reco = zeros(length(μs), Int64(ceil(budget / freqHist)));
        real_Nruns = zeros(1, Int64(ceil(budget / freqHist)));
        times = [i*freqHist for i in eachindex(real_Nruns)];
        times[end] = Int64(ceil(budget));

        results = data[r, :];
        for x in results
            if length(x[4]) > length(real_Nruns)
                _recos = x[4][1:length(real_Nruns)];
            else
                _recos = x[4];
            end
            for (t, a) in enumerate(_recos)
                counts_reco[a, t] += 1;
            end
            real_Nruns[1:length(_recos)] .+= 1;
        end

        means_reco = zeros(length(μs), Int64(ceil(budget / freqHist)));
        varperarm_reco = zeros(length(μs), Int64(ceil(budget / freqHist)));
        for a in 1:length(μs)
            means_reco[a,:] = counts_reco[a,:] * (maximum(μs) - μs[a]);
            varperarm_reco[a,:] = counts_reco[a,:] .* (1 .- counts_reco[a,:] ./ real_Nruns[1,:]) * (maximum(μs) - μs[a])^2;
        end
        covar_reco = zeros(1, Int64(ceil(budget / freqHist)));
        for a in 1:length(μs)
            for b in 1:length(μs)
                if b > a
                    covar_reco[1,:] += means_reco[a,:] .* means_reco[b,:];
                end
            end
        end
        sregrets = sum(means_reco, dims=1) ./ real_Nruns;
        std_sregrets = sqrt.(sum(varperarm_reco, dims=1) .- 2 * covar_reco ./ real_Nruns) ./ real_Nruns

        if r == 1
            plot(times[start:end], sregrets[start:end], ribbon=std_sregrets[start:end], label=xname, color=xcolor, linewidth=2, xlabel=L"n", ylabel=L"E_{\nu}[\mu^\star - \mu_{\hat \imath_{n}}]", linestyle=:auto, fillalpha=.1);
            #ylims!((0.,0.09));
            #plot(times[start:end], sregrets[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        else
            plot!(times[start:end], sregrets[start:end], ribbon=std_sregrets[start:end], label=xname, linestyle=:auto, fillalpha=.1, color=xcolor, linewidth=2);
            #plot!(times[start:end], sregrets[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        end
    end

    savefig(file);
end


function plot_cregret(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_cregret(iss[i], iss_to_keep), 1:length(iss));
    _plot_cregret(pep, μs, Tstar, wstar, δ, iss[keep], data[keep,:], N, file, iss_to_keep, name_iss);
end

function filter_plot_cregret(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_cregret(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    # Hard-coded default parameter
    freqHist = 50;
    history = "partial";

    med_list = Float64[];
    is_dict = Dict();
    for r in eachindex(iss)
        results = data[r, :];
        push!(med_list, median([sum(x[2]) for x in results]));
    end
    #budget = minimum(med_list);
    #budget = 2000;
    #budget = 1000;
    budget = 500;
    start = 1;

    for r in eachindex(iss)
        xname = name_algos_plots_bis[abbrev(iss[r][1]) * "-" * abbrev(iss[r][2])];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];

        sregrets = zeros(1, Int64(ceil(budget / freqHist)));
        squared_sregrets = zeros(1, Int64(ceil(budget / freqHist)));
        real_Nruns = zeros(1, Int64(ceil(budget / freqHist)));
        times = [i*freqHist for i in eachindex(sregrets)];
        times[end] = Int64(ceil(budget));

        results = data[r, :];
        for x in results
            if length(x[4]) > length(sregrets)
                _recos = x[4][1:length(sregrets)];
            else
                _recos = x[4];
            end
            sregrets[1:length(_recos)] .+= maximum(μs) .- μs[_recos];
            real_Nruns[1:length(_recos)] .+= 1;
        end
        sregrets ./= real_Nruns;
        cregrets = cumsum(sregrets[start:end])

        if r == 1
            #plot(times[start:end], sregrets[start:end], ribbon=std_regrets[start:end], label="", linestyle=:auto, fillalpha=.1, color=dict_colors[xname]);
            plot(times[start:end], cregrets, label=xname, color=xcolor, xlabel=L"T", ylabel=L"\sum_{n \le T} E_{\nu}[\mu^\star - \mu_{\hat \imath_{n}}]", linestyle=:auto, linewidth=3, legend=:topleft);
        else
            #plot!(times[start:end], sregrets[start:end], ribbon=std_regrets[start:end], label="", linestyle=:auto, fillalpha=.1, color=dict_colors[xname]);
            plot!(times[start:end], cregrets, label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        end
    end

    savefig(file);
end

function plot_arms(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_arms(iss[i], iss_to_keep), 1:length(iss));
    _plot_arms(pep, μs, Tstar, wstar, δ, iss[keep], data[keep,:], N, file, iss_to_keep, name_iss);
end

function filter_plot_arms(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_arms(pep, μs, Tstar, wstar, δ, iss, data, N, file, iss_to_keep, name_iss)
    # xlabel for identification strategies
    _xs = permutedims(vcat(collect((abbrev(sr) * "-" * abbrev(rsp) for (sr, rsp) in iss))...));
    xs = [name_algos_plots[x] for x in _xs];

    K = length(μs);

    ps = Any[];
    for i in 1:K
        # Mean empirical sampling proportions for arm i
        means = sum(getindex.(getindex.(data,2), i) ./ sum.(getindex.(data,2)), dims=2)/N;

        p = boxplot(
            xs,
            map(x -> x[2][i] / sum(x[2]), data)',  # the ' calls the adjoint
            label="",
            notch=true,
            outliers=true);

        # plot means
        plot!(p, xs, means', marker=(:star4,10,:black), label="", ylabel="w_" * string(i));

        push!(ps, p);
    end

    plot(ps..., layout=K);
    savefig(file);
end

function plot_cput(pep, μs, Tstar, wstar, δ, iss, data, N, file)
    keep = filter(i -> filter_plot_cput(iss[i]), 1:length(iss));
    _plot_cput(pep, μs, Tstar, wstar, δ, iss[keep], data[keep,:], N, file);
end

function filter_plot_cput(is)
    sr, rsp = is;
    if false
        return false;
    else
        return true;
    end
end

function _plot_cput(pep, μs, Tstar, wstar, δ, iss, data, N, file)
    # xlabel for identification strategies
    _xs = permutedims(vcat(collect((abbrev(sr) * "-" * abbrev(rsp) for (sr, rsp) in iss))...));
    xs = [name_algos_plots[x] for x in _xs];

    # Median CPU time
    medians = median(getindex.(data, 3), dims=2)/1e6;

    boxplot(
        xs,
        map(x -> x[3]/1e6, data)',  # the ' calls the adjoint
        label=:none,
        yaxis=:log,
        notch=true,
        outliers=true);

    savefig(file);
end

function plot_rand_samp(iss, data, iss_index, eps, Nruns, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_rand_samp(iss[i], iss_to_keep), 1:length(iss));
    _plot_rand_samp(iss[keep], data, iss_index, eps, Nruns, file, iss_to_keep, name_iss);
end

function filter_plot_rand_samp(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_rand_samp(iss, data, iss_index, eps, Nruns, file, iss_to_keep, name_iss)
    med_τs = zeros(length(iss));
    is_dict = Dict();
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        med_τs[i] = median([sum(x[2]) for x in results]);
        is_dict[abbrev(is[1]) * "-" * abbrev(is[2])] = is;
    end
    #med_τ = 8 * median(med_τs);
    #med_τ = 8000;
    #med_τ = 15000;
    #med_τ = 20000;
    #med_τ = 35000;
    for (i, is_idx) in enumerate(iss_to_keep)
        is = is_dict[is_idx];
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        τs = [sum(x[2]) for x in results];
        xname = name_iss[i];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];
        if i == 1
            boxplot(
                [xname],
                [τs],
                label="",
                notch=true,
                #ylim=(0, med_τ),
                #ylabel=L"E_{\nu}[\tau_{0.05,0.01}]",
                #ylabel=L"E_{\nu}[\tau_{0,\delta}]",
                ylabel=L"E_{\nu}[\tau_{\varepsilon,\delta}]",
                #ylim=(0, 35000),
                xtickfontsize=8,
                ytickfontsize=12,
                color=xcolor,
                outliers=true);
        else
            boxplot!(
                [xname],
                [τs],
                label="",
                notch=true,
                #ylim=(0, med_τ),
                #ylim=(0, 35000),
                xtickfontsize=8,
                ytickfontsize=12,
                color=xcolor,
                outliers=true);
        end
        plot!([xname], [mean(τs)], marker=(:star4,10,:black), label="");
    end

    savefig(file);
end

function plot_rand_perror(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_rand_perror(iss[i], iss_to_keep), 1:length(iss));
    _plot_rand_perror(iss[keep], data, iss_index, Nruns, file, iss_to_keep, name_iss);
end

function filter_plot_rand_perror(is, iss_to_keep)
    if abbrev(is[1]) * "-" * abbrev(is[2]) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_rand_perror(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    # Hard-coded default parameter
    freqHist = 50;
    history = "partial";

    med_list = Float64[];
    is_dict = Dict();
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        push!(med_list, median([sum(x[2]) for x in results]));
        is_dict[abbrev(is[1]) * "-" * abbrev(is[2])] = is;
    end
    #budget = minimum(med_list);
    #budget = 10000;
    #budget = 5000;
    budget = 2000;

    α = 0.05;
    k = Distributions.quantile(Distributions.Normal(), 1 - α / 2);
    start = 1;

    for (i, is_idx) in enumerate(iss_to_keep)
        is = is_dict[is_idx];
        r = iss_index[is];
        xname = name_iss[i];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];
        results = getindex.(data[r, :], 1);
        instances = getindex.(data[r, :], 2);

        # Compute probability of errors: Mean and Wilson CI
        perrors = zeros(1, Int64(ceil(budget / freqHist)));
        real_Nruns = zeros(1, Int64(ceil(budget / freqHist)));
        times = [i*freqHist for i in eachindex(perrors)];
        times[end] = Int64(ceil(budget));

        for (x, inst) in zip(results, instances)
            μs = inst[1];
            eps = inst[2];
            opt = inst[3];
            if length(x[4]) > length(perrors)
                _recos = x[4][1:length(perrors)];
            else
                _recos = x[4];
            end
            if opt == "add"
                perrors[1:length(_recos)] .+= maximum(μs) .- μs[_recos] .> eps;
            else
                perrors[1:length(_recos)] .+= maximum(μs) .- μs[_recos] .> eps * maximum(μs);
            end
            real_Nruns[1:length(_recos)] .+= 1;
        end

        p0_errors = (perrors .+ k^2) ./ (real_Nruns .+ k^2);
        perrors ./= real_Nruns;
        CIs_errors = k * sqrt.(real_Nruns .* perrors .* (1 .- perrors) .+ k^2 / 4) ./ (real_Nruns .+ k^2);

        if i == 1
            #plot(times[start:end], p0_errors[start:end], ribbon=CIs_errors[start:end], label="", linecolor = invisible(), fillalpha=.1, color=dict_colors[xname]);
            plot(times[start:end], perrors[start:end], label=xname, color=xcolor, xlabel=L"n", ylabel=L"P_{\nu}[\mu_{\hat \imath_{n}} < \mu^\star - $(eps)]", linestyle=:auto, linewidth=3);
        else
            #plot!(times[start:end], p0_errors[start:end], ribbon=CIs_errors[start:end], label="", linecolor = invisible(), fillalpha=.1, color=dict_colors[xname]);
            plot!(times[start:end], perrors[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        end
    end
    savefig(file);
end


function plot_rand_sregret(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_rand_sregret(iss[i], iss_to_keep), 1:length(iss));
    _plot_rand_sregret(iss[keep], data, iss_index, Nruns, file, iss_to_keep, name_iss);
end

function filter_plot_rand_sregret(is, iss_to_keep)
    if abbrev(is[1]) * "-" * abbrev(is[2]) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_rand_sregret(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    # Hard-coded default parameter
    freqHist = 50;
    history = "partial";

    med_list = Float64[];
    is_dict = Dict();
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        push!(med_list, median([sum(x[2]) for x in results]));
        is_dict[abbrev(is[1]) * "-" * abbrev(is[2])] = is;
    end
    #budget = minimum(med_list);
    #budget = 10000;
    #budget = 5000;
    budget = 1000;
    start = 1;

    for (i, is_idx) in enumerate(iss_to_keep)
        is = is_dict[is_idx];
        r = iss_index[is];
        xname = name_iss[i];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];
        results = getindex.(data[r, :], 1);
        instances = getindex.(data[r, :], 2);

        sregrets = zeros(1, Int64(ceil(budget / freqHist)));
        squared_sregrets = zeros(1, Int64(ceil(budget / freqHist)));
        real_Nruns = zeros(1, Int64(ceil(budget / freqHist)));
        times = [i*freqHist for i in eachindex(sregrets)];
        times[end] = Int64(ceil(budget));

        for (x, inst) in zip(results, instances)
            μs = inst[1];
            if length(x[4]) > length(sregrets)
                _recos = x[4][1:length(sregrets)];
            else
                _recos = x[4];
            end
            sregrets[1:length(_recos)] .+= maximum(μs) .- μs[_recos];
            squared_sregrets[1:length(_recos)] .+= (maximum(μs) .- μs[_recos]).^2;
            real_Nruns[1:length(_recos)] .+= 1;
        end
        sregrets ./= real_Nruns;
        std_sregrets = sqrt.(squared_sregrets ./ real_Nruns .- (sregrets ./ real_Nruns).^2);

        if i == 1
            plot(times[start:end], sregrets[start:end], ribbon=std_sregrets[start:end], label="", xlabel=L"n", ylabel=L"E_{\nu}[\mu^\star - \mu_{\hat \imath_{n}}]", linestyle=:auto, fillalpha=.1, color=xcolor);
            #plot(times[start:end], sregrets[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        else
            plot!(times[start:end], sregrets[start:end], ribbon=std_sregrets[start:end], label="", linestyle=:auto, fillalpha=.1, color=xcolor);
            #plot!(times[start:end], sregrets[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        end
    end
    savefig(file);
end


function plot_rand_cregret(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_rand_cregret(iss[i], iss_to_keep), 1:length(iss));
    _plot_rand_cregret(iss[keep], data, iss_index, Nruns, file, iss_to_keep, name_iss);
end

function filter_plot_rand_cregret(is, iss_to_keep)
    if abbrev(is[1]) * "-" * abbrev(is[2]) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_rand_cregret(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    # Hard-coded default parameter
    freqHist = 50;
    history = "partial";

    med_list = Float64[];
    is_dict = Dict();
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        push!(med_list, median([sum(x[2]) for x in results]));
        is_dict[abbrev(is[1]) * "-" * abbrev(is[2])] = is;
    end
    #budget = minimum(med_list);
    #budget = 10000;
    #budget = 5000;
    budget = 2000;
    start = 1;

    for (i, is_idx) in enumerate(iss_to_keep)
        is = is_dict[is_idx];
        r = iss_index[is];
        xname = name_iss[i];
        #xcolor = dict_colors[xname];
        xcolor = dict_colors_Bench_All[xname];
        results = getindex.(data[r, :], 1);
        instances = getindex.(data[r, :], 2);

        cregrets = zeros(1, Int64(ceil(budget / freqHist)));
        squared_cregrets = zeros(1, Int64(ceil(budget / freqHist)));
        real_Nruns = zeros(1, Int64(ceil(budget / freqHist)));
        times = [i*freqHist for i in eachindex(cregrets)];
        times[end] = Int64(ceil(budget));

        for (x, inst) in zip(results, instances)
            μs = inst[1];
            if length(x[4]) > length(cregrets)
                _recos = x[4][1:length(cregrets)];
            else
                _recos = x[4];
            end
            #cregret = cumsum(maximum(μs) .- μs[_recos]) ./ collect(1:length(_recos));
            cregret = cumsum(maximum(μs) .- μs[_recos]);
            cregrets[1:length(_recos)] .+= cregret;
            squared_cregrets[1:length(_recos)] .+= cregret.^2;
            real_Nruns[1:length(_recos)] .+= 1;
        end
        cregrets ./= real_Nruns;
        std_cregrets = sqrt.(squared_cregrets ./ real_Nruns .- cregrets.^2);

        if i == 1
            #plot(times[start:end], cregrets[start:end], ribbon=std_cregrets[start:end], label=xname, linestyle=:auto, fillalpha=.1, color=xcolor);
            plot(times[start:end], cregrets[start:end], label=xname, color=xcolor, xlabel=L"T", ylabel=L"\sum_{n \le T} E_{\nu}[\mu^\star - \mu_{\hat \imath_{n}}]", linestyle=:auto, legend=:topleft, linewidth=3);
        else
            #plot!(times[start:end], cregrets[start:end], ribbon=std_cregrets[start:end], label=xname, linestyle=:auto, fillalpha=.1, color=xcolor);
            plot!(times[start:end], cregrets[start:end], label=xname, color=xcolor, linestyle=:auto, linewidth=3);
        end
    end
    savefig(file);
end


function plot_rand_cput(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    keep = filter(i -> filter_plot_rand_cput(iss[i], iss_to_keep), 1:length(iss));
    _plot_rand_cput(iss[keep], data, iss_index, Nruns, file, iss_to_keep, name_iss);
end

function filter_plot_rand_cput(is, iss_to_keep)
    sr, rsp = is;
    if abbrev(sr) * "-" * abbrev(rsp) in iss_to_keep
        return true;
    else
        return false;
    end
end

function _plot_rand_cput(iss, data, iss_index, Nruns, file, iss_to_keep, name_iss)
    med_τs = zeros(length(iss));
    is_dict = Dict();
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        med_τs[i] = median([sum(x[3]) for x in results]);
        is_dict[abbrev(is[1]) * "-" * abbrev(is[2])] = is;
    end
    med_τ = 8 * median(med_τs);
    for (i, is_idx) in enumerate(iss_to_keep)
        is = is_dict[is_idx];
        r = iss_index[is];
        results = getindex.(data[r, :], 1);
        τs = [sum(x[3]) for x in results];
        x = name_iss[i];
        println("Algo: ", x);
        println("Mean cput: ", mean(τs)/1e6);
        println("Std cput: ", std(τs)/1e6);
        println("----")
        if i == 1
            boxplot(
                [x],
                [τs],
                label="",
                notch=true,
                #ylim=(0, med_τ),
                #ylim=(0, 35000),
                xtickfontsize=8,
                ytickfontsize=12,
                color=dict_colors[x],
                outliers=true);
        else
            boxplot!(
                [x],
                [τs],
                label="",
                notch=true,
                #ylim=(0, med_τ),
                #ylim=(0, 35000),
                xtickfontsize=8,
                ytickfontsize=12,
                color=dict_colors[x],
                outliers=true);
        end
        plot!([x], [mean(τs)], marker=(:star4,10,:black), label="");
    end

    savefig(file);
end


function plot_instances_dists(data, iss, iss_index, file_plot_charactimes, file_plot_allocbest)
    # Get the instances for the first iss since they are all identical
    r = iss_index[iss[1]];
    instances = getindex.(data[r, :], 2);

    Tstars = getindex.(instances, 2);
    Tstar_betas = getindex.(instances, 4);

    exp_ratios_time = exp.(Tstar_betas ./ Tstars);
    histogram(exp_ratios_time, bins=:scott);
    savefig(file_plot_charactimes);

    wstars = getindex.(instances, 3);
    wstars_astar = getindex.(wstars, 1);
    histogram(wstars_astar, bins=:scott);
    savefig(file_plot_allocbest);
end


