MasterMind/include/scripts.h

#ifndef SCRIPTS_H_INCLUDED
#define SCRIPTS_H_INCLUDED

#include <vector>

#include "analysis.h"
#include "ParameterOptimizer.h"

using namespace std;

class ProbaOptimizer : public ParameterOptimizer {
    public:
        ProbaOptimizer(int initial_n, std::default_random_engine& randomizer) : ParameterOptimizer(randomizer), m_initial_n{initial_n} {};
        SearchHeuristic* get_context(unsigned int i, double parameter, std::default_random_engine& randomizer) {
            EA* algo = new EA(1, 1, PLUS, parameter, randomizer);
            MasterMind* mm = new MasterMind(m_initial_n+i, m_initial_n+i, randomizer);
            algo->set_problem(mm);
            return algo;
        }
        virtual std::string display_context(unsigned int context) {
           std::stringstream ss;
           ss << "n=" << context+m_initial_n;
           return ss.str();
        }
        std::string get_context_for_csvfile(unsigned int num_context) {
           return std::to_string(m_initial_n+num_context);
        }

    private:
        int m_initial_n;
};

void optimal_proba_vmath3(int n_initial, int n_final, std::default_random_engine& randomizer, double p_init=0.05, double p_final=0.5, double p_step=0.01, vector<int> nb_tests_by_depth={}, vector<double> stq_by_depth={}) {
    ProbaOptimizer po(n_initial, randomizer);
    po.set_nb_contexts(n_final-n_initial);
    po.set_study_interval(p_init, p_final, p_step);
    for(int depth=0; depth<nb_tests_by_depth.size(); depth++)
        po.add_test_wave(nb_tests_by_depth[depth], stq_by_depth[depth], depth);
    po.run();
}

void script_optimal_proba(std::default_random_engine& randomizer) {
    std::string namefile = "optimal-proba.csv";
    {
        int n_init=3, n_final=16;
        double p_init=0.05, p_final=0.5, p_step=0.01;
        ProbaOptimizer po(n_init, randomizer);
        po.set_nb_contexts(n_final-n_init);
        po.set_study_interval(p_init, p_final, p_step);
        po.add_test_wave(20, 2.3362421597, 0);
        po.add_test_wave(200, 2.3451370823, 1);
        po.add_test_wave(2000, 2.2953981367, 2);
        po.save_run(namefile);
        po.run();
    }
    {
        int n_init=16, n_final=28;
        double p_init=0.03, p_final=0.2, p_step=0.01;
        ProbaOptimizer po(n_init, randomizer);
        po.set_nb_contexts(n_final-n_init);
        po.set_study_interval(p_init, p_final, p_step);
        po.add_test_wave(20, 2.6742086844, 0);
        po.add_test_wave(100, 2.4755765532, 1);
        po.add_test_wave(1000, 2.4365960683, 2);
        po.save_run(namefile);
        po.run();
    }
    {
        int n_init=28, n_final=46;
        double p_init=0.01, p_final=0.1, p_step=0.005;
        ProbaOptimizer po(n_init, randomizer);
        po.set_nb_contexts(n_final-n_init);
        po.set_study_interval(p_init, p_final, p_step);
        po.add_test_wave(20, 2.6742086844, 0);
        po.add_test_wave(100, 2.4755765532, 1);
        po.add_test_wave(1000, 2.4365960683, 2);
        po.save_run(namefile);
        po.run();
    }
}

vector<ConfidenceInterval> optimal_proba_vmath2(int n_initial, int n_final, std::default_random_engine& randomizer, double p_init=0.05, double p_final=0.5, double p_step=0.01, int nb_tests=10, double student_law_quartile=-1, bool displaying=true) {
    int nb_diff_proba = floor((p_final-p_init)/p_step);

    // Initialize result file
    ofstream probafile("probas-vmath.csv");
    probafile << "n;";
    for(int i=0; i<nb_diff_proba; i++)
        probafile << p_init+i*p_step << ";";
    probafile << endl;

    vector<ConfidenceInterval> res;

    for(int n=n_initial; n<n_final; n++) {
        auto start_n = start_chrono();

        Statistic stats[nb_diff_proba];
        double thresold = 0;
        // Realize the experiences
        for(int ind=0; ind<nb_diff_proba; ind++) {
            double p = p_init+ind*p_step;

            EA algo(1, 1, PLUS, p, randomizer);
            MasterMind mm(n, n, randomizer);
            algo.set_problem(&mm);
            stats[ind] = avg_calls(&algo, nb_tests, -1, student_law_quartile);
            double potential_thresold = stats[ind].average + stats[ind].alea;
            if(ind == 0 || thresold > potential_thresold)
                thresold = potential_thresold;

            auto stop_loop = stop_chrono();
            int elapsed_milliseconds_loop = interval_chrono(start_n, stop_loop);
            if(displaying) cout << '\r' << "n=" << n << ", p=" << p << "... executed in " << elapsed_milliseconds_loop/1000. << " s";
        }

        // Calculate the confidence interval of 95%
        if(displaying) cout << '\r' << "n=" << n << ". Thresold=" << thresold << ". Keep best... ";
        int min_possible = -1;
        int max_possible = -1;
        for(int i=0; i<nb_diff_proba; i++) {
            double min_bound = stats[i].average - stats[i].alea;

            if(min_possible < 0 && min_bound < thresold)
                min_possible = i;
            if(min_bound < thresold)
                max_possible = i;
        }
        double pmin = p_init+min_possible*p_step;
        double pmax = p_init+max_possible*p_step;

        // Store the results in a file
        probafile << n << ";";
        for(int i=0; i<nb_diff_proba; i++) {
            probafile << "("<<pmin<<","<<pmax<<")" << ";";
        }
        probafile << endl;

        auto stop_n = stop_chrono();
        int elapsed_milliseconds_n = interval_chrono(start_n, stop_n);
        if(displaying) cout << "Interval=[" << pmin << ", " << pmax << "]" << endl;
        if(displaying) cout << "... executed in " << elapsed_milliseconds_n/1000. << " seconds ";
        if(displaying) cout << "with " << nb_tests*nb_diff_proba << " resolutions" << endl;

        ConfidenceInterval ci(pmin, -1, pmax);
        res.push_back(ci);
    }

    return res;
}

vector<ConfidenceInterval> optimal_proba_vmath2_rec(int n_initial, int n_final, std::default_random_engine& randomizer, double base_p_init=0.05, double base_p_final=0.5, double p_step=0.01, vector<int> nb_tests_by_depth={}, vector<double> stq_by_depth={}) {
    vector<ConfidenceInterval> res;
    for(int n=n_initial; n<n_final; n++) {
        cout << "STUDY N=" << n << ":" << endl;
        vector<ConfidenceInterval> cis;
        int nb_tests = 0;
        for(unsigned int depth=0; depth<nb_tests_by_depth.size(); depth++) {
            double p_init = (depth == 0) ? base_p_init : cis[0].down-p_step/pow(2, depth-1);
            double p_final = (depth == 0) ? base_p_final : cis[0].up+p_step/depth;

            cout << " - depth=" << depth+1 << ": ";
            cis = optimal_proba_vmath2(n, n+1, randomizer, p_init, p_final, p_step/pow(2, depth), nb_tests_by_depth[depth], stq_by_depth[depth], false);
            nb_tests += nb_tests_by_depth[depth];
            cout << "["<< cis[0].down <<","<< cis[0].up <<"]" << endl;
        }
        res.push_back(cis[0]);

        ofstream probafile("probas-vmath-rec.csv", ios::app);
        probafile << n << ";" << cis[0].down << ";" << cis[0].up << ";" << 1./cis[0].up << ";" << 1./cis[0].down << ";" << nb_tests << endl;
    }

    return res;
}

double calculate_up_bound(double x, int nb) {
    return (x + 2./nb + sqrt(4./nb * (x - pow(x, 2) + 1./nb))) / (1. + 4./nb);
}
double calculate_down_bound(double x, int nb) {
    return (x + 2./nb - sqrt(4./nb * (x - pow(x, 2) + 1./nb))) / (1. + 4./nb);
}

void optimal_proba_vmath(int n_initial, int n_final, std::default_random_engine& randomizer, double p_init=0.05, double p_final=0.5, double p_step=0.01, int nb_loops=1000) {
    int nb_diff_proba = floor((p_final-p_init)/p_step);
    int nb_tests = 10;

    // Initialize result file
    ofstream probafile("probas-vmath.csv");
    probafile << "n;";
    for(int i=0; i<nb_diff_proba; i++)
        probafile << p_init+i*p_step << ";";
    probafile << endl;

    for(int n=n_initial; n<n_final; n++) {
        auto start_n = start_chrono();

        int nb_best[nb_diff_proba] = {0};
        int nb_solved_problems = 0;

        // Realize the experiences
        int old_min_ind = 0;
        int max_min_ind = 0;
        for(int num_loop=1; num_loop<nb_loops; num_loop++) {
            double min_avg = -1;
            int min_ind;
            // Realize the experience once
            for(int ind2=-1; ind2<nb_diff_proba; ind2++) {
                if(ind2 == old_min_ind) continue;
                int ind = (ind2 >= 0) ? ind2 : old_min_ind;
                double p = p_init+ind*p_step;

                EA algo(1, 1, PLUS, p, randomizer);
                MasterMind mm(n, n, randomizer);
                algo.set_problem(&mm);
                Statistic stat = avg_calls(&algo, nb_tests, min_avg+1);
                nb_solved_problems += stat.nb;

                double avg_nb_calls = stat.average;
                if(min_avg < 0 || min_avg > avg_nb_calls) {
                    min_avg = avg_nb_calls;
                    min_ind = ind;
                }
            }

            // Store the result of the experience
            nb_best[min_ind]++;
            if(nb_best[min_ind] > nb_best[max_min_ind])
                max_min_ind = min_ind;
            old_min_ind = min_ind;

            // Display temporary results
            auto stop_loop = stop_chrono();
            int elapsed_milliseconds_loop = interval_chrono(start_n, stop_loop);
            cout << "\r";
            cout << "n=" << n << ", num_loop=" << num_loop << ": ";
            cout << p_init+max_min_ind*p_step << " (" << nb_best[max_min_ind] << "), ";
            cout << "executed in " <<elapsed_milliseconds_loop/1000. << " seconds [" << nb_solved_problems << "]";
            cout << "    ";
        }

        // Display the global distribution
        cout << '\r' << "===== Graph ====                     " << endl;
        for(int j=9; j>=0; j--) {
            for(int i=0; i<nb_diff_proba; i++) {
                double pourcent = (double) nb_best[i]/nb_best[max_min_ind];
                cout << ((pourcent >= j/10.) ? '#' : ' ');
            }
            cout << endl;
        }
        cout << endl;

        // Calculate the confidence interval of 95%
        double pbest = p_init+max_min_ind*p_step;
        double x = (double) nb_best[max_min_ind]/nb_loops;
        cout << "x=" << x << ',' << nb_loops << ',' << nb_best[max_min_ind] << endl;
        double up_bound = calculate_up_bound(x, nb_loops);
        double down_bound = calculate_down_bound(x, nb_loops);
        cout << "\\mu_{p_best} \\in [" << down_bound << ", " << up_bound << "]" << endl;
        cout << "Keep best..." << endl;

        int min_possible = -1;
        int max_possible = -1;
        for(int i=0; i<nb_diff_proba; i++) {
            double measured = (double) nb_best[i]/nb_loops;
            double max_bound = calculate_up_bound(measured, nb_loops);

            if(min_possible < 0 && max_bound >= down_bound)
                min_possible = i;
            if(max_bound >= down_bound)
                max_possible = i;
        }
        double pmin = p_init+min_possible*p_step;
        double pmax = p_init+max_possible*p_step;

        // Store the results in a file
        probafile << n << ";";
        for(int i=0; i<nb_diff_proba; i++) {
            probafile << "("<<pmin<<","<<pbest<<","<<pmax<<")" << ";";
        }
        probafile << endl;

        auto stop_n = stop_chrono();
        int elapsed_milliseconds_n = interval_chrono(start_n, stop_n);
        cout << "n=" << n << ": We are sure at 95% that it is in [" << pmin << ", " << pmax << "], with p_best=" << pbest << "." << endl;
        cout << "... executed in " << elapsed_milliseconds_n/1000. << " seconds ";
        cout << "with " << nb_solved_problems << " resolutions";

        for(int u=0; u<3; u++)
            Beep(1568, 200);

        //delete[] nb_best;
    }
}

/// Bellow, it is just some tests of algorithm...

struct Weight {
    Weight() : Weight(0) {};
    Weight(double w) {
        weight = w;
        nb_executed = 0;
        nb_best = 0;
    }

    double weight;
    int nb_executed;
    int nb_best;
};

void optimal_proba_v2(int n_initial, int n_final, std::default_random_engine& randomizer, double p_init=0.05, double p_final=0.5, double p_step=0.01, int depth_initial=8) {
    int nb_diff_proba = floor((p_final-p_init)/p_step);
    double initial_value = 50;
    int nb_test = 20;

    double zero_limit = 2/3.;
    double pa = -1./(pow(zero_limit, 2)-2*zero_limit+1);
    double pb = -2*pa;
    double pc = 1+pa;

    ofstream probafile("probas.csv");
    for(int i=0; i<nb_diff_proba; i++)
        probafile << p_init+i*p_step << ";";
    probafile << endl;

    for(int n=n_initial; n<n_final; n++) {
        auto start_n = start_chrono();

        Weight* tab = new Weight[nb_diff_proba];
        for(int i=0; i<nb_diff_proba; i++)
            tab[i].weight = initial_value;
        double max_weight = initial_value;
        int max_weight_ind_proba = 0;
        double second_weight = initial_value;
        int second_weight_ind_proba = 1;

        int nb_solved_problems = 0;
        for(int depth=depth_initial; depth>0; depth/=2) {
            double start_with_weight = max_weight;
            for(int num_loop=1; second_weight + 2*initial_value > max_weight && max_weight < 20*initial_value+start_with_weight; num_loop++) {
                int nb_new_problems = 0;
                for(int k=0; k<20; k++) {
                    double min_avg = -1;
                    double min_proba;
                    for(int ind2=-2; ind2<nb_diff_proba; ind2++) {
                        if(ind2 == max_weight_ind_proba || ind2 == second_weight_ind_proba) continue;
                        int ind = (ind2 >= 0) ? ind2 : ((ind2 == -2) ? max_weight_ind_proba: second_weight_ind_proba);
                        if(ind % depth != 0) continue;

                        double p = p_init+ind*p_step;

                        double weight = tab[ind].weight;
                        double execution_proba = weight/max_weight;
                        if(ind-depth >= 0 && ind+depth < nb_diff_proba)
                            execution_proba = ((tab[ind-depth].weight+weight+tab[ind+depth].weight)/3)/max_weight;
                        else if(ind-depth >= 0)
                            execution_proba = ((tab[ind-depth].weight + weight)/2)/max_weight;
                        else if(ind+depth < nb_diff_proba)
                            execution_proba = ((weight + tab[ind+depth].weight)/2)/max_weight;

                        execution_proba = (2*execution_proba > 1) ? pa*pow(execution_proba, 2)+pb*execution_proba+pc : 0;
                        std::bernoulli_distribution distribution(execution_proba);

                        if(distribution(randomizer)) {
                            tab[ind].nb_executed++;
                            EA algo(1, 1, PLUS, p, randomizer);
                            MasterMind mm(n, n, randomizer);
                            algo.set_problem(&mm);
                            Statistic stat = avg_calls(&algo, nb_test, min_avg+1);
                            double avg = stat.average;

                            if(min_avg < 0 || min_avg > avg) {
                                min_avg = avg;
                                min_proba = p;
                            }

                            nb_new_problems += stat.nb;
                        }
                    }
                    int ind_p = round((min_proba - p_init)/p_step);
                    tab[ind_p].weight++;

                    if(tab[ind_p].weight > max_weight) {
                        max_weight = tab[ind_p].weight;
                        max_weight_ind_proba = ind_p;
                    }

                    tab[max_weight_ind_proba].nb_best++;
                }

                auto stop_loop = stop_chrono();
                int elapsed_milliseconds_loop = interval_chrono(start_n, stop_loop);

                second_weight = -1;
                for(int i=0; i<nb_diff_proba; i++)
                    if(i != max_weight_ind_proba)
                        if(second_weight < 0 || second_weight < tab[i].weight) {
                            second_weight = tab[i].weight;
                            second_weight_ind_proba = i;
                        }

                nb_solved_problems += nb_new_problems;

                cout << "\r";
                cout << "n=" << n << ", depth=" << depth << ", num_loop=" << num_loop << ": ";
                cout << p_init+max_weight_ind_proba*p_step << " (" << max_weight << "), ";
                cout << p_init+second_weight_ind_proba*p_step << " (" << second_weight << "), ";
                cout << "executed in " <<elapsed_milliseconds_loop/1000. << " seconds [.+" << nb_new_problems << "=" << nb_solved_problems << "]";
                cout << "    ";
            }

            int fdepth = depth / 2;
            if(fdepth) {
                for(int ind=0; ind<nb_diff_proba; ind++) {
                    if(ind%depth != 0 && ind%fdepth == 0) {
                        if(ind+fdepth<nb_diff_proba)
                            tab[ind].weight = max(tab[ind-fdepth].weight, tab[ind+fdepth].weight);
                        else
                            tab[ind].weight = tab[ind-fdepth].weight;
                    }
                }
                second_weight = max_weight;
            }
        }

        double p_best = p_init+max_weight_ind_proba*p_step;

        double mean = 0;
        for(int i=0; i<nb_diff_proba; i++)
            mean += tab[i].nb_executed ? ((double) tab[i].nb_best/tab[i].nb_executed) * (p_init+i*p_step) : 0;

        double stddev = 0;
        for(int i=0; i<nb_diff_proba; i++)
            stddev += tab[i].nb_executed ? ((double) tab[i].nb_best/tab[i].nb_executed) * pow(p_init+i*p_step, 2) : 0;
        stddev -= pow(mean, 2);
        cout << '%' << mean << "," << stddev  << endl;
        stddev = sqrt(stddev);


        for(int i=0; i<nb_diff_proba; i++)
            probafile << (tab[i].nb_executed ? ((double) tab[i].nb_best/tab[i].nb_executed) * pow(p_init+i*p_step, 2) : 0) << ";";
        probafile << endl;

        auto stop_n = stop_chrono();
        int elapsed_milliseconds_n = interval_chrono(start_n, stop_n);
        cout << '\r' << "n=" << n << ", p_best=" << p_best << " (" << stddev << "), ";
        cout << "executed in " << elapsed_milliseconds_n/1000. << " seconds ";
        cout << "with " << nb_solved_problems << " resolutions";
        cout << "                           " << endl;

//        for(int u=0; u<3; u++)
//            Beep(1568, 200);

        delete[] tab;
    }
}

void optimal_proba(int n_initial, int n_final, std::default_random_engine& randomizer, double p_init=0.05, double p_final=0.5, double p_step=0.01) {
    for(int n=n_initial; n<n_final; n++) {
        vector<double> avg_list;
        vector<double> probas;

        auto start = start_chrono();
        for(double p=p_init; p<p_final; p+=p_step) {
            EA algo(1, 1, PLUS, p, randomizer);
            MasterMind mm(n, n, randomizer);
            algo.set_problem(&mm);
            Statistic stat = avg_calls(&algo);
            double avg = stat.average;
            avg_list.push_back(avg);
            probas.push_back(p);
        }
        auto stop = stop_chrono();
        int elapsed_milliseconds = interval_chrono(start, stop);

        int min_ind = 0;
        for(unsigned int i=1; i<avg_list.size(); i++)
            if(avg_list[min_ind] > avg_list[i])
                min_ind = i;
        double p_best = probas[min_ind];
        cout << "n=" << n << ": " << p_best << " in " << elapsed_milliseconds/1000. << " seconds " << "(" << 100*floor((p_final-p_init)/p_step) << " problems)" << endl;
    }
}

#endif // SCRIPTS_H_INCLUDED