/*****************************************************************//**
 * \file   SubsetRowInequalities.h
 * \brief  Finds and adds subset row inequalities
 * 
 * \author Matteo Petris
 * \date   March 2023
 *********************************************************************/
#ifndef SUBSETROWINEQUALITIES_H
#define SUBSETROWINEQUALITIES_H

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

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

protected:

    const vector<Vertex>& vertices;
    vector<vector<int>> subsets;
    vector<vector<int>> subsets_bit;
    int cardinality;
    vector<int> multipliers;
    int gcd;
    int maxNumberPerIteration;
    bool isSRNotFound;

public:


    /**
     * Constructor with pricing problem.
     *
     * \param instance: Instance&
     */
    SubsetRowInequalities(const vector<Vertex>& vertices, const int& cardinality, const vector<int>& multipliers, const int& gcd, const int& maxNumberPerIteration) : 
        CutGeneratorAbstract(), vertices(vertices) {
        this->cardinality = cardinality;
        this->multipliers = multipliers;
        this->maxNumberPerIteration = maxNumberPerIteration;
        this->gcd = gcd;
        this->isSRNotFound = false;
    }

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

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

    /**
    * \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 == "SUBSETROWINEQUALITIES" || request == "WHATEVER";
    }

    /**
    * \brief build the subsets of agents 
    *
    * \param request : string
    * \return True if the concrete algorithm can find the cuts, false otherwise.
    */
    //void buildSubsets(const Instance& instance) {

    //    for (int i = 1; i < instance.getVertexNumber(); i++) {
    //        const Vertex* v_i = &instance.getVertices().at(i);
    //        for (int j = i + 1; j < instance.getVertexNumber(); j++) {
    //            const Vertex* v_j = &instance.getVertices().at(j);
    //            int are_i_j_in_exchange = instance.can_j_be_in_an_exchange_after_i(i, j) || instance.can_j_be_in_an_exchange_after_i(j, i);
    //            for (int k = j + 1; k < instance.getVertexNumber(); k++) {
    //                const Vertex* v_k = &instance.getVertices().at(k);
    //                int are_j_k_in_exchange = instance.can_j_be_in_an_exchange_after_i(j, k) || instance.can_j_be_in_an_exchange_after_i(k, j);
    //                int are_i_k_in_exchange = instance.can_j_be_in_an_exchange_after_i(i, k) || instance.can_j_be_in_an_exchange_after_i(k, i);
    //                if (are_i_j_in_exchange + are_j_k_in_exchange + are_i_k_in_exchange >= 2) {
    //                    //cout << "ij " << are_i_j_in_exchange << " jk " << are_j_k_in_exchange << " ik " << are_i_k_in_exchange << "\n";
    //                    vector<int> subset = { i,j,k };
    //                    vector<int> subset_bit(instance.getWordsNumber(), 0);
    //                    subset_bit.at(v_i->getWord()) += v_i->getBitKey();
    //                    subset_bit.at(v_j->getWord()) += v_j->getBitKey();
    //                    subset_bit.at(v_k->getWord()) += v_k->getBitKey();
    //                    this->subsets.push_back(subset);
    //                    this->subsets_bit.push_back(subset_bit);
    //                    //::appendToFile("subsets.log", to_string(i) + "\t" + to_string(j) + "\t" + to_string(k) + "\n");
    //                }
    //            }
    //        }
    //    }

    //    /*for (int i = 1; i < instance.getVertexNumber(); i++) {
    //        const Vertex* v_i = &instance.getVertices().at(i);
    //        for (int j = i + 1; j < instance.getVertexNumber(); j++) {
    //            const Vertex* v_j = &instance.getVertices().at(j);
    //            for (int k = j + 1; k < instance.getVertexNumber(); k++) {
    //                const Vertex* v_k = &instance.getVertices().at(k);
    //                vector<int> subset = { i,j,k };
    //                vector<int> subset_bit(instance.getWordsNumber(), 0);
    //                subset_bit.at(v_i->getWord()) += v_i->getBitKey();
    //                subset_bit.at(v_j->getWord()) += v_j->getBitKey();
    //                subset_bit.at(v_k->getWord()) += v_k->getBitKey();
    //                this->subsets.push_back(subset);
    //                this->subsets_bit.push_back(subset_bit);
    //            }
    //        }
    //    }*/

    //}

    /**
     * Tries to find subset row 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
            clock_t start = clock();

            //node.get_restricted_master_problem()->write_solution_value_into_text(0);

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

            /****************************************DETECT FRACTIONAL COLUMNS AND VERTICES TO CONSIDER IN THE SSI*************************************************************************/
            vector<int> positiveColumns;
            vector<int> fractional_vertices;
            unordered_set<int> fractional_vertices_set;
            for (int k = 0; k < node.get_restricted_master_problem()->get_columns().size(); k++) {
                double val = node.get_restricted_master_problem()->get_solution().at(k);
                if (val > ::EPSILON &&
                    ::strictly_less(val, 1.0)) {
                    positiveColumns.push_back(k);
                    const Exchange* exchange = &node.get_restricted_master_problem()->get_columns().at(k);
                    if (val > 0.3) {
                        for (int i = 1; i < exchange->getSequence().size(); i++) {
                            int id = exchange->getSequence().at(i);
                            if (fractional_vertices_set.find(id) == fractional_vertices_set.end()) {
                                fractional_vertices.push_back(id);
                                fractional_vertices_set.insert(id);
                            }
                        }
                    }
                }
            }

            /*******************************************************************BUILD SSI SUBSETS*************************************************************************/
            vector<vector<int>> subsets;
            vector<int> subsets_pos;
            for (int i = 0; i < fractional_vertices.size(); i++) {
                const Vertex* v_i = &this->vertices.at(fractional_vertices.at(i));
                for (int j = i + 1; j < fractional_vertices.size(); j++) {
                    const Vertex* v_j = &this->vertices.at(fractional_vertices.at(j));
                    for (int k = j + 1; k < fractional_vertices.size(); k++) {
                        const Vertex* v_k = &this->vertices.at(fractional_vertices.at(k));
                        vector<int> subset = { v_i->getId(),v_j->getId(),v_k->getId() };
                        /*vector<int> subset_bit(wordsNumber, 0);
                        subset_bit.at(v_i->getWord()) += v_i->getBitKey();
                        subset_bit.at(v_j->getWord()) += v_j->getBitKey();
                        subset_bit.at(v_k->getWord()) += v_k->getBitKey();*/
                        subsets_pos.push_back(subsets.size());
                        subsets.push_back(subset);
                        //subsets_bit.push_back(subset_bit);
                    }
                }
            }

            //cout << "subset size " << subsets.size() << " fract vert " << fractional_vertices.size() << endl;

            /*******************************************************************SEPARATE SSI BY ENUMERATION*************************************************************************/
            unsigned seed = chrono::system_clock::now().time_since_epoch().count();
            shuffle(subsets_pos.begin(), subsets_pos.end(), default_random_engine(::generator));
            //cout << "SUBSET NUMBER " << subsets.size() << endl;
            vector<pair<double, int>> violated_subset_row;
            int violation_above = 0;
            clock_t start1 = clock();
            clock_t end1 = clock();
            double duration1 = (float)(end1 - start1) / CLOCKS_PER_SEC;
            int counter = 0;
            int maxSubsetNumber = min(subsets_pos.size(), ::SSIMAXSUBSETNUMBER);
            while ((duration1 < ::MAXTIMESEPARATIONROUND) && (violation_above < this->maxNumberPerIteration) && counter < maxSubsetNumber) {

                int sId = subsets_pos.at(counter);
                counter++;

                double lhs = 0.0;
                for (int k = 0; k < positiveColumns.size(); k++) {
                    int col = positiveColumns.at(k);
                    const Exchange* exchange = &node.get_restricted_master_problem()->get_columns().at(col);
                    //double rows_contribution = compute_subsetRow_contribution(exchange, this->subsets.at(i));
                    //cout << exchange->toString() << endl;
                    double rows_contribution = 0.0;
                    for (int j = 0; j < subsets.at(sId).size(); j++) {
                        int id = subsets.at(sId).at(j);
                        rows_contribution += (exchange->getVisit(id) * this->multipliers.at(j));
                        /*rows_contribution += static_cast<double>(exchange->getVisit(id)) *
                            static_cast<double>(this->multipliers.at(j));*/
                    }
                    rows_contribution /= this->gcd;
                    rows_contribution = floor(rows_contribution + ::EPSILON);

                    lhs += (rows_contribution * node.get_restricted_master_problem()->get_solution().at(col));

                    /*cout << col << "\t" << rows_contribution << "\t" << node.get_restricted_master_problem()->get_solution().at(col) << "\n";
                    if (node.get_restricted_master_problem()->get_solution().at(col) > EPSILON && rows_contribution > EPSILON) {
                        cout << exchange->toString() << endl;
                    }*/
                }
                /*if (lhs > ::EPSILON) {
                    for (int j = 0; j < this->subsets.at(i).size(); j++) {
                        cout << this->subsets.at(i).at(j) << " ";
                    }
                    cout << "\t" << lhs << "\n";
                    for (int j = 0; j < cols.size(); j++) {
                        cout << node.get_restricted_master_problem()->get_columns().at(cols.at(j)).toString() << "\n";
                        cout << "CONTR " << contrib.at(j)  << "\n";
                    }
                    cout << "\n";
                }*/
                double violation = lhs - rhs;
                if (::greater(violation, ::VIOLATIONEPSILON)) {
                    violated_subset_row.push_back(make_pair(violation, sId));
                    violation_above++;
                }


                end1 = clock();
                duration1 = (float)(end1 - start1) / CLOCKS_PER_SEC;

            }

            ///*******************************************************************BUILD SSI SUBSETS*************************************************************************/
            //vector<pair<double, vector<int>>> subsets;
            //vector<int> subsets_pos;
            //for (int i = 0; i < fractional_vertices.size(); i++) {
            //    const Vertex* v_i = &this->vertices.at(fractional_vertices.at(i));
            //    for (int j = i + 1; j < fractional_vertices.size(); j++) {
            //        const Vertex* v_j = &this->vertices.at(fractional_vertices.at(j));

            //        double lhs = 0.0;
            //        for (int k = 0; k < positiveColumns.size(); k++) {
            //            int col = positiveColumns.at(k);
            //            const Exchange* exchange = &node.get_restricted_master_problem()->get_columns().at(col);
            //            lhs += (exchange->getVisit(v_i->getId()) * exchange->getVisit(v_j->getId()) * node.get_restricted_master_problem()->get_solution().at(col));
            //        }

            //        subsets_pos.push_back(subsets.size());
            //        vector<int> subset = { v_i->getId(),v_j->getId() };
            //        subsets.push_back(make_pair(lhs, subset));

            //    }
            //}

            ///*******************************************************************SEPARATE SSI BY ENUMERATION*************************************************************************/
            //sort(subsets.rbegin(), subsets.rend());
           
            //vector<pair<double, vector<int>>> violated_subset_row;
            //int violation_above = 0;
            //clock_t start1 = clock();
            //clock_t end1 = clock();
            //double duration1 = (float)(end1 - start1) / CLOCKS_PER_SEC;
            //int counter = 0;

            //vector<int> covered(this->vertices.size(), 0);

            //for (int j = 0; j < subsets.size(); j++) {

            //    int id1 = subsets.at(j).second.at(0);
            //    int id2 = subsets.at(j).second.at(1);
            //    



            //    if (covered.at(id1) == 0 && covered.at(id2) == 0 && (duration1 < ::MAXTIMESEPARATIONROUND)) {

            //        for (int i = 0; i < fractional_vertices.size(); i++) {

            //            vector<int> sub = subsets.at(j).second;
            //            int id3 = fractional_vertices.at(i);
            //            sub.push_back(id3);
            //            double lhs = 0.0;


            //            for (int k = 0; k < positiveColumns.size(); k++) {
            //                int col = positiveColumns.at(k);
            //                const Exchange* exchange = &node.get_restricted_master_problem()->get_columns().at(col);
            //                //double rows_contribution = compute_subsetRow_contribution(exchange, this->subsets.at(i));
            //                //cout << exchange->toString() << endl;
            //                double rows_contribution = 0.0;
            //                for (int l = 0; l < sub.size(); l++) {
            //                    int id = sub.at(l);
            //                   // cout << id << " ";
            //                    rows_contribution += (exchange->getVisit(id) * this->multipliers.at(l));
            //                    /*rows_contribution += static_cast<double>(exchange->getVisit(id)) *
            //                        static_cast<double>(this->multipliers.at(j));*/
            //                }
            //                rows_contribution /= this->gcd;
            //                rows_contribution = floor(rows_contribution + ::EPSILON);
            //               // cout << " f\n";

            //                lhs += (rows_contribution * node.get_restricted_master_problem()->get_solution().at(col));

            //                /*cout << col << "\t" << rows_contribution << "\t" << node.get_restricted_master_problem()->get_solution().at(col) << "\n";
            //                if (node.get_restricted_master_problem()->get_solution().at(col) > EPSILON && rows_contribution > EPSILON) {
            //                    cout << exchange->toString() << endl;
            //                }*/
            //            }

            //            double violation = lhs - rhs;
            //            if (::greater(violation, ::VIOLATIONEPSILON)) {
            //                violated_subset_row.push_back(make_pair(violation, sub));
            //                violation_above++;
            //            }

            //        }

            //    }

            //    covered.at(id1)++;
            //    covered.at(id2)++;

            //    end1 = clock();
            //    duration1 = (float)(end1 - start1) / CLOCKS_PER_SEC;

            //}

            //cout << "NEW SEPARATION ROUND\n";

            //vector<pair<double, int>> violated_subset_row;
            //int violation_above = 0;
            //for (int i = 0; i < subsets.size(); i++) {
            //    double lhs = 0.0;
            //    /*vector<int> cols;
            //    vector<double> contrib;*/
            //    /*bool check = false;
            //    if (this->subsets.at(i).at(0) == 78 && this->subsets.at(i).at(1) == 100 && this->subsets.at(i).at(2) == 122) {
            //        check = true;
            //    }*/
            //   //cout << "SUBSET\t" << subsets.at(i).at(0) << "\t" << subsets.at(i).at(1) << "\t" << subsets.at(i).at(2) << endl;
            //    for (int k = 0; k < positiveColumns.size(); k++) {
            //        int col = positiveColumns.at(k);
            //        const Exchange* exchange = &node.get_restricted_master_problem()->get_columns().at(col);
            //        //double rows_contribution = compute_subsetRow_contribution(exchange, this->subsets.at(i));
            //        //cout << exchange->toString() << endl;
            //        double rows_contribution = 0.0;
            //        for (int j = 0; j < subsets.at(i).size(); j++) {
            //            int id = subsets.at(i).at(j);
            //            rows_contribution += (exchange->getVisit(id) * this->multipliers.at(j));
            //            /*rows_contribution += static_cast<double>(exchange->getVisit(id)) *
            //                static_cast<double>(this->multipliers.at(j));*/
            //        }
            //        rows_contribution /= this->gcd;
            //        rows_contribution = floor(rows_contribution + ::EPSILON);

            //        lhs += (rows_contribution * node.get_restricted_master_problem()->get_solution().at(col));
            //        
            //        /*cout << col << "\t" << rows_contribution << "\t" << node.get_restricted_master_problem()->get_solution().at(col) << "\n";
            //        if (node.get_restricted_master_problem()->get_solution().at(col) > EPSILON && rows_contribution > EPSILON) {
            //            cout << exchange->toString() << endl;
            //        }*/
            //    }
            //    //if (lhs - rhs > ::EPSILON) {
            //    //    for (int j = 0; j < subsets.at(i).size(); j++) {
            //    //        cout << subsets.at(i).at(j) << " ";
            //    //    }
            //    //    cout << "\t" << lhs << "\n";
            //    //    /*for (int j = 0; j < cols.size(); j++) {
            //    //        cout << node.get_restricted_master_problem()->get_columns().at(cols.at(j)).toString() << "\n";
            //    //        cout << "CONTR " << contrib.at(j)  << "\n";
            //    //    }*/
            //    //    cout << "\n";
            //    //}
            //    double violation = lhs - rhs;
            //    if (::greater(violation, ::VIOLATIONEPSILON)) {
            //        violated_subset_row.push_back(make_pair(violation, i));
            //        violation_above++;
            //    }
            //    if (violation_above >= this->maxNumberPerIteration) {
            //        break;
            //    }
            //}

            sort(violated_subset_row.rbegin(), violated_subset_row.rend());
            //Add the cut constraint to the node as well as to the RMP
            ofstream subsetrow_violated;
            string file_name = ::solutionPath + ::instanceName + "_subsetrow_violated.txt";
            subsetrow_violated.open(file_name, ios::app);
            int insertedSSI = 0;
            double minViolation = ::OMEGA;
            double maxViolation = -::OMEGA;
            int wordsNumber = this->vertices.at(1).getNgSet().size();
            for (int i = 0; i < violated_subset_row.size(); i++) {
                int pos = violated_subset_row.at(i).second;
                string subsetToString;
                for (int j = 0; j < subsets.at(pos).size(); j++) {
                    subsetToString += "_" + to_string(subsets.at(pos).at(j));
                    subsetrow_violated << subsets.at(pos).at(j) << " ";
                }
                /*for (int j = 0; j < violated_subset_row.at(i).second.size(); j++) {
                    subsetToString += "_" + to_string(violated_subset_row.at(i).second.at(j));
                    subsetrow_violated << violated_subset_row.at(i).second.at(j) << " ";
                }*/
                double viol = violated_subset_row.at(i).first;
                subsetrow_violated << "\t" << viol << "\n";
                minViolation = min(minViolation, viol);
                maxViolation = max(maxViolation, viol);

                //build the subset encoded in bits
                vector<int> subset_bit(wordsNumber, 0);
                for (int j = 0; j < subsets.at(pos).size(); j++) {
                    const Vertex* v = &this->vertices.at(subsets.at(pos).at(j));
                    subset_bit.at(v->getWord()) += v->getBitKey();
                }
                /*for (int j = 0; j < violated_subset_row.at(i).second.size(); j++) {
                    const Vertex* v = &this->vertices.at(violated_subset_row.at(i).second.at(j));
                    subset_bit.at(v->getWord()) += v->getBitKey();
                }*/

                node.add_cut_constraint(new SubsetRowConstraint<GRBConstr>(rhs, "<=", subsets.at(pos), subset_bit, this->multipliers, this->gcd, subsetToString, "SUBSETROW"));
                //node.add_cut_constraint(new SubsetRowConstraint<GRBConstr>(rhs, "<=", violated_subset_row.at(i).second, subset_bit, this->multipliers, this->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);
                insertedSSI++;
                if (insertedSSI >= this->maxNumberPerIteration) {
                    break;
                }
            }
            cut_counter += insertedSSI;

            if (insertedSSI == 0) {
                this->isSRNotFound = true;
            }

            /*if (::isToPrint) {
                cout << "\n\nCUTS ADDED\t\t" << setprecision(9) << minViolation << "\t" << setprecision(9) << maxViolation << endl;
            }*/
            subsetrow_violated << "MIN VIOL " << minViolation << " MAX VIOL " << maxViolation << "\n\n";
            subsetrow_violated << "\n\n\n";
            
            clock_t end = clock();
            double duration = (float)(end - start) / CLOCKS_PER_SEC;
            separation_time += duration;
            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() + insertedSSI));
            subsetrow_violated << "TIME SEPARATION ROUND\t" << duration << endl;
            subsetrow_violated.close();

            //stats collection
            information.push_back(new CutFamilyInformation("SubsetRowInequalities", node.get_id(), cut_iteration, insertedSSI, node.get_dual_bound(), duration, minViolation, maxViolation));
            node.get_restricted_master_problem()->attach(information.back());

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

    }

    double compute_subsetRow_contribution(const Exchange* e, const vector<int>& subset) const {
        double rows_contribution = 0.0;
        for (int j = 0; j < subset.size(); j++) {
            int id = subset.at(j);
            int prod = e->getVisit(id)* this->multipliers.at(j);
            rows_contribution += prod;
            /*rows_contribution += static_cast<double>(exchange->getVisit(id)) *
                static_cast<double>(this->multipliers.at(j));*/
        }
        rows_contribution /= static_cast<double>(this->gcd);
        rows_contribution = floor(rows_contribution + ::EPSILON);
        return rows_contribution;
    }

};

#endif //SUBSETROWINEQUALITIES_H