using Random, LinearAlgebra, Pickle;
using StatsPlots, ColorSchemes;
include("../binary_search.jl");
include("boundeddists.jl");

## Getting identification strategy

function everybody(expe, wstar)
    if expe == "all"
        iss = [(TopTwo(0.5, "EB", "TC"), GLRT("TT")),
               (TopTwo(0.5, "EB", "TCI-log"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TC"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TCI-log"), GLRT("TT")),
               (NPGame(CTracking, true), GLRT("TT")),
               (LUCB("KL"), GLRT("TT")),
               (FixedWeights(wstar), GLRT("TT")),
               (RoundRobin(), GLRT("TT"))];
    elseif expe == "all_GK16"
        iss = [(TopTwo(0.5, "EB", "TC"), GLRT("GK16")),
                (TopTwo(0.5, "EB", "RS"), GLRT("GK16")),
                (TopTwo(0.5, "EB", "TCI-log"), GLRT("GK16")),
                (TopTwo(0.5, "TS", "TC"), GLRT("GK16")),
                (TopTwo(0.5, "TS", "RS"), GLRT("GK16")),
                (TopTwo(0.5, "TS", "TCI-log"), GLRT("GK16")),
                (NPGame(CTracking, true), GLRT("GK16")),
                (LUCB("KL"), GLRT("GK16")),
                (FixedWeights(wstar), GLRT("GK16")),
                (RoundRobin(), GLRT("GK16"))];
    elseif expe == "random"
        iss = [(TopTwo(0.5, "EB", "TC"), GLRT("TT")),
               (TopTwo(0.5, "EB", "TCI-log"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TC"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TCI-log"), GLRT("TT")),
               (NPGame(CTracking, true), GLRT("TT")),
               (LUCB("KL"), GLRT("TT")),
               (RoundRobin(), GLRT("TT"))];
   elseif expe == "random_GK16"
       iss = [(TopTwo(0.5, "EB", "TC"), GLRT("GK16")),
              (TopTwo(0.5, "EB", "RS"), GLRT("GK16")),
              (TopTwo(0.5, "EB", "TCI-log"), GLRT("GK16")),
              (TopTwo(0.5, "TS", "TC"), GLRT("GK16")),
              (TopTwo(0.5, "TS", "RS"), GLRT("GK16")),
              (TopTwo(0.5, "TS", "TCI-log"), GLRT("GK16")),
              (NPGame(CTracking, true), GLRT("GK16")),
              (LUCB("KL"), GLRT("GK16")),
              (RoundRobin(), GLRT("GK16"))];
    elseif expe == "fail_equ"
        iss = [(TopTwo(0.5, "EB", "TC"), GLRT("TT")),
               (TopTwo(0.5, "EB", "TCI-log"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TC"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TCI-log"), GLRT("TT")),
               (NPGame(CTracking, true), GLRT("TT")),
               (LUCB("KL"), GLRT("TT")),
               (FixedWeights(wstar), GLRT("TT")),
               (RoundRobin(), GLRT("TT"))];
     elseif expe == "fail_equ_GK16"
         iss = [(TopTwo(0.5, "EB", "TC"), GLRT("GK16")),
                (TopTwo(0.5, "EB", "RS"), GLRT("GK16")),
                (TopTwo(0.5, "EB", "TCI-log"), GLRT("GK16")),
                (TopTwo(0.5, "TS", "TC"), GLRT("GK16")),
                (TopTwo(0.5, "TS", "RS"), GLRT("GK16")),
                (TopTwo(0.5, "TS", "TCI-log"), GLRT("GK16")),
                (NPGame(CTracking, true), GLRT("GK16")),
                (LUCB("KL"), GLRT("GK16")),
                (FixedWeights(wstar), GLRT("GK16")),
                (RoundRobin(), GLRT("GK16"))];
    elseif expe == "dsat_laz"
        iss = [(TopTwo(0.5, "EB", "TC"), LazyGLRT("TT", 10, 0.2)),
               (TopTwo(0.5, "EB", "TCI-log"), LazyGLRT("TT", 10, 0.2)),
               (TopTwo(0.5, "TS", "TC"), LazyGLRT("TT", 10, 0.2)),
               (TopTwo(0.5, "TS", "TCI-log"), LazyGLRT("TT", 10, 0.2)),
               (LUCB("Ki"), LazyGLRT("TT", 10, 0.2)),
               (FixedWeights(wstar), LazyGLRT("TT", 10, 0.2)),
               (RoundRobin(), LazyGLRT("TT", 10, 0.2))];
    elseif expe == "dsat_exp"
        iss = [(TopTwo(0.5, "EB", "TCI-log"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TC"), GLRT("TT")),
               (TopTwo(0.5, "TS", "TCI-log"), GLRT("TT")),
               (NPGame(CTracking, true), GLRT("TT")),
               (LUCB("Ki"), GLRT("TT")),
               (FixedWeights(wstar), GLRT("TT")),
               (RoundRobin(), GLRT("TT"))];
    elseif expe == "dsat_exp_GK16"
        iss = [(TopTwo(0.5, "EB", "TC"), GLRT("GK16")),
               (TopTwo(0.5, "EB", "RS"), GLRT("GK16")),
               (TopTwo(0.5, "EB", "TCI-log"), GLRT("GK16")),
               (TopTwo(0.5, "TS", "TC"), GLRT("GK16")),
               (TopTwo(0.5, "TS", "RS"), GLRT("GK16")),
               (TopTwo(0.5, "TS", "TCI-log"), GLRT("GK16")),
               (NPGame(CTracking, true), GLRT("GK16")),
               (LUCB("Ki"), GLRT("GK16")),
               (FixedWeights(wstar), GLRT("GK16")),
               (RoundRobin(), GLRT("GK16"))];
    elseif expe == "test"
        iss = [(TopTwo(0.5, "TS", "RS"), GLRT("GK16")),
               (TopTwo(0.5, "TS", "TC"), GLRT("GK16"))];
    else
        @error "Not implemented";
    end
    return iss;
end


## Instances

function get_instance_experiment(param_inst)
    if param_inst["inst"] == "ber_easy"
        μs = [0.7, 0.6, 0.5, 0.4, 0.3];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif param_inst["inst"] == "ber_hard"
        μs = [0.7, 0.65, 0.6, 0.55, 0.5];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif param_inst["inst"] == "ber_triv"
        μs = [0.7, 0.5, 0.4, 0.3, 0.2];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif param_inst["inst"] == "ber_eq3rd"
        μs = [0.7, 0.6, 0.5, 0.5];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif param_inst["inst"] == "fail_eq2nd_easy"
        μs = [0.5, 0.4, 0.4];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif param_inst["inst"] == "fail_eq2nd_okay"
        μs = [0.5, 0.425, 0.425];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif param_inst["inst"] == "fail_eq2nd_hard"
        μs = [0.5, 0.45, 0.45];
        dists = [BernoulliGen(μ, 1.) for μ in μs];
        return μs, dists, 1.;
    elseif occursin("dsat", inst)
        # Load
        dsat_file = param_inst["data_dir"] * "20210501_214815_dssat_samples_N_1000000_ST_6_D_15_SL_HC_GEN0027_E_1.pkl";
        dsat_data = Pickle.npyload(dsat_file);
        _samples = dsat_data["samples"];

        # We rescale the data
        scale = ceil(maximum(_samples)/100) * 100;
        if split(inst, "_")[2] == "7"
            samples = [_samples[a,:] / scale for a in 1:param_inst["nK"]];
        elseif split(inst, "_")[2] == "4"
            samples = [_samples[a,:] / scale for a in [3, 2, 5, 1]];
        elseif split(inst, "_")[2] == "3"
            samples = [_samples[a,:] / scale for a in [3, 5, 1]];
        elseif split(inst, "_")[2] == "5"
            samples = [_samples[a,:] / scale for a in [3, 6, 2, 5, 1]];
        elseif split(inst, "_")[2] == "2"
            samples = [_samples[a,:] / scale for a in [3, 1]];
        elseif split(inst, "_")[2] == "6"
            samples = [_samples[a,:] / scale for a in [3, 4, 6, 2, 5, 1]];
        else
            @error "Not Implemented";
        end
        μs = sum.(samples) / size(_samples)[2];
        dists = [HistogramGen(μs[a], samples[a], 1.) for a in 1:param_inst["nK"]];
        return μs, dists, 1.;
    else
        @error "Not Implemented";
    end
end

## Summary

function print_summary(pep, dists, μs, Tstar, wstar, δs, iss, datas, Nruns, file)
    nK = nanswers(pep);
    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_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];
        τs = [sum(x[2]) for x in data[r, :]];
        Eτ = mean(τs);
        στ = std(τs);
        err = sum(x->x[1] .!= 1, data[r, :])/Nruns;
        tim = sum(x->x[3], data[r, :])/Nruns;

        for io in ios
            println(io, @sprintf("%-30s", abbrev(is[1]) * "-" * abbrev(is[2])),
                        @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[1]) * "-" * abbrev(is[2]);
        end
    end
    for io in ios
        println(io, rule);
    end
end

## Plots

name_algos_plots = Dict("EB-TC-5-G-GK16" => "EB-TC",
                        "EB-RS-5-G-GK16" => "EB-RS",
                        "EB-TCI-log-5-G-GK16" => "EB-TCI",
                        "EB-TCI-sqrt-5-G-GK16" => "EB-TCIs",
                        "TS-TC-5-G-GK16" => "TS-TC",
                        "TS-RS-5-G-GK16" => "TS-RS",
                        "TS-TCI-log-5-G-GK16" => "TS-TCI",
                        "TS-TCI-sqrt-5-G-GK16" => "TS-TCIs",
                        "opt-G-GK16" => "Fixed",
                        "NPGame-C-G-GK16" => "WF-DKM",
                        "F-NPGame-C-G-GK16" => "Kinf-DKM",
                        #"F-NPGame-C-G-GK16" => "Ki-DKM",
                        "Ho-LUCB-G-GK16" => "Ho-LUCB",
                        "KL-LUCB-G-GK16" => "KL-LUCB",
                        "Ki-LUCB-G-GK16" => "Kinf-LUCB",
                        #"Ki-LUCB-G-GK16" => "Ki-LUCB",
                        "RR-G-GK16" => "Uniform",
                        "EB-TC-5-LG-GK16" => "EB-TC",
                        "EB-RS-5-LG-GK16" => "EB-RS",
                        "EB-TCI-log-5-LG-GK16" => "EB-TCI",
                        "EB-TCI-sqrt-5-LG-GK16" => "EB-TCIs",
                        "TS-TC-5-LG-GK16" => "TS-TC",
                        "TS-RS-5-LG-GK16" => "TS-RS",
                        "TS-TCI-log-5-LG-GK16" => "TS-TCI",
                        "TS-TCI-sqrt-5-LG-GK16" => "TS-TCIs",
                        "opt-LG-GK16" => "Fixed",
                        "NPGame-C-LG-GK16" => "WF-DKM",
                        "F-NPGame-C-LG-GK16" => "Kinf-DKM",
                        "Ho-LUCB-LG-GK16" => "Ho-LUCB",
                        "Ki-LUCB-LG-GK16" => "Kinf-LUCB",
                        "KL-LUCB-LG-GK16" => "KL-LUCB",
                        "RR-LG-GK16" => "Uniform",
                        "EB-TC-5-G-TT" => "EB-TC",
                        "EB-RS-5-G-TT" => "EB-RS",
                        "EB-TCI-log-5-G-TT" => "EB-TCI",
                        "EB-TCI-sqrt-5-G-TT" => "EB-TCIs",
                        "TS-TC-5-G-TT" => "TS-TC",
                        "TS-RS-5-G-TT" => "TS-RS",
                        "TS-TCI-log-5-G-TT" => "TS-TCI",
                        "TS-TCI-sqrt-5-G-TT" => "TS-TCIs",
                        "opt-G-TT" => "Fixed",
                        "NPGame-C-G-TT" => "WF-DKM",
                        "F-NPGame-C-G-TT" => "Kinf-DKM",
                        "Ho-LUCB-G-TT" => "Ho-LUCB",
                        "KL-LUCB-G-TT" => "KL-LUCB",
                        "Ki-LUCB-G-TT" => "Kinf-LUCB",
                        "RR-G-TT" => "Uniform",
                        "EB-TC-5-G-FTT" => "EB-TC",
                        "EB-RS-5-G-FTT" => "EB-RS",
                        "EB-TCI-log-5-G-FTT" => "EB-TCI",
                        "EB-TCI-sqrt-5-G-FTT" => "EB-TCIs",
                        "TS-TC-5-G-FTT" => "TS-TC",
                        "TS-RS-5-G-FTT" => "TS-RS",
                        "TS-TCI-log-5-G-FTT" => "TS-TCI",
                        "TS-TCI-sqrt-5-G-FTT" => "TS-TCIs",
                        "opt-G-FTT" => "Fixed",
                        "NPGame-C-G-FTT" => "WF-DKM",
                        "F-NPGame-C-G-FTT" => "Kinf-DKM",
                        "Ho-LUCB-G-FTT" => "Ho-LUCB",
                        "KL-LUCB-G-FTT" => "KL-LUCB",
                        "Ki-LUCB-G-FTT" => "Kinf-LUCB",
                        "RR-G-FTT" => "Uniform",
                        "EB-TC-5-LG-TT" => "EB-TC",
                        "EB-RS-5-LG-TT" => "EB-RS",
                        "EB-TCI-log-5-LG-TT" => "EB-TCI",
                        "EB-TCI-sqrt-5-LG-TT" => "EB-TCIs",
                        "TS-TC-5-LG-TT" => "TS-TC",
                        "TS-RS-5-LG-TT" => "TS-RS",
                        "TS-TCI-log-5-LG-TT" => "TS-TCI",
                        "TS-TCI-sqrt-5-LG-TT" => "TS-TCIs",
                        "opt-LG-TT" => "Fixed",
                        "NPGame-C-LG-TT" => "WF-DKM",
                        "F-NPGame-C-LG-TT" => "Kinf-DKM",
                        "Ho-LUCB-LG-TT" => "Ho-LUCB",
                        "Ki-LUCB-LG-TT" => "Kinf-LUCB",
                        "KL-LUCB-LG-TT" => "KL-LUCB",
                        "RR-LG-TT" => "Uniform");


pal_colors = palette(:tab10);
dict_colors = Dict("EB-TC" => pal_colors[2],
                   "EB-TCI" => pal_colors[1],
                   "EB-RS" => pal_colors[8],
                   "TS-TC" => pal_colors[3],
                   "TS-TCI" => pal_colors[4],
                   "TS-RS" => pal_colors[7],
                   "Kinf-DKM" => pal_colors[5],
                   "Ki-DKM" => pal_colors[5],
                   "Kinf-LUCB" => pal_colors[6],
                   "Ki-LUCB" => pal_colors[6],
                   "KL-LUCB" => pal_colors[6],
                   "Fixed" => pal_colors[10],
                   "Uniform" => pal_colors[9]);

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

function filter_plot_samp(is)
    sr, rsp = is;
    #if occursin("TCI-sqrt", abbrev(sr)) || occursin("EB-TC-", abbrev(sr)) || occursin("-RS-", abbrev(sr))
    if occursin("TCI-sqrt", abbrev(sr))
        return false;
    else
        return true;
    end
end

function _plot_samp(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];

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

    for r in eachindex(iss)
        τs = [sum(x[2]) for x in data[r, :]];
        x = name_algos_plots[abbrev(iss[r][1]) * "-" * abbrev(iss[r][2])];
        if r == 1
            boxplot(
                [x],
                [τs],
                label="",
                notch=true,
                xtickfontsize=8,
                ytickfontsize=12,
                #ylim=(0, 110000),
                #ylim=(0, 10 * lbd),
                ylim=(0, 25 * lbd),
                color=dict_colors[x],
                outliers=true);
        else
            boxplot!(
                [x],
                [τs],
                label="",
                notch=true,
                xtickfontsize=8,
                ytickfontsize=12,
                #ylim=(0, 110000),
                #ylim=(0, 10 * lbd),
                ylim=(0, 25 * lbd),
                color=dict_colors[x],
                outliers=true);
        end
        plot!([x], [mean(τs)], marker=(:star4,10,:black), label="");
    end

    #plot!([lbd], seriestype=:hline, legend=:topleft, label="lbd");
    plot!([lbd], seriestype=:hline, legend=:topright, label="lbd");
    #plot!([lbd], seriestype=:hline, legend=:topleft, label="");

    savefig(file);
end

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

function filter_plot_arms(is)
    sr, rsp = is;
    #if occursin("TCI-sqrt", abbrev(sr)) || occursin("EB-TC-", abbrev(sr)) || occursin("-RS-", abbrev(sr))
    if occursin("TCI-sqrt", abbrev(sr))
        return false;
    else
        return true;
    end
end

function _plot_arms(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];

    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 occursin("TCI-sqrt", abbrev(sr)) || occursin("EB-TC-", abbrev(sr)) || occursin("-RS-", abbrev(sr))
    if occursin("TCI-sqrt", abbrev(sr))
        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, Nruns, file)
    keep = filter(i -> filter_plot_rand_samp(iss[i]), 1:length(iss));
    _plot_rand_samp(iss[keep], data, iss_index, Nruns, file);
end

function filter_plot_rand_samp(is)
    sr, rsp = is;
    #if occursin("TCI-sqrt", abbrev(sr)) || occursin("EB-TC-", abbrev(sr))
    if occursin("TCI-sqrt", abbrev(sr))
        return false;
    else
        return true;
    end
end

function _plot_rand_samp(iss, data, iss_index, Nruns, file)
    med_τs = zeros(length(iss));
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        med_τs[i] = median([sum(x[2]) for x in data[r, :]]);
    end
    #med_τ = 8 * median(med_τs);
    med_τ = 9 * median(med_τs);
    for (i, is) in enumerate(iss)
        r = iss_index[is];
        τs = [sum(x[2]) for x in data[r, :]];
        x = name_algos_plots[abbrev(is[1]) * "-" * abbrev(is[2])];
        if i == 1
            boxplot(
                [x],
                [τs],
                label="",
                notch=true,
                ylim=(0, med_τ),
                xtickfontsize=8,
                ytickfontsize=12,
                color=dict_colors[x],
                outliers=true);
        else
            boxplot!(
                [x],
                [τs],
                label="",
                notch=true,
                ylim=(0, med_τ),
                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_histogram(dists, file)
    nK = length(dists);
    samples = Any[];
    for dist in dists
        if typeof(dist) == HistogramGen
            samples_a = dist.histogram;
        elseif typeof(dist) == BernoulliGen
            rng = MersenneTwister(42);
            samples_a = [sample(rng, dist::BernoulliGen) for i in 1:1000000]
        else
            @warn "Not implemented"
        end
        push!(samples, samples_a);
    end
    colors = palette(:tol_bright);
    density(samples[1], color=colors[1], linewidth=4, label="Arm 1", yticks=:none, legend=:topright, legendfontsize=12);
    vline!([dists[1].μ], color=colors[1], line=:dashdot, linewidth=2, label="");
    for a in 2:nK
        density!(samples[a], color=colors[a], linewidth=4, label="Arm $(a)");
        vline!([dists[a].μ], color=colors[a], line=:dashdot, linewidth=2, label="");
    end
    savefig(file);
end



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

function filter_plot_table_timeout(is)
    sr, rsp = is;
    if occursin("TCI-sqrt", abbrev(sr))
        return false;
    else
        return true;
    end
end

function _plot_table_timeout(pep, μs, Tstar, wstar, δ, iss, data, N, file)
    io = open(file, "a");
    rule = repeat("-", 90);
    println(io, "");
    println(io, rule);
    println(io, @sprintf("%10s", "δ"), "    ", @sprintf("%10.3f", δ));
    println(io, rule);

    for r in eachindex(iss)
        x = name_algos_plots[abbrev(iss[r][1]) * "-" * abbrev(iss[r][2])];

        err = sum(x->x[1] .!= 1, data[r, :]) / N;
        τs = [sum(x[2]) for x in data[r, :]];
        taumax = maximum(τs);
        timeout = sum(τs .== taumax) / N;
        println(io, x);
        println(io, err);
        println(io, timeout);
    end
    println(io, rule);
end
