/*****************************************************************//**
 * \file   CliqueInequalities.h
 * \brief  Finds and adds cliques inequalities to the node on the conflicting columns containing at least two subjects in common.
 * 
 * \author talmefta
 * \date   June 2022
 *********************************************************************/
#ifndef CLIQUEINEQUALITIES_H
#define CLIQUEINEQUALITIES_H

#include "Header.h"
#include "CutGeneratorAbstract.h"
#include "RestrictedMasterProblem.h"
#include "PricingProblem.h"
#include "ColumnGeneration.h"

/**
 * \class CliqueInequalities.
 * \brief  Finds and adds cliques inequalities to the node on the conflicting columns containing at least two subjects in common.
 */
class CliqueInequalities: public CutGeneratorAbstract<RestrictedMasterProblem, PricingProblem, ColumnGeneration<RestrictedMasterProblem, PricingProblem>, GRBConstr> {

protected:

    const vector<Vertex>& vertices;
    int maxNumberPerIteration;
    bool isCliqueNotFound;

    //stats
    double avg_length;

    /**
     * Tries to greedy grow the clique in parameter by other neighboring vertices in the graph.
     * \details The function is a pure implementation of the procedure GreedyGrow in https://doi.org/10.1007/s00500-019-03769-y
     * 
     * \param conflict_graph: map<int, vector<int>>&
     * \param weights: map<int, double>&
     * \param clique: const vector<int>&
     * \return Augmented clique if any: vector<int>
     */
    vector<int> GreedyGrow(map<int, vector<int>>& conflict_graph, map<int, double>& weights, const vector<int>& clique)
    {
        /* Construct the graph induced by the neighbors of all vertices included in the clique in the same time.*/
        map<int, vector<int>> induced_graph;
        vector<vector<int>> neighbors_of_clique_nodes;

        //Construct for every node in clique, the set of its neighbors in the conflict graph
        for (auto& clique_node : clique)
        {
            neighbors_of_clique_nodes.push_back({});
            for (auto& conflicting : conflict_graph[clique_node])
                if (find(clique.begin(), clique.end(), conflicting) == clique.end())
                    neighbors_of_clique_nodes.back().push_back(conflicting);
        }
        //Constrcut the set of nodes that are neighbors to all of the clique in the same time.
        vector<int> neighbors;
        
        if (!neighbors_of_clique_nodes.empty())
        {
            neighbors = neighbors_of_clique_nodes.front();
            sort(neighbors.begin(), neighbors.end());
            vector<int>::iterator it;
            vector<int> common(neighbors.size());
            for (int k = 1; k < neighbors_of_clique_nodes.size(); k++)
            {
                sort(neighbors_of_clique_nodes[k].begin(), neighbors_of_clique_nodes[k].end());
                it = set_intersection(neighbors.begin(), neighbors.end(),
                    neighbors_of_clique_nodes[k].begin(), neighbors_of_clique_nodes[k].end(),
                    common.begin());
                common.resize(it - common.begin());
                neighbors = common;
            }
        }
        sort(neighbors.begin(), neighbors.end());
        neighbors.erase(unique(neighbors.begin(), neighbors.end()), neighbors.end());

        //Construct the set of edges involving the set of neighbors
        for (auto& neighbor : neighbors)
        {
            for (auto& conflicting_of_neighbors : conflict_graph[neighbor])
            {
                if (find(neighbors.begin(), neighbors.end(), conflicting_of_neighbors) != neighbors.end() ||
                    find(clique.begin(), clique.end(), conflicting_of_neighbors) != clique.end())
                {
                    induced_graph[neighbor].push_back(conflicting_of_neighbors);
                }
            }
        }

        /* Initialize the result clique with the clique given in parameters.*/
        vector<int> grown_clique = clique;

        /* Loop on vertices V_bar and add them to the resulting clique as soon as the V_bar is not empty.*/
        while (!neighbors.empty())
        {
            //Search the argmax of wights in S (the node that has the maximum weight among S)
            vector<int>::iterator argmax;
            double max = -INFINITY;
            for (vector<int>::iterator element = neighbors.begin(); element != neighbors.end(); element++)
            {
                if (weights[*element] > max)
                {
                    max = weights[*element];
                    argmax = element;
                }
            }
            //Grow clique by argmax
            grown_clique.push_back(*argmax);
            int j = *argmax;
            //Remove argmax and all the nodes that are not neighbors to argmax
            neighbors.erase(argmax);
            vector<int> to_erase;
            for (auto& potentiollay_not_neighbor_of_arg_max : neighbors)
            {
                if (find(induced_graph[j].begin(), induced_graph[j].end(), potentiollay_not_neighbor_of_arg_max) == induced_graph[j].end())
                    to_erase.push_back(potentiollay_not_neighbor_of_arg_max);
            }
            for (auto& erase : to_erase)
            {
                neighbors.erase(find(neighbors.begin(), neighbors.end(), erase));
            }
        }


        return grown_clique;
    }

public:
        

    /**
     * Constructor with pricing problem.
     *
     * \param instance: Instance&
     */
    CliqueInequalities(const vector<Vertex>& vertices, const int& maxNumberPerIteration): CutGeneratorAbstract(), vertices(vertices) {
        this->maxNumberPerIteration = maxNumberPerIteration;
        this->avg_length = 0.0;
        this->isCliqueNotFound = false;
    }

    /**
     * \brief Destructor.
     *
     */
    virtual ~CliqueInequalities(){}

    const double& get_avg_length() const {
        return this->avg_length;
    }
    /**
    * \brief Evaluates if the handler (the algorithm) can handle (solve) the cutting procedure based on the code of request.
    *
    * \param request : string
    * \return True if the concrete algorithm can find the cuts, false otherwise.
    */
    bool can_cut(const string& request) const override {
        return request == "CLIQUEINEQUALITIES" || request == "WHATEVER";
    }

    void reset_found_inequality() override {
        this->isCliqueNotFound = false;
    }

    struct sort_decr_by_size {
        inline bool operator() (const pair<int, unordered_set<int>>& p1, const pair<int, unordered_set<int>>& p2) {
            return p1.second.size() > p2.second.size();
        }
    };

    bool is_associated_sri_violated(Node<RestrictedMasterProblem, PricingProblem, ColumnGeneration<RestrictedMasterProblem, PricingProblem>, GRBConstr>& node,
        const vector<int>& clique, const double& violation_clique, vector<vector<int>>& violated_subset_row, double& minimal_violation_sri, double& maximal_violation_sri) {

        int wordsNumber = this->vertices.at(1).getNgSet().size();
        ofstream cliques_violated;
        string file_name = ::solutionPath + ::instanceName + "_cliques_violated.txt";
        cliques_violated.open(file_name, ios::app);

        //TO REMOVE
        cliques_violated << "CLIQUE: ";
        for (int j = 0; j < clique.size(); j++) {
            cliques_violated << clique.at(j) << " ";
            //cout << clique.at(j) << " ";
        }
        //cout << "\n";
        cliques_violated << "\t" << violation_clique << "\n";
        for (int j = 0; j < clique.size(); j++) {
            cliques_violated << node.get_restricted_master_problem()->get_columns().at(clique.at(j)).toString() << "\t" <<
                node.get_restricted_master_problem()->get_solution().at(clique.at(j)) << "\n";
        }

        //for each vertex determine the covered exchanges
        vector<pair<int, unordered_set<int>>> conflict_set;        //for each vertex the set of exchanges covered by the vertex
        unordered_map<int, int> vertices_in_set_dict;
        for (int j = 0; j < clique.size() - 1; j++) {
            const Exchange* e1 = &node.get_restricted_master_problem()->get_columns().at(clique.at(j));
            for (int k = j + 1; k < clique.size(); k++) {
                const Exchange* e2 = &node.get_restricted_master_problem()->get_columns().at(clique.at(k));
                unordered_set<int> common = e1->vertices_in_common(*e2);
                for (unordered_set<int>::const_iterator it = common.begin(); it != common.end(); it++) {
                    if (vertices_in_set_dict.find(*it) == vertices_in_set_dict.end()) {
                        vertices_in_set_dict.insert(make_pair(*it, conflict_set.size()));
                        unordered_set<int> covered_exchange = { e1->getId(),e2->getId() };
                        conflict_set.push_back(make_pair(*it, covered_exchange));
                    }
                    else {
                        int pos = vertices_in_set_dict.at(*it);
                        conflict_set.at(pos).second.insert(e1->getId());
                        conflict_set.at(pos).second.insert(e2->getId());
                    }
                }
            }
        }

        std::sort(conflict_set.begin(), conflict_set.end(), sort_decr_by_size());

        vector<int> subset;
        set<pair<int, int>> cumulative_cover;
        for (int j = 0; j < conflict_set.size(); j++) {

            bool is_merged = false;
            for (unordered_set<int>::const_iterator it1 = conflict_set.at(j).second.begin(); it1 != prev(conflict_set.at(j).second.end()); it1++) {
                for (unordered_set<int>::const_iterator it2 = it1; it2 != conflict_set.at(j).second.end(); it2++) {
                    int minimum = min(*it1, *it2);
                    int maximum = max(*it1, *it2);
                    if (*it2 != *it1 && cumulative_cover.find(make_pair(minimum, maximum)) == cumulative_cover.end()) {
                        cumulative_cover.insert(make_pair(minimum, maximum));
                        is_merged = true;
                    }
                }
            }
            if (is_merged) {
                subset.push_back(conflict_set.at(j).first);
            }

        }

        cliques_violated << "SUBSET ROW: ";
        for (int j = 0; j < subset.size(); j++) {
            cliques_violated << subset.at(j) << " ";
        }
        cliques_violated << "\n";

        //check if the associated sri is violated
        vector<int> multipliers;
        int gcd = -1;
        if (subset.size() == 3) {
            multipliers = { 1,1,1 };
            gcd = 2;
        }
        else if (subset.size() == 4) {
            multipliers = { 2,2,2,2 };
            gcd = 3;
        }
        else if (subset.size() == 5) {
            multipliers = { 1,1,1,1,1 };
            gcd = 2;
        }


        if (gcd != -1) {

            double rhs = 0.0;
            double lhs = 0.0;
            for (int j = 0; j < subset.size(); j++) {
                rhs += multipliers.at(j);
            }
            rhs /= static_cast<double>(gcd);
            rhs += ::EPSILON;
            rhs = floor(rhs);

            for (int i = 0; i < clique.size(); i++) {
                const Exchange* e = &node.get_restricted_master_problem()->get_columns().at(clique.at(i));
                double rows_contribution = 0.0;
                for (int j = 0; j < subset.size(); j++) {
                    int id = subset.at(j);
                    rows_contribution += (e->getVisit(id) * multipliers.at(j));
                    /*rows_contribution += static_cast<double>(exchange->getVisit(id)) *
                        static_cast<double>(this->multipliers.at(j));*/
                }
                rows_contribution /= gcd;
                rows_contribution = floor(rows_contribution + ::EPSILON);
                lhs += (rows_contribution * node.get_restricted_master_problem()->get_solution().at(clique.at(i)));
            }

            double violation = lhs - rhs;
            cliques_violated << "VIOLATION SRI " << violation << "\n\n";
            if (violation > EPSILON) {

                string subsetToString;
                vector<int> subset_bit(wordsNumber, 0);
                for (int j = 0; j < subset.size(); j++) {
                    int id = subset.at(j);
                    subsetToString += "_" + to_string(id);
                    const Vertex* v = &this->vertices.at(id);
                    subset_bit.at(v->getWord()) += v->getBitKey();
                }

                node.add_cut_constraint(new SubsetRowConstraint<GRBConstr>(rhs, "<=", subset, subset_bit, multipliers, gcd, subsetToString, "SUBSETROW"));
                node.get_cut_constraints().back()->set_birth(node.get_restricted_master_problem()->get_solving_rmp_iteration());
                node.get_cut_constraints().back()->set_become_effective(node.get_restricted_master_problem()->get_solving_rmp_iteration());
                node.get_restricted_master_problem()->add_cut_constraint(node.get_cut_constraints().back(), node.get_cut_constraints().size() - 1, true);
                cliques_violated.close();
                minimal_violation_sri = min(minimal_violation_sri, violation);
                maximal_violation_sri = max(maximal_violation_sri, violation);
                this->avg_length += subset.size();
                return true;
            }
            else {
                cliques_violated << "\n";
            }

        }
        
        //call the mip to see if coefficients can be found
        GRBEnv* env = new GRBEnv();	/**< Gurobi environment.*/
        env->start();
        GRBModel* model = new GRBModel(*env);	/**< Gurobi model.*/
        model->set(GRB_IntParam_LogToConsole, 0);
        model->set(GRB_IntParam_Method, 1);
        model->set(GRB_IntParam_Threads, 1);
        /*model->set(GRB_DoubleParam_OptimalityTol, ::EPSILON);
        model->set(GRB_DoubleParam_FeasibilityTol, ::EPSILON);*/
        model->set(GRB_DoubleParam_IntFeasTol, 1e-7);
        vector<GRBVar> m;	/**< one variable per multiplier.*/
        //vector<GRBVar> dev;	/**< one variable per multiplier deviation from 0.5.*/
        vector<GRBVar> y_int;	/**< one variable per exchange in the clique int.*/
        vector<GRBVar> y_fract;	/**< one variable per exchange in the clique fractional part.*/
        string name = "rhs_int";
        double ub_model = 0.99;
        GRBVar rhs_int = model->addVar(0.0, GRB_INFINITY, 0.0, GRB_INTEGER, name);
        name = "rhs_fract";
        GRBVar rhs_fract = model->addVar(0.0, ub_model, 0.0, GRB_CONTINUOUS, name);
        name = "mini";
        GRBVar mini = model->addVar(0.0, GRB_INFINITY, 0.0, GRB_CONTINUOUS, name);	/**< one variable per multiplier deviation from 0.5.*/
        name = "maxi";
        GRBVar maxi = model->addVar(0.0, GRB_INFINITY, 0.0, GRB_CONTINUOUS, name);	/**< one variable per multiplier deviation from 0.5.*/

        for (int j = 0; j < clique.size(); j++) {
            name = "y_int_" + to_string(clique.at(j));
            y_int.push_back(model->addVar(0.0, GRB_INFINITY, 0.0, GRB_INTEGER, name));
            name = "y_fract_" + to_string(clique.at(j));
            y_fract.push_back(model->addVar(0.0, ub_model, 0.0, GRB_CONTINUOUS, name));
        }
        for (int j = 0; j < subset.size(); j++) {
            name = "m_" + to_string(subset.at(j));
            m.push_back(model->addVar(0.0, 1.0, 0.0, GRB_CONTINUOUS, name));
            //name = "dev_" + to_string(subset.at(j));
            //dev.push_back(model->addVar(0.0, GRB_INFINITY, 0.0, GRB_CONTINUOUS, name));
        }

        GRBLinExpr objective1 = 0;	/**< Gurobi expression of objective.*/
        for (int j = 0; j < clique.size(); j++) {
            objective1 -= y_int.at(j) * node.get_restricted_master_problem()->get_solution().at(clique.at(j));
        }
        objective1 += rhs_int;
        GRBLinExpr objective2 = maxi - mini;	/**< Gurobi expression of objective.*/
        //GRBLinExpr objective2 = 0;	/**< Gurobi expression of objective.*/
        /*for (int j = 0; j < subset.size(); j++) {
            objective2 += dev.at(j);
        }*/

        //hierachical multi-objective
        model->set(GRB_IntAttr_ModelSense, GRB_MINIMIZE);
        name = "obj_viol";
        model->setObjectiveN(objective1, 0, 1000, 1.0, 0.0, 0.0, name);
        name = "obj_dev";
        model->setObjectiveN(objective2, 1, 1, 1.0, 0.0, 0.0, name);

        //model->setObjective(objective1, GRB_MAXIMIZE);

        for (int j = 0; j < subset.size(); j++) {
            model->addConstr(mini - m.at(j) <= 0, "mini_" + to_string(subset.at(j)));
            model->addConstr(maxi - m.at(j) >= 0, "maxi_" + to_string(subset.at(j)));
            /*model->addConstr(dev.at(j) - m.at(j) + 0.5 >= 0, "dev_1_" + to_string(subset.at(j)));
            model->addConstr(dev.at(j) + m.at(j) - 0.5 >= 0, "dev_2_" + to_string(subset.at(j)));*/
        }
        GRBLinExpr sum_rhs = 0;
        for (int j = 0; j < subset.size(); j++) {
            sum_rhs += m.at(j);
        }
        sum_rhs -= rhs_fract;
        sum_rhs -= rhs_int;
        model->addConstr(sum_rhs == 0, "rhs_floor");
        for (int j = 0; j < clique.size(); j++) {
            const Exchange* e = &node.get_restricted_master_problem()->get_columns().at(clique.at(j));
            GRBLinExpr sum = y_int.at(j);
            for (int i = 0; i < subset.size(); i++) {
                sum -= e->getVisit(subset.at(i)) * m.at(i);
            }
            sum += y_fract.at(j);
            model->addConstr(sum == 0, "exchange_floor_" + to_string(clique.at(j)));
        }

        //model->write("MIP_clique.lp");
        model->optimize();
        int solCount = model->get(GRB_IntAttr_SolCount);
        if (solCount > 0) {

            //cliques_violated << "VIOL MIP " << model->get(GRB_DoubleAttr_ObjVal) << "\n";
            model->set(GRB_IntParam_ObjNumber, 0);
            double violation_mip = -model->get(GRB_DoubleAttr_ObjNVal);
            cliques_violated << "VIOL MIP " << -model->get(GRB_DoubleAttr_ObjNVal);
            model->set(GRB_IntParam_ObjNumber, 1);
            cliques_violated << "\tDEVIATION " << model->get(GRB_DoubleAttr_ObjNVal) << "\n";

            if (violation_mip > EPSILON) {

                for (int j = 0; j < clique.size(); j++) {
                    cliques_violated << "lhs_var_" << to_string(clique.at(j)) << "\t" << y_int.at(j).get(GRB_DoubleAttr_X) << "\t" << y_fract.at(j).get(GRB_DoubleAttr_X) << endl;
                }

                //BUILD THE SRI DATA
                string subsetToString;
                vector<int> subset_bit(wordsNumber, 0);
                vector<int> multipliers_mip;
                int rhsInt = rhs_int.get(GRB_DoubleAttr_X);
                for (int j = 0; j < subset.size(); j++) {
                    int id = subset.at(j);
                    subsetToString += "_" + to_string(id);
                    const Vertex* v = &this->vertices.at(id);
                    subset_bit.at(v->getWord()) += v->getBitKey();

                    double mult = m.at(j).get(GRB_DoubleAttr_X) * ::OMEGA;
                    //cout << m.at(j).get(GRB_DoubleAttr_X) << " " << mult << " " << ceil(mult) << endl;
                    mult = ceil(mult);
                    multipliers_mip.push_back(mult);
                    cliques_violated << "mult_var_" << to_string(subset.at(j)) << "\t" << m.at(j).get(GRB_DoubleAttr_X) << "\t" << mult << endl;

                }

                cliques_violated << "RHS " << rhs_int.get(GRB_DoubleAttr_X) << "\t" << rhs_fract.get(GRB_DoubleAttr_X) << "\t" << rhsInt << "\n";

                gcd = static_cast<int>(::OMEGA);
                cliques_violated << "GCD " << gcd << "\n\n";

                //Sanity check for numerical instability
                double viol_check = 0.0;
                for (int i = 0; i < clique.size(); i++) {
                    double coeff = 0.0;
                    for (int j = 0; j < subset.size(); j++) {
                        coeff += (multipliers_mip.at(j) * node.get_restricted_master_problem()->get_columns().at(clique.at(i)).getVisit(subset.at(j)));
                    }
                    coeff /= gcd;
                    coeff += EPSILON;
                    coeff = floor(coeff);
                    viol_check += coeff * node.get_restricted_master_problem()->get_solution().at(clique.at(i));
                }
                viol_check -= rhsInt;

                if (viol_check <= ::EPSILON) {
                    cliques_violated << "SANITY CHECK VIOLATION: " << viol_check << endl;
                    delete model;
                    delete env;
                    cliques_violated.close();
                    return false;
                }

                node.add_cut_constraint(new SubsetRowConstraint<GRBConstr>(rhsInt, "<=", subset, subset_bit, multipliers_mip, gcd, subsetToString, "SUBSETROW"));
                node.get_cut_constraints().back()->set_birth(node.get_restricted_master_problem()->get_solving_rmp_iteration());
                node.get_cut_constraints().back()->set_become_effective(node.get_restricted_master_problem()->get_solving_rmp_iteration());
                node.get_restricted_master_problem()->add_cut_constraint(node.get_cut_constraints().back(), node.get_cut_constraints().size() - 1, true);
                cliques_violated.close();
                minimal_violation_sri = min(minimal_violation_sri, violation_mip);
                maximal_violation_sri = max(maximal_violation_sri, violation_mip);
                this->avg_length += subset.size();
                delete model;
                delete env;
                return true;
            }

        }
        else {
            delete model;
            delete env;
            cliques_violated << "MIP INFEASIBLE\n";
        }

        
        cliques_violated.close();
        return false;
    }

    /**
     * Tries to find clique inequality cuts to add to the node. Stores the new cuts if any in the node's cut_constraints.
     *
     * \param node: Node<class RestrictedMasterProblem, class PricingSubProblem, class ColumnGeneration>&
     */
    bool cut(Node<RestrictedMasterProblem, PricingProblem, ColumnGeneration<RestrictedMasterProblem, PricingProblem>, GRBConstr>& node) override {

        if (can_cut(node.get_restricted_master_problem()->get_cut_request()) && node.get_depth() <= ::MAXLEVELTREESEPARATION && node.get_depth() >= ::MINLEVELTREESEPARATION &&
            cut_counter < ::MAXCUTSNUMBER && separation_time < ::MAXCUTSTIME) {
            ////Consider the conflict_graph induced by the solution variables (having y_k > 0)
            ////Construct the induced conflict graph restricted to nodes with 0 < y_k < 1
            cut_iteration++;
            clock_t start = clock();

            //node.get_restricted_master_problem()->write_solution_value_into_text(0);
           
           
            ofstream cliques_violated;
            string file_name = ::solutionPath + ::instanceName + "_cliques_violated.txt";
            cliques_violated.open(file_name, ios::app);
            cliques_violated << "BaB node " << node.get_id() << " LEVEL " << node.get_depth() << "\n";
            cliques_violated.close();

            vector<int> indexes_conflict_graph;
            vector<int> indexes_restricted_conflict_graph;
            for (int k = 0; k < node.get_restricted_master_problem()->get_columns().size(); k++)
                if (node.get_restricted_master_problem()->get_solution()[k] > 0)
                {
                    if (node.get_restricted_master_problem()->get_solution()[k] < 1)
                        indexes_restricted_conflict_graph.push_back(k);
                    indexes_conflict_graph.push_back(k);
                }


            map<int, vector<int>> conflicting_for_each_variable;
            map<int, vector<int>> restricted_conflicting_for_each_variable;

            for (int k = 0; k < indexes_conflict_graph.size(); k++)
            {
                conflicting_for_each_variable[indexes_conflict_graph[k]] = {};
                for (int l = 0; l < indexes_conflict_graph.size(); l++)
                    if (l != k)
                        if (node.get_restricted_master_problem()->get_columns()[indexes_conflict_graph[k]].conflicts_with(node.get_restricted_master_problem()->get_columns()[indexes_conflict_graph[l]])) {
                            conflicting_for_each_variable[indexes_conflict_graph[k]].push_back(indexes_conflict_graph[l]);
                        }
            }


            for (int k = 0; k < indexes_restricted_conflict_graph.size(); k++)
            {
                restricted_conflicting_for_each_variable[indexes_restricted_conflict_graph[k]] = {};
                for (int l = 0; l < indexes_restricted_conflict_graph.size(); l++)
                    if (l != k)
                        if (node.get_restricted_master_problem()->get_columns()[indexes_restricted_conflict_graph[k]].conflicts_with(node.get_restricted_master_problem()->get_columns()[indexes_restricted_conflict_graph[l]]))
                            restricted_conflicting_for_each_variable[indexes_restricted_conflict_graph[k]].push_back(indexes_restricted_conflict_graph[l]);
            }

            //Construct for every node its weight that is the size of its neighbors in the conflict graph times solution value in RMP
            map<int, double> weights;
            for (auto variable : conflicting_for_each_variable)
            {
                if (node.get_restricted_master_problem()->get_solution()[variable.first] < 1)
                    weights.insert({ variable.first, node.get_restricted_master_problem()->get_solution()[variable.first] * variable.second.size() });
                else
                    weights.insert({ variable.first, -::OMEGA });

            }

            ////Sort nodes by nondecreasing value of f(i)
            vector<pair<int, double>> pairs;
            for (auto itr = weights.begin(); itr != weights.end(); ++itr)
                pairs.push_back(*itr);

            sort(pairs.begin(), pairs.end(), [=](std::pair<int, double>& a, std::pair<int, double>& b)
                {
                    return a.second > b.second;
                }
            );

            //Variables sorted by increasing value of weights
            vector<int> ordered_nodes;
            for (auto& pair : pairs)
                ordered_nodes.push_back(pair.first);


            vector<int> S;
            vector<int> neighbors;
            double UB_viol;
            double sum;
            vector<vector<int>> violated_subset_row;
            int number_per_it = 0;
            vector<vector<int>> violated_cliques;
            vector<double> violations;

            double minimal_violation = 1e5;
            double maximal_violation = 0;
            double minimal_violation_sri = 1e5;
            double maximal_violation_sri = 0;
            //The main loop
            for (int& node_i : ordered_nodes) {
                //construct the set S of neighbors j of node_i in restricted conflict graph with f(j) > 0
                S.clear();
                for (auto& neighbor : restricted_conflicting_for_each_variable[node_i])
                    if (weights[neighbor] > 0)
                        S.push_back(neighbor);

                //Compute upper bound of the violation based on the values of f(j), j in the neighbors of node_i in conflict graph
                sum = 0;
                for (auto& neighbor : conflicting_for_each_variable[node_i])
                    sum += node.get_restricted_master_problem()->get_solution()[neighbor];

                UB_viol = node.get_restricted_master_problem()->get_solution()[node_i] + sum - 1.0;

                //Search for violated cliques containing node_i
                while (UB_viol > VIOLATIONEPSILON && !S.empty()) {
                    //Search the argmax of weights in S (the node that has the maximum weight among S)
                    int argmax;
                    double max = -INFINITY;
                    for (auto& element : S)
                    {
                        if (weights[element] > max)
                        {
                            max = weights[element];
                            argmax = element;
                        }
                    }

                    //Try to grow the clique {node_i, argmax}
                    vector<int> clique = GreedyGrow(conflicting_for_each_variable, weights, { node_i, argmax });

                    //S = S\Clique

                    for (vector<int>::iterator node = clique.begin(); node != clique.end(); node++)
                    {
                        vector<int>::iterator found = find(S.begin(), S.end(), *node);
                        if (found != S.end())
                            S.erase(found);
                    }
                    //E' = E'\E'(Clique)
                    for (auto& element : clique)
                    {
                        for (auto& node : restricted_conflicting_for_each_variable)
                        {
                            vector<int>::iterator element_found = find(node.second.begin(), node.second.end(), element);
                            if (element_found != node.second.end())
                                node.second.erase(element_found);
                        }
                    }

                    double sum_weights_clique_except_node_i = 0;
                    double sum_weights_clique = 0;
                    for (auto& element : clique)
                    {

                        if (element != node_i)
                            sum_weights_clique_except_node_i += node.get_restricted_master_problem()->get_solution()[element];
                        sum_weights_clique += node.get_restricted_master_problem()->get_solution()[element];
                    }

                    UB_viol -= sum_weights_clique_except_node_i;
                    if (sum_weights_clique > 1 + VIOLATIONEPSILON) {
                        violated_cliques.push_back(clique);
                        violations.push_back(sum_weights_clique - 1);

                        if (is_associated_sri_violated(node, clique, violations.back(), violated_subset_row, minimal_violation_sri, maximal_violation_sri)) {
                            number_per_it++;
                        }

                        if (sum_weights_clique - 1 < minimal_violation)
                            minimal_violation = sum_weights_clique - 1;

                        if (sum_weights_clique - 1 > maximal_violation)
                            maximal_violation = sum_weights_clique - 1;
                        if (number_per_it >= this->maxNumberPerIteration) {
                            break;
                        }
                    }
                }

                if (number_per_it >= this->maxNumberPerIteration) {
                    break;
                }

            }

            //Add the cut constraint (clique inequality) to the node as well as to the RMP
            ofstream cliques_violated1;
            cliques_violated1.open(file_name, ios::app);
            cliques_violated1 << "CLIQUES: MIN VIOL " << minimal_violation << " MAX VIOL " << maximal_violation << "\n";
            cliques_violated1 << "SRI: MIN VIOL " << minimal_violation_sri << " MAX VIOL " << minimal_violation_sri << "\n";
           
            if (number_per_it == 0) {
                this->isCliqueNotFound = true;
            }
            cut_counter += number_per_it;

            clock_t end = clock();
            double duration = (float)(end - start) / CLOCKS_PER_SEC;
            separation_time += duration;
            cliques_violated1 << "SEPARATION TIME " << duration << "\n\n\n";
            cliques_violated1.close();
            node.set_cut_generation_time_spent(node.get_cut_generation_time_spent() + duration);
            node.set_nb_cuts_added(static_cast<int>(node.get_nb_cuts_added() + number_per_it));
            
            //stats collection
            information.push_back(new CutFamilyInformation("CliqueInequalities", node.get_id(), cut_iteration, number_per_it, node.get_dual_bound(), duration, minimal_violation_sri, maximal_violation_sri));
            node.get_restricted_master_problem()->attach(information.back());

            if (number_per_it > 0) {
                return (number_per_it > 0); //Return true if the algorithm found violated cliques
            }
            else
                return CutGeneratorAbstract<RestrictedMasterProblem, PricingProblem, ColumnGeneration<RestrictedMasterProblem, PricingProblem>, GRBConstr>::cut(node);
        }
        else
            return CutGeneratorAbstract<RestrictedMasterProblem, PricingProblem, ColumnGeneration<RestrictedMasterProblem, PricingProblem>, GRBConstr>::cut(node);
    }

};

#endif //CLIQUEINEQUALITIES_H