#include "RestrictedMasterProblem.h"

RestrictedMasterProblem::RestrictedMasterProblem(Instance& instance) :
	RestrictedMasterProblemInterface<Instance, vector<double>, DualSolution, vector<Exchange>, GRBConstr>(instance), env(nullptr), model(nullptr) {

	this->objective_value = GRB_INFINITY;
	this->objective_value_int = -::OMEGA;
	this->env = new GRBEnv();

	//this->env->set("LogFile", "RestrictedMasterProblem.log");

	this->env->start();
	this->model = new GRBModel(*env);

	this->model->set(GRB_DoubleParam_MIPGap, 0);
	this->model->set(GRB_IntParam_LogToConsole, 0);
	this->model->set(GRB_IntParam_Method, 1);
	this->model->set(GRB_IntParam_Threads, 1);
}

RestrictedMasterProblem::~RestrictedMasterProblem() {
}

void RestrictedMasterProblem::construct() {
	try {
		/****** Variables ******/
		addLambdaVariables();

		/***** Constraints *****/
		addPackingConstraints();
		addUpperBoundExchangeNumber();

		/****** Objective function ******/
		objective_function();

		/****** Export model ******/
		//export_model();
	}
	catch (GRBException e) {
		cout << "[construct_model()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[construct_model()] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::solve() {
	try {
		clock_t		start_solving_rmp = clock();
		//this->export_model();
		this->model->optimize();
		clock_t		end_solving_rmp = clock();
		double solving_rmp_duration = (float)(end_solving_rmp - start_solving_rmp) / CLOCKS_PER_SEC;
		this->solving_time = solving_rmp_duration;
		//cout << "TL\t" << timeLimit << endl;
	}
	catch (GRBException e) {
		cout << "[RestrictedMasterProblem::solve()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[RestrictedMasterProblem::solve()] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::extract_solution_value() {
	bool status;
	try {

		this->solution = vector<double>(this->columns.size(), 0);
		status = is_infeasible();
		for (int k = 0; k < this->columns.size(); k++) {
			this->solution.at(k) = this->lambda.at(k).get(GRB_DoubleAttr_X);
		}
		this->objective_value = this->model->get(GRB_DoubleAttr_ObjVal);

		/*vector<double> summ(this->instance.getVertices().size(), 0.0);
		for (int i = 1; i < this->instance.getVertices().size(); i++) {
			for (int k = 0; k < this->columns.size(); k++) {
				summ.at(i) += this->columns.at(k).getVisit(i) * this->solution.at(k);
			}
			if (::strictly_greater(summ.at(i), 1.0)) {
				cout << "INFEASIBLE SOL\n";
			}
		}

		double ss = 0.0;
		for (int k = 0; k < this->solution.size(); k++) {
			ss += this->solution.at(k);
		}
		if (::strictly_greater(ss, this->instance.getVertices().size() / 2.0)) {
			cout << "INFEASIBLE SOL\n";
		}

		if (::equal(vaal, this->objective_value)) {
			cout << "REC " << vaal << " STORED " << this->objective_value << endl;
			cout << "\n";
			//exit(EXIT_FAILURE);
		}*/
	}
	catch (GRBException e) {
		cout << "[RestrictedMasterProblem::extract_solution_value()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
		cout << (is_infeasible() ? " Infeasible " : " Something else ") << endl;
	}
	catch (...) {
		cout << "[RestrictedMasterProblem::extract_solution_value()] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::extract_dual_solution_value() {
	try {

		//cout << "duals\n";
		/*if (isToPrint) {
			::appendToFile("fileCheck.log", "DUALS:");
		}*/
		this->dualSolutionPi = vector<double>(this->instance.getVertexNumber(), 0.0);
		this->dualSolutionGamma = vector<double>(this->branching_constraints.size(), 0.0);
		this->dualSolutionSigma = vector<double>(this->cut_constraints.size(), 0.0);

		//cout << "PACKING\n";
		for (int i = 0; i < this->instance.getVertices().size(); i++) {
			const Vertex* v = &this->instance.getVertices().at(i);
			if (v->getType() != 'D') {
				this->dualSolutionPi.at(v->getId()) = this->packingConstraints.at(v->getId()).get(GRB_DoubleAttr_Pi);
				//if (this->dualSolutionPi.at(v->getId()) > 0 && ::isToPrint) {
				//	string check = this->packingConstraints.at(v->getId()).get(GRB_StringAttr_ConstrName) + "\t" + to_string(v->getId()) + "\t" + to_string(this->dualSolutionPi.at(v->getId())) + "\n";
				//	::appendToFile("fileCheck.log", check);
				//	//cout << this->packingConstraints.at(v->getId()).get(GRB_StringAttr_ConstrName) << "\t" << v->getId() << "\t" << this->dualSolutionPi.at(v->getId()) << endl;
				//}
			}
		}
		/*for (int i = 0; i < this->packingConstraints.size(); i++) {
			this->dualSolutionPi.at(i + 1) = this->packingConstraints.at(i).get(GRB_DoubleAttr_Pi);
			if (this->dualSolutionPi.at(i + 1) > 0) {
				cout << this->packingConstraints.at(i).get(GRB_StringAttr_ConstrName) << "\t" << i+1 << "\t" << this->dualSolutionPi.at(i+1) << endl;
			}
		}*/

		//cout << "BRANCH\n";
		for (int i = 0; i < this->branching_constraints.size(); i++) {
			this->dualSolutionGamma.at(i) = this->branching_constraints.at(i)->get_real_constraint()->get(GRB_DoubleAttr_Pi);
			/*if (this->dualSolutionGamma.at(i) != 0.0 && ::isToPrint) {
				int id1 = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_items().first;
				int id2 = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_items().second;
				string check = "branch " + to_string(id1) + "\t" + to_string(id2) + "\t" + to_string(this->dualSolutionGamma.at(i)) + "\n";
				::appendToFile("fileCheck.log", check);
			}*/
			//cout << i << "\t" << this->dualSolutionGamma.at(i) << endl;
		}

		for (int i = 0; i < this->cut_constraints.size(); i++) {
			this->dualSolutionSigma.at(i) = this->cut_constraints.at(i)->get_real_constraint()->get(GRB_DoubleAttr_Pi);
			const SubsetRowConstraint<GRBConstr>* ssi = dynamic_cast<SubsetRowConstraint<GRBConstr>*>(this->cut_constraints.at(i));
			//if (this->dualSolutionSigma.at(i) > EPSILON && ::isToPrint) {
			//	::appendToFile("fileCheck.log", ssi->toString() + "\t" + to_string(this->dualSolutionSigma.at(i)) + "\n");
			//	//cout << ssi->toString() << "\t" << this->dualSolutionSigma.at(i) << endl;
			//}
			//cout << i << "\t" << this->dualSolutionGamma.at(i) << endl;
		}
		
		this->dualSolutionBeta = this->upperBoundExchangeNumber.get(GRB_DoubleAttr_Pi);
		/*if (::isToPrint) {
			::appendToFile("fileCheck.log", "BOUND\t" + to_string(this->dualSolutionBeta) + "\n");
		}*/
		//cout << "BOUND " << this->dualSolutionBeta << endl;
		this->dual_variables = DualSolution(this->dualSolutionPi, this->dualSolutionGamma, this->dualSolutionSigma, this->dualSolutionBeta);
		//cout << "DUAL SOLUTION:\n" << this->dual_variables.toString() << endl;

	}
	catch (GRBException e) {
		cout << "[RestrictedMasterProblem::extract_dual_soltion_value()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[RestrictedMasterProblem::extract_dual_soltion_value()] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::write_solution_value() {
	try {
		cout << "Objective   value: " << setprecision(10) << objective_value << endl;
		//double c_k_part = 0;
		//double sum_alpha_i_beta_ij = 0;

		//for (int i = 0; i < instance.size(); i++)
		//{
		//	sum_alpha_i_beta_ij += dual_variables.get_alpha_i(i);
		//	//for (int j = i + 1; j < instance.size(); j++)
		//	//	sum_alpha_i_beta_ij += dual_variables.get_beta_ij(i, j);
		//}
		//cout << " 2 * sum_alpha_i= " <<  2 * sum_alpha_i_beta_ij << endl;
		//cout << "sum_alpha_i_beta_ij = " << sum_alpha_i_beta_ij << endl;
		int ones = 0;
		int frac = 0;
		cout << "-------Les variables y non nulles " << endl;
		for (int k = 0; k < solution.size(); k++) {
			if (solution[k] > ::EPSILON) {
				if (solution.at(k) > 1.0 - EPSILON) {
					ones++;
				}
				else {
					frac++;
				}
				cout << "exchange k = " << k << " {" << columns[k].toString() << "} : " << solution[k] << "\n\n";
			}
		}
		cout << "COL " << solution.size() <<  "\tONES " << ones << "\tFRAC " << frac << endl;
		cout << endl;

	}
	catch (GRBException e) {
		cout << "RestrictedMasterProblem::write_solution_value()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "RestrictedMasterProblem::write_solution_value()] - Error during optimization" << endl;
	}
}

string RestrictedMasterProblem::write_solution_value_in_string() {
	try {
		stringstream mystream;
		mystream << "Objective   value: " << setprecision(10) << objective_value << endl;
		//double c_k_part = 0;
		//double sum_alpha_i_beta_ij = 0;

		//for (int i = 0; i < instance.size(); i++)
		//{
		//	sum_alpha_i_beta_ij += dual_variables.get_alpha_i(i);
		//	//for (int j = i + 1; j < instance.size(); j++)
		//	//	sum_alpha_i_beta_ij += dual_variables.get_beta_ij(i, j);
		//}
		//cout << " 2 * sum_alpha_i= " <<  2 * sum_alpha_i_beta_ij << endl;
		//cout << "sum_alpha_i_beta_ij = " << sum_alpha_i_beta_ij << endl;

		mystream << "-------Les variables y non nulles " << endl;
		for (int k = 0; k < solution.size(); k++)
			if (solution[k] > 0)
			{

				mystream << "exchange k = " << k << " {" << columns[k].toString() << "} : " << solution[k] << "\n\n";
			}
		mystream << endl;

		return mystream.str();

	}
	catch (GRBException e) {
		cout << "RestrictedMasterProblem::write_solution_value()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "RestrictedMasterProblem::write_solution_value()] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::write_solution_value_into_text(const int& iter) {
	ofstream myfile(::instanceName + "_RMP_iter_" + to_string(iter) + ".txt");
	if (myfile.is_open())
	{
		myfile << "Objective   value : " << setprecision(6) << objective_value << endl;
		myfile << "-------Les variables y non nulles " << endl;
		for (int k = 0; k < solution.size(); k++)
			if (solution[k] > 0)
			{
				myfile << "exchange k = " << k << " {" << columns[k].toString() << "} : " << solution[k] << endl;
			}

		myfile << endl;

		myfile.close();
	}
	else
		cout << "Unable to open file";
}

void RestrictedMasterProblem::write_solution_info_into_text(const int& file_name) {

	vector<double> packing(this->instance.getVertices().size(), 0.0);
	vector<vector<double>> arcsUsage(this->instance.getVertices().size(), vector<double>(this->instance.getVertices().size(), 0.0));
	for (int k = 0; k < this->columns.size(); k++) {
		double val = this->solution.at(k);
		if (val > EPSILON) {
			const Exchange* ex = &this->columns.at(k);
			for (int j = 1; j < ex->getSequence().size(); j++) {
				packing.at(ex->getSequence().at(j)) += val;
			}
			for (int j = 0; j < ex->getSequence().size() - 1; j++) {
				arcsUsage.at(ex->getSequence().at(j)).at(ex->getSequence().at(j + 1)) += val;
			}
		}
	}

	ofstream myfile(::instanceName + "_RMP_iter_" + to_string(file_name) + ".txt");
	if (myfile.is_open()) {
		myfile << "Objective   value : " << setprecision(6) << objective_value << endl;
		myfile << "-------Les variables y non nulles " << endl;
		for (int k = 0; k < solution.size(); k++) {
			if (solution[k] > 0) {
				myfile << "exchange k = " << k << " {" << columns[k].toString() << "} : " << solution[k] << endl;
			}
		}
		myfile << "\n\nVERTEX VAL:\n";
		for (int j = 0; j < packing.size(); j++) {
			if (packing.at(j) > EPSILON) {
				myfile << this->instance.getVertices().at(j).toString(false) << "\t" << packing.at(j) << "\n";
			}
		}
		myfile << "\n\nARCS VAL:\n";
		for (int j = 0; j < arcsUsage.size(); j++) {
			for (int i = 0; i < arcsUsage.at(j).size(); i++) {
				if (arcsUsage.at(j).at(i) > EPSILON) {
					myfile << "ARC " << j << "\t" << i << "\t" << arcsUsage.at(j).at(i) << "\n";
				}
			}
		}
		myfile << "\n\nCONFLICTS:\n";
		for (int k = 0; k < this->columns.size() - 1; k++) {
			if (solution[k] > 0) {
				for (int h = k + 1; h < this->columns.size(); h++) {
					if (solution[h] > 0 && this->columns.at(k).conflicts_with(this->columns.at(h))) {
						myfile << "exchange k = " << k << " and " << h << " are in conflict\n";
					}
				}
			}
		}

		myfile << endl;

		myfile.close();
	}
	else
		cout << "Unable to open file";

}

void RestrictedMasterProblem::write_dual_solution_value() {
	try {
		cout << "Dual variables alpha:" << endl;
		for (int i = 0; i < this->dualSolutionPi.size(); i++)
			if (fabs(this->dualSolutionPi.at(i)) > EPSILON)
				cout << "alpha_" << i << " = " << this->dualSolutionPi.at(i) << endl;

		cout << endl;

	}
	catch (GRBException e) {
		cout << "[RestrictedMasterProblem::write_dual_solution_value()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[RestrictedMasterProblem::write_dual_solution_value()] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::insert_columns(const vector<Exchange>& new_columns) {
	vector<Exchange>::iterator found;
	for (int i = 0; i < new_columns.size(); i++)
		insert_column(new_columns[i]);

	//this->export_model();
}

void RestrictedMasterProblem::insert_column(const Exchange& column) {
	try {

		this->columns.push_back(column);
		this->number_columns = static_cast<int>(this->columns.size());

		//cout << column.toString() << "\n";

		//Add group k to the constraints
		GRBColumn col;

		for (int i = 0; i < this->instance.getVertices().size(); i++) {
			const Vertex* v = &this->instance.getVertices().at(i);
			if (v->getType() != 'D') {
				col.addTerm(column.getVisit(v->getId()), this->packingConstraints.at(v->getId()));
			}
		}

		//for (int i = 0; i < this->packingConstraints.size(); i++) {
		//	//cout << this->packingConstraints.at(i).get(GRB_StringAttr_ConstrName) << "\t" << column.getVisit(i + 1) << "\n";
		//	col.addTerm(column.getVisit(i + 1), this->packingConstraints.at(i));
		//}
		col.addTerm(1, this->upperBoundExchangeNumber);

		for (int i = 0; i < this->cut_constraints.size(); i++) {
			if (this->cut_constraints.at(i)->get_name() == "SUBSETROW") {
				double coeff = 0.0;
				const SubsetRowConstraint<GRBConstr>* ssi = dynamic_cast<SubsetRowConstraint<GRBConstr>*>(this->cut_constraints.at(i));
				//cout << i<< " "<< ssi->toString() << endl;
				for (int j = 0; j < ssi->get_subset().size(); j++) {
					int id = ssi->get_subset().at(j);
					const Vertex* v = &this->instance.getVertices().at(id);
					coeff += (column.getVisit(id) * ssi->get_multipliers().at(j));
				}
				coeff /= ssi->getGCD();
				coeff += ::EPSILON;
				coeff = floor(coeff);
				//cout << "\tCOEFF " << coeff << endl;
				col.addTerm(coeff, *this->cut_constraints.at(i)->get_real_constraint());
			}
		}

		for (int i = 0; i < this->branching_constraints.size(); i++) {
			if (this->branching_constraints.at(i)->get_name() == "ARCBRANCH") {
				int id1 = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_items().first;
				int id2 = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_items().second;
				//cout << "ARC " << id1 << "\t" << id2 << "\t" << column.getTraverses(id1, id2) << endl;
				col.addTerm(column.getTraverses(id1, id2), *this->branching_constraints.at(i)->get_real_constraint());
			}
			else if (this->branching_constraints.at(i)->get_name() == "VISITBRANCH") {
				int id = dynamic_cast<VisitBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_item();
				col.addTerm(column.getVisit(id), *this->branching_constraints.at(i)->get_real_constraint());
			}
			else if (this->branching_constraints.at(i)->get_name() == "EXCHANGENUMBERBRANCH") {
				col.addTerm(1, *this->branching_constraints.at(i)->get_real_constraint());
			}
		}

		this->lambda.push_back(this->model->addVar(0, GRB_INFINITY, column.getWeight(), GRB_CONTINUOUS, col, "lambdaPP_" + to_string(this->lambda.size())));
		
	}
	catch (GRBException e) {
		cout << "[RestrictedMasterProblem::insert_column()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[RestrictedMasterProblem::extract_solution_value()] - Error during optimization" << endl;
	}

}

bool RestrictedMasterProblem::check_solution() const {
	bool isSolutionFeasible = true;
	vector<int> positiveColumns;
	double obj = 0.0;
	for (int k = 0; k < this->columns.size(); k++) {
		if (this->solution.at(k) > ::EPSILON) {
			positiveColumns.push_back(k);
			obj += this->columns.at(k).getWeight() * this->solution.at(k);
		}
	}
	//obj
	if (!::equal(obj, this->objective_value)) {
		cout << "WRONG OBJ " << obj << "\tSTORED\t" << this->objective_value << "\n";
		isSolutionFeasible = false;
	}
	//dual bound
	double sum = 0.0;
	for (int k = 0; k < positiveColumns.size(); k++) {
		sum += this->solution.at(positiveColumns.at(k));
	}
	if (::strictly_greater(sum, this->upperBoundExchangeNumber.get(GRB_DoubleAttr_RHS))) {
		cout << "WRONG UB EXCHANGE NUMBER " << sum << "\tRHS\t" << this->upperBoundExchangeNumber.get(GRB_DoubleAttr_RHS) << "\n";
		isSolutionFeasible = false;
	}
	//packing
	for (int i = 0; i < this->instance.getVertices().size(); i++) {
		const Vertex* v = &this->instance.getVertices().at(i);
		if (v->getType() != 'D') {
			sum = 0.0;
			for (int k = 0; k < positiveColumns.size(); k++) {
				sum += this->columns.at(positiveColumns.at(k)).getVisit(v->getId()) * this->solution.at(positiveColumns.at(k));
			}
			if (::strictly_greater(sum, 1.0)) {
				cout << "WRONG PACKING " << v->getId() << "\tLHS " << sum << "\n";
				isSolutionFeasible = false;
			}
		}
	}
	//branch
	for (int i = 0; i < this->branching_constraints.size(); i++) {
		double rhs = this->branching_constraints.at(i)->get_rhs();
		sum = 0.0;
		if (this->branching_constraints.at(i)->get_name() == "ARCBRANCH") {
			int id1 = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_items().first;
			int id2 = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_items().second;
			for (int k = 0; k < positiveColumns.size(); k++) {
				sum += this->columns.at(positiveColumns.at(k)).getTraverses(id1,id2) * this->solution.at(positiveColumns.at(k));
			}
		}
		else if (this->branching_constraints.at(i)->get_name() == "VISITBRANCH") {
			int id = dynamic_cast<VisitBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i))->get_item();
			for (int k = 0; k < positiveColumns.size(); k++) {
				sum += this->columns.at(positiveColumns.at(k)).getVisit(id) * this->solution.at(positiveColumns.at(k));
			}
		}
		else if (this->branching_constraints.at(i)->get_name() == "EXCHANGENUMBERBRANCH") {
			for (int k = 0; k < positiveColumns.size(); k++) {
				sum += this->solution.at(positiveColumns.at(k));
			}
		}
		if (this->branching_constraints.at(i)->get_operator() == "<=") {
			if (::strictly_greater(sum, rhs)) {
				cout << "WRONG BRANCHING " << i << "\tLHS " << sum << "\tRHS" << rhs << "\n";
				isSolutionFeasible = false;
			}
		}
		else if (this->branching_constraints.at(i)->get_operator() == ">=") {
			if (::strictly_less(sum, rhs)) {
				cout << "WRONG BRANCHING " << i << "\tLHS " << sum << "\tRHS" << rhs << "\n";
				isSolutionFeasible = false;
			}
		}
		else {
			if (!::equal(sum, rhs)) {
				cout << "WRONG BRANCHING " << i << "\tLHS " << sum << "\tRHS" << rhs << "\n";
				isSolutionFeasible = false;
			}
		}
	}
	//subset row
	for (int i = 0; i < this->cut_constraints.size(); i++) {
		if (this->cut_constraints.at(i)->get_name() == "SUBSETROW") {
			const SubsetRowConstraint<GRBConstr>* ssi = dynamic_cast<SubsetRowConstraint<GRBConstr>*>(this->cut_constraints.at(i));
			double lhs = 0.0;
			double rhs = ssi->get_rhs();
			for (int k = 0; k < positiveColumns.size(); k++) {
				double colContrib = 0.0;
				for (int l = 0; l < ssi->get_subset().size(); l++) {
					int vId = ssi->get_subset().at(l);
					colContrib += this->columns.at(positiveColumns.at(k)).getVisit(vId) * ssi->get_multipliers().at(l);
				}
				colContrib /= static_cast<double>(ssi->getGCD());
				colContrib += ::EPSILON;
				colContrib = floor(colContrib);
				lhs += (colContrib * this->solution.at(positiveColumns.at(k)));
			}
			if (::strictly_greater(lhs, rhs)) {
				cout << "WRONG SSI " << i << "\tLHS " << setprecision(6) << lhs << "\tRHS " << setprecision(6) << rhs << "\n";
				isSolutionFeasible = false;
			}
		}
	}

	return isSolutionFeasible;

}

bool RestrictedMasterProblem::check_opt_sol_violation() const {
	bool isSolutionFeasible = true;

	vector<Exchange> optSol;
	optSol.push_back(Exchange("chain", 873, 12, { 0,4,154,259,248,128,188,69,181,210,180,95,42 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 809, 13, { 0,6,84,166,174,263,30,91,144,108,171,113,203,54 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 573, 9, { 0,5,196,256,151,125,41,82,37,215 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 762, 12, { 0,1,167,236,160,157,62,67,26,194,70,23,60 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 748, 13, { 0,8,29,100,161,88,205,75,233,96,250,14,18,118 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 715, 11, { 0,13,172,121,229,226,99,245,146,184,204,55 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 814, 12, { 0,12,105,45,198,119,78,111,260,230,123,156,223 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 878, 13, { 0,10,104,202,56,57,76,72,68,87,131,197,253,138 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 788, 13, { 0,11,24,222,234,59,243,206,241,77,155,189,140,249 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 911, 13, { 0,9,136,94,64,133,135,150,231,109,43,255,195,147 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 780, 12, { 0,7,110,116,220,15,182,169,132,79,63,227,107 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 797, 13, { 0,2,46,61,217,103,170,200,25,152,240,74,232,237 }, this->instance.getVertexNumber()));
	optSol.push_back(Exchange("chain", 682, 11, { 0,3,228,106,201,173,27,143,153,31,81,51 }, this->instance.getVertexNumber()));

	//dual bound
	double sum = 0.0;
	for (int k = 0; k < optSol.size(); k++) {
		//cout << optSol.at(k).toString() << endl;
		sum += 1.0;
	}
	if (::strictly_greater(sum, this->upperBoundExchangeNumber.get(GRB_DoubleAttr_RHS))) {
		cout << "WRONG UB EXCHANGE NUMBER " << sum << "\tRHS\t" << this->upperBoundExchangeNumber.get(GRB_DoubleAttr_RHS) << "\n";
		isSolutionFeasible = false;
	}
	//packing
	for (int i = 0; i < this->instance.getVertices().size(); i++) {
		const Vertex* v = &this->instance.getVertices().at(i);
		if (v->getType() != 'D') {
			sum = 0.0;
			for (int k = 0; k < optSol.size(); k++) {
				sum += optSol.at(k).getVisit(v->getId());
			}
			if (::strictly_greater(sum, 1.0)) {
				cout << "WRONG PACKING " << v->getId() << "\tLHS " << sum << "\n";
				isSolutionFeasible = false;
			}
		}
	}
	
	//subset row
	for (int i = 0; i < this->cut_constraints.size(); i++) {
		if (this->cut_constraints.at(i)->get_name() == "SUBSETROW") {
			const SubsetRowConstraint<GRBConstr>* ssi = dynamic_cast<SubsetRowConstraint<GRBConstr>*>(this->cut_constraints.at(i));
			double lhs = 0.0;
			double rhs = ssi->get_rhs();
			for (int k = 0; k < optSol.size(); k++) {
				double colContrib = 0.0;
				for (int l = 0; l < ssi->get_subset().size(); l++) {
					int vId = ssi->get_subset().at(l);
					colContrib += optSol.at(k).getVisit(vId) * ssi->get_multipliers().at(l);
				}
				colContrib /= static_cast<double>(ssi->getGCD());
				colContrib += ::EPSILON;
				colContrib = floor(colContrib);
				lhs += (colContrib);
			}
			if (::strictly_greater(lhs, rhs)) {
				cout << "WRONG SSI " << i << "\tLHS " << setprecision(6) << lhs << "\tRHS " << setprecision(6) << rhs << "\n";
				isSolutionFeasible = false;
			}
		}
	}

	return isSolutionFeasible;

}

void RestrictedMasterProblem::change_variable_type_to_binary() {
	for (int k = 0; k < this->lambda.size(); k++) {
		this->lambda.at(k).set(GRB_CharAttr_VType, GRB_BINARY);
	}
}

void RestrictedMasterProblem::change_variable_type_to_continuous() {
	for (int k = 0; k < this->lambda.size(); k++) {
		this->lambda.at(k).set(GRB_CharAttr_VType, GRB_CONTINUOUS);
	}
}

void RestrictedMasterProblem::attach(RMPNodeObserver* observer) {
	this->list_observers.push_back(observer);
}

void RestrictedMasterProblem::detach(RMPNodeObserver* observer) {
	this->list_observers.remove(observer);
}

//void RestrictedMasterProblem::notify_nodes_arc_branching_constraints(const vector<pair<int, pair< int, int>>>& new_column_indexes_and_involved_items) {
//	list<RMPNodeObserver*>::iterator iterator = this->list_observers.begin();
//	while (iterator != this->list_observers.end()) {
//		(*iterator)->update_arc_branching_constraints(new_column_indexes_and_involved_items);
//		++iterator;
//	}
//}

void RestrictedMasterProblem::notify_nodes_clique_inequalities(const vector<int>& new_column_indexes) {
	list<RMPNodeObserver*>::iterator iterator = this->list_observers.begin();
	while (iterator != this->list_observers.end()) {
		(*iterator)->update_clique_inequalities(new_column_indexes);
		++iterator;
	}
}

void RestrictedMasterProblem::update_ineffective_cuts() {
	/*double sum;
	vector<Constraint<GRBConstr>*>::iterator cut = this->cut_constraints.begin();
	while (cut != this->cut_constraints.end()) {
		sum = 0;
		for (auto& k : (*cut)->get_variables_indexes())
			sum += get_solution()[k];

		if ((*cut)->get_rhs() - sum > INEFFICACITYEPSILON) {
			if ((*cut)->get_become_ineffective() == 0)
				(*cut)->set_become_ineffective(this->solving_rmp_iteration);
			(*cut)->set_ineffective_counter((*cut)->get_ineffective_counter() + 1);
			if ((*cut)->get_ineffective_counter() >= MAX_NB_ITER_INEFFECTIVE_CUT) {
				add_ineffective_cut_constraint(*cut);
				this->model->remove(*(*cut)->get_real_constraint());
				cut = this->cut_constraints.erase(cut);
			}
			else
				cut++;
		}
		else
			cut++;
	}*/
}

void RestrictedMasterProblem::attach(RMPCutInformationObserver* observer) {
	this->list_cut_family_information_observer.push_back(observer);
}

void RestrictedMasterProblem::detach(RMPCutInformationObserver* observer) {
	this->list_cut_family_information_observer.remove(observer);
}

void RestrictedMasterProblem::attach(RMPBaPObserver* observer) {
	this->bapObserver = observer;
}

void RestrictedMasterProblem::detach_BaP_observer() {
	this->bapObserver = nullptr;
}

void RestrictedMasterProblem::notify_cut_information_observer() {
	list<RMPCutInformationObserver*>::iterator iterator = this->list_cut_family_information_observer.begin();
	vector<double> information;
	information.push_back(get_objective_value());
	information.push_back(are_all_variables_binary() ? 1 : 0);
	information.push_back(get_number_new_columns());
	information.push_back(get_solving_time());
	while (iterator != this->list_cut_family_information_observer.end()) {
		(*iterator)->update_information(information);
		iterator = this->list_cut_family_information_observer.erase(iterator);
	}
}

bool RestrictedMasterProblem::detect_conflicting_branching_rules() const {
	for (int i = 0; i < this->branching_constraints.size(); i++) {
		if (this->branching_constraints[i]->get_name() == "ARCBRANCH") {
			for (int j = i + 1; j < this->branching_constraints.size(); j++) {
				if (this->branching_constraints[i]->get_name() == "ARCBRANCH") {
					if (dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints[i])->get_items() ==
						dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints[j])->get_items() &&
						((this->branching_constraints[i]->get_operator() == "<=" && this->branching_constraints[j]->get_operator() == ">=") ||
						 (this->branching_constraints[i]->get_operator() == ">=" && this->branching_constraints[j]->get_operator() == "<="))
						)
						return true;
				}
			}
		}
		else if (this->branching_constraints[i]->get_name() == "VISITBRANCH") {
			for (int j = i + 1; j < this->branching_constraints.size(); j++) {
				if (this->branching_constraints[i]->get_name() == "VISITBRANCH") {
					if (dynamic_cast<VisitBranchConstraint<GRBConstr>*>(this->branching_constraints[i])->get_item() ==
						dynamic_cast<VisitBranchConstraint<GRBConstr>*>(this->branching_constraints[j])->get_item() &&
						((this->branching_constraints[i]->get_operator() == "<=" && this->branching_constraints[j]->get_operator() == ">=") ||
							(this->branching_constraints[i]->get_operator() == ">=" && this->branching_constraints[j]->get_operator() == "<="))
						)
						return true;
				}
			}
		}
	}
	return false;
}

void RestrictedMasterProblem::notify_nodes_for_ineffective_cuts()
{
	list<RMPNodeObserver*>::iterator iterator = list_observers.begin();
	while (iterator != list_observers.end()) {
		(*iterator)->update_ineffective_cuts();
		++iterator;
	}
}

void RestrictedMasterProblem::notify_BaP_for_new_incumbent() {
	this->bapObserver->update_best_incumbent();
}

void RestrictedMasterProblem::addLambdaVariables() {
	this->lambda.clear();
	for (int k = 0; k < this->columns.size(); k++) {
		string name = "lambda_" + to_string(k);
		this->lambda.push_back(this->model->addVar(0.0, GRB_INFINITY, 0.0, GRB_CONTINUOUS, name));
	}
}

void RestrictedMasterProblem::call_restricted_master_heuristic_new_MIP(const double& timeLimit, const unordered_set<int>& warmStart, const bool& fix_variables) {

	try {

		if (this->lambda.empty()) {
			return;
		}

		GRBEnv* envMIP = new GRBEnv();
		GRBModel* modelMIP = new GRBModel(*envMIP);

		/****** Variables ******/
		vector<GRBVar> lambdas;
		unordered_map<int, int> dictionary;
		for (int i = 0; i < this->lambda.size(); i++) {
			string name = "lambda_" + to_string(i);
			lambdas.push_back(modelMIP->addVar(0.0, 1.0, 0.0, GRB_INTEGER, name));
		}

		/***** Constraints *****/
		vector<GRBConstr> packingConstraints(this->instance.getVertexNumber());
		size_t size = this->columns.size();
		for (int i = 0; i < this->instance.getVertexNumber(); i++) {
			const Vertex* v = &this->instance.getVertices().at(i);
			int id = v->getId();
			if (v->getType() != 'D') {
				GRBLinExpr sum = 0;
				//vector<double> coeffs(this->columns.size());
				for (int k = 0; k < size; k++) {
					//sum += this->columns.at(k).getVisit(v->getId()) * this->lambda.at(k);
					//coeffs.at(k) = this->columns.at(k).getVisit(id);
					double coeff = this->columns.at(k).getVisit(i);
					if (coeff > ::EPSILON) {
						sum.addTerms(&coeff, &lambdas.at(k), 1);
					}
				}
				//sum.addTerms(&coeffs[0], &this->lambda[0], size);
				packingConstraints.at(id) = modelMIP->addConstr(sum <= 1, "packing_" + to_string(id));
			}
		}
		if (fix_variables && !warmStart.empty()) {
			for (int k = 0; k < this->columns.size(); k++) {
				int id = this->columns.at(k).getId();
				if (warmStart.find(id) != warmStart.end()) {
					modelMIP->addConstr(lambdas.at(k) == 1, "fix_exchange_" + to_string(k));
				}
				else {
					modelMIP->addConstr(lambdas.at(k) == 0, "fix_exchange_" + to_string(k));
				}
			}
		}

		/****** Objective function ******/
		GRBLinExpr objective = 0;
		for (int k = 0; k < this->columns.size(); k++) {
			objective += this->columns.at(k).getWeight() * lambdas.at(k);
		}

		modelMIP->setObjective(objective, GRB_MAXIMIZE);

		/********Warm-start*************/
		if (!warmStart.empty()) {
			for (int k = 0; k < this->columns.size(); k++) {
				int id = this->columns.at(k).getId();
				if (warmStart.find(id) != warmStart.end()) {
					lambdas.at(k).set(GRB_DoubleAttr_Start, 1);
				}
			}
		}

		/****** Export model ******/
		//export_model();

		//modelMIP->set(GRB_DoubleParam_TimeLimit, timeLimit);
		modelMIP->set(GRB_DoubleParam_TimeLimit, timeLimit);
		modelMIP->set(GRB_IntParam_LogToConsole, 0);
		modelMIP->set(GRB_IntParam_Threads, 1);
		modelMIP->set(GRB_IntParam_MIPFocus, 1);
		modelMIP->set(GRB_DoubleParam_MemLimit, 16.0);
		//modelMIP->set(GRB_DoubleParam_Heuristics, 1);
		//this->export_model();
		modelMIP->optimize();

		int status = modelMIP->get(GRB_IntAttr_Status);
		int solCount = modelMIP->get(GRB_IntAttr_SolCount);
		::time_final_MIP = modelMIP->get(GRB_DoubleAttr_Runtime);
		::gap_final_MIP = -1.0;
		
		if (solCount > 0) {

			::gap_final_MIP = modelMIP->get(GRB_DoubleAttr_MIPGap);
			double candidate_value = modelMIP->get(GRB_DoubleAttr_ObjVal);

			if (::strictly_greater(candidate_value, this->objective_value_int)) {
				this->objective_value_int = candidate_value;

				this->solution_columns_int.clear();
				for (int k = 0; k < this->columns.size(); k++) {
					if (lambdas.at(k).get(GRB_DoubleAttr_X) > 0.5) {
						this->solution_columns_int.insert(k);
					}
				}

			}

		}
		delete modelMIP;
		delete envMIP;

	}
	catch (GRBException e) {
		cout << "[call_restricted_master_heuristic_new_MIP()] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[call_restricted_master_heuristic_new_MIP()] - Error during optimization" << endl;
	}

}

string RestrictedMasterProblem::to_string_integer_solution() {
	string res = "OBJECTIVE VALUE:\t" + to_string(this->objective_value_int) + "\n";
	res.append("COLUMNS:\n");
	for (const auto& col : this->solution_columns_int) {
		res.append("\t" + this->columns.at(col).toString() + "\n\n");
	}
	return res;
}

double RestrictedMasterProblem::compute_reduced_cost(const Exchange& e) {
	double rd = e.getWeight() - this->dualSolutionBeta;
	for (int i = 0; i < e.getSequence().size(); i++) {
		rd -= this->dualSolutionPi.at(e.getSequence().at(i));
	}
	for (int i = 0; i < this->branching_constraints.size(); i++) {
		double coeff = -1;
		if (this->branching_constraints.at(i)->get_name() == "ARCBRANCH") {
			const ArcBranchConstraint<GRBConstr>* ab = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i));
			int id1 = ab->get_items().first;
			int id2 = ab->get_items().second;
			coeff = e.getTraverses(id1, id2);
		}
		else if (this->branching_constraints.at(i)->get_name() == "VISITBRANCH") {
			const VisitBranchConstraint<GRBConstr>* vb = dynamic_cast<VisitBranchConstraint<GRBConstr>*>(this->branching_constraints.at(i));
			int id = vb->get_item();
			coeff = e.getVisit(id);
		}
		else if (this->branching_constraints.at(i)->get_name() == "EXCHANGENUMBERBRANCH") {
			coeff = 1;
		}
		rd -= (coeff * this->dualSolutionGamma.at(i));
	}
	for (int i = 0; i < this->cut_constraints.size(); i++) {
		double coeff = -1;
		if (this->cut_constraints.at(i)->get_name() == "SUBSETROW") {
			const SubsetRowConstraint<GRBConstr>* ssi = dynamic_cast<SubsetRowConstraint<GRBConstr>*>(this->cut_constraints.at(i));
			coeff = 0.0;
			for (int j = 0; j < ssi->get_subset().size(); j++) {
				coeff += (ssi->get_multipliers().at(j) * e.getVisit(ssi->get_subset().at(j)));
			}
			coeff /= ssi->getGCD();
			coeff += EPSILON;
			coeff = floor(coeff);
		}
		rd -= (coeff * this->dualSolutionSigma.at(i));
	}
	return rd;
}

double RestrictedMasterProblem::compute_lagrangian_bound_RMP_obj_term() const {
	double term = 0.0;
	for (int i = 0; i < this->dualSolutionPi.size(); i++) {
		term += this->dualSolutionPi.at(i);
	}
	for (int i = 0; i < this->dualSolutionGamma.size(); i++) {
		term += (this->dualSolutionGamma.at(i) * this->branching_constraints.at(i)->get_rhs());
	}
	for (int i = 0; i < this->dualSolutionSigma.size(); i++) {
		term += (this->dualSolutionSigma.at(i) * this->cut_constraints.at(i)->get_rhs());
	}
	return term;
}

bool RestrictedMasterProblem::is_infeasible_dummy_variable() const {
	return false;
}

void RestrictedMasterProblem::addPackingConstraints() {
	//Every subject appears in two columns.
	this->packingConstraints = vector<GRBConstr>(this->instance.getVertexNumber());
	size_t size = this->columns.size();
	for (int i = 0; i < this->instance.getVertexNumber(); i++) {
		const Vertex* v = &this->instance.getVertices().at(i);
		int id = v->getId();
		if (v->getType() != 'D') {
			GRBLinExpr sum = 0;
			//vector<double> coeffs(this->columns.size());
			for (int k = 0; k < size; k++) {
				//sum += this->columns.at(k).getVisit(v->getId()) * this->lambda.at(k);
				//coeffs.at(k) = this->columns.at(k).getVisit(id);
				double coeff = this->columns.at(k).getVisit(i);
				if (coeff > ::EPSILON) {
					sum.addTerms(&coeff, &this->lambda.at(k), 1);
				}
			}
			//sum.addTerms(&coeffs[0], &this->lambda[0], size);
			this->packingConstraints.at(id) = this->model->addConstr(sum <= 1, "packing_" + to_string(id));
		}
	}
}

void RestrictedMasterProblem::addUpperBoundExchangeNumber() {
	GRBLinExpr sum = 0;
	for (int k = 0; k < this->columns.size(); k++)
		sum += this->lambda.at(k);
	this->upperBoundExchanges = static_cast<double>(this->instance.getVertexNumber() - 1) / 2.0;
	this->upperBoundExchangeNumber = this->model->addConstr(sum <= this->upperBoundExchanges, "UBexchanges");
}

void RestrictedMasterProblem::objective_function() {
	this->objective = 0;
	for (int k = 0; k < this->columns.size(); k++)
		this->objective += this->columns.at(k).getWeight() * this->lambda.at(k);

	this->model->setObjective(objective, GRB_MAXIMIZE);
}

void RestrictedMasterProblem::free_memory() {
	delete model;
	delete env;
}

void RestrictedMasterProblem::export_model() {
	this->model->write("RestrictedMasterModel.lp");
}

void RestrictedMasterProblem::export_model(const int& id) {
	this->model->write("RestrictedMasterModel" + to_string(id) + ".lp");
}

bool RestrictedMasterProblem::are_all_variables_binary() const {
	//double int_part;
	//cout << "are_all_variables_binary\n";
	for (vector<double>::const_iterator value = this->solution.begin(); value != this->solution.end(); value++) {
		//cout << "VALUE " << *value << "\t" << modf(*value, &int_part) << endl;
		if (*value > ::EPSILON) {
			double fractionalPart = computeFractionalPart(*value);
			if (fractionalPart > ::EPSILON) {
				return false;
			}
		}
	}

	return true;
}

void RestrictedMasterProblem::update_integer_solution() {

	if (::strictly_greater(this->objective_value, this->objective_value_int)) {
		this->objective_value_int = this->objective_value;
		this->solution_columns_int.clear();
		for (int i = 0; i < this->solution.size(); i++) {
			if (this->solution.at(i) > 0.5) {
				this->solution_columns_int.insert(i);
			}
		}
	}

}

void RestrictedMasterProblem::sanity_check_integer_solution() const {
	unordered_map<int, int> checkDuplicate;
	int solVal = 0;
	for (const auto& col_id : this->solution_columns_int) {
		solVal += this->columns.at(col_id).getWeight();
		for (int k = 1; k < this->columns.at(col_id).getSequence().size(); k++) {
			int id = this->columns.at(col_id).getSequence().at(k);
			if (checkDuplicate.find(id) == checkDuplicate.end()) {
				checkDuplicate.insert(make_pair(id, 1));
			}
			else {
				checkDuplicate.at(id)++;
			}
		}
	}

	bool isFail = false;
	if (!::equal(round(solVal), round(this->objective_value_int))) {
		cout << "\n\n\nRECOMPUTED " << solVal << " STORED " << this->objective_value_int << " !!!\n\n\n";
		isFail = true;
	}
	for (unordered_map<int, int>::const_iterator it = checkDuplicate.begin(); it != checkDuplicate.end(); it++) {
		if (it->second > 1) {
			cout << "\n\n\nVERTEX " << it->first << " VISITED " << it->second << " TIMES!!!\n\n\n";
			isFail = true;
		}
	}
	if (isFail) {
		exit(EXIT_FAILURE);
	}
}

void RestrictedMasterProblem::set_warm_start(const unordered_set<int>& warmStart) {
	for (int k = 0; k < this->columns.size(); k++) {
		int id = this->columns.at(k).getId();
		if (warmStart.find(id) != warmStart.end()) {
			this->lambda.at(k).set(GRB_DoubleAttr_Start, 1);
		}
	}
}

bool RestrictedMasterProblem::is_infeasible() const {
	int optimstatus = this->model->get(GRB_IntAttr_Status);
	int solCount = this->model->get(GRB_IntAttr_SolCount);
	return (optimstatus != GRB_OPTIMAL);
}

bool RestrictedMasterProblem::is_feasible_MIP() const {
	int solCount = this->model->get(GRB_IntAttr_SolCount);
	return (solCount > 0);
}

int RestrictedMasterProblem::get_model_status() const {
	return this->model->get(GRB_IntAttr_Status);
}

double RestrictedMasterProblem::get_model_MIP_gap() const {
	return this->model->get(GRB_DoubleAttr_MIPGap);
}

void RestrictedMasterProblem::initialize() {
	RestrictedMasterProblemInterface<Instance, vector<double>, DualSolution, vector<Exchange>, GRBConstr>::initialize();
	if (this->initializer->get_cost() < this->objective_value)
		this->objective_value = this->initializer->get_cost();
	set_initial_solution_value(this->objective_value);
}

void RestrictedMasterProblem::add_branching_constraint(Constraint<GRBConstr>* new_constraint, const int& pos) {
	try {
		GRBLinExpr sum;
		//build the left hand-side
		string name = "branch_" + to_string(pos) + "_";
		if (new_constraint->get_name() == "ARCBRANCH") {
			const ArcBranchConstraint<GRBConstr>* ab = dynamic_cast<ArcBranchConstraint<GRBConstr>*>(new_constraint);
			int id1 = ab->get_items().first;
			int id2 = ab->get_items().second;
			name += "arc_" + to_string(id1) + "_" + to_string(id2);
			for (int i = 0; i < this->columns.size(); i++) {
				int coeff = this->columns.at(i).getTraverses(id1, id2);
				sum += coeff * this->lambda.at(i);
			}
		}
		else if (new_constraint->get_name() == "VISITBRANCH") {
			const VisitBranchConstraint<GRBConstr>* vb = dynamic_cast<VisitBranchConstraint<GRBConstr>*>(new_constraint);
			int id = vb->get_item();
			name += "visit_" + to_string(id);
			for (int i = 0; i < this->columns.size(); i++) {
				int coeff = this->columns.at(i).getVisit(id);
				sum += coeff * this->lambda.at(i);
			}
		}
		else if (new_constraint->get_name() == "EXCHANGENUMBERBRANCH") {
			name += "exchange";
			for (int i = 0; i < this->columns.size(); i++) {
				sum += this->lambda.at(i);
			}
		}
		GRBConstr constraint;
		if (new_constraint->get_operator() == "<=") {
			this->branching_constraints.push_back(new_constraint);
			new_constraint->set_real_constraint(new GRBConstr(this->model->addConstr(sum <= new_constraint->get_rhs(), name)));
		}
		else if (new_constraint->get_operator() == ">=") {
			this->branching_constraints.push_back(new_constraint);
			new_constraint->set_real_constraint(new GRBConstr(this->model->addConstr(sum >= new_constraint->get_rhs(), name)));
		}
		else {
			this->branching_constraints.push_back(new_constraint);
			new_constraint->set_real_constraint(new GRBConstr(this->model->addConstr(sum == new_constraint->get_rhs(), name)));
		}
	}
	catch (GRBException e) {
		cout << "[add_branching_constraint] - Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "[add_branching_constraint] - Error during optimization" << endl;
	}
}

void RestrictedMasterProblem::add_branching_constraints(vector<Constraint<GRBConstr>*>& new_branching_constraints) {
	for (int i = 0; i < new_branching_constraints.size(); i++) {
		add_branching_constraint(new_branching_constraints.at(i), i);
	}
	/*for (vector<Constraint<GRBConstr>*>::iterator branching = new_branching_constraints.begin(); branching != new_branching_constraints.end(); branching++)
		add_branching_constraint(*branching);*/
}

void RestrictedMasterProblem::add_cut_constraint(Constraint<GRBConstr>* new_constraint, const int& pos, const bool& isAddedInList) {
	GRBLinExpr sum;
	string name = "Cut" + to_string(pos) + "_";
	if (new_constraint->get_name() == "SUBSETROW") {
		const SubsetRowConstraint<GRBConstr>* ssi = dynamic_cast<SubsetRowConstraint<GRBConstr>*>(new_constraint);
		name += "ssi" + ssi->get_subsetToString();
		for (int i = 0; i < this->columns.size(); i++) {
			double coeff = 0.0;
			for (int j = 0; j < ssi->get_subset().size(); j++) {
				coeff += (ssi->get_multipliers().at(j) * this->columns.at(i).getVisit(ssi->get_subset().at(j)));
			}
			coeff /= ssi->getGCD();
			coeff += EPSILON;
			coeff = floor(coeff);
			sum += coeff * this->lambda.at(i);
		}
	}
	if (isAddedInList) {
		this->cut_constraints.push_back(new_constraint);
	}
	if (new_constraint->get_operator() == "<=") {
		new_constraint->set_real_constraint(new GRBConstr(this->model->addConstr(sum <= new_constraint->get_rhs(), name)));
	}
	else if (new_constraint->get_operator() == ">=") {
		new_constraint->set_real_constraint(new GRBConstr(this->model->addConstr(sum >= new_constraint->get_rhs(), name)));
	}
	else {
		new_constraint->set_real_constraint(new GRBConstr(this->model->addConstr(sum == new_constraint->get_rhs(), name)));
	}
	//export_model();
}

void RestrictedMasterProblem::add_cut_constraints(vector<Constraint<GRBConstr>*>& new_cut_constraints) {
	for (int i = 0; i < new_cut_constraints.size(); i++) {
		add_cut_constraint(new_cut_constraints.at(i), i, true);
	}
	/*for(vector<Constraint<GRBConstr>*>::iterator cut = new_cut_constraints.begin(); cut != new_cut_constraints.end(); cut++)
		add_cut_constraint(*cut);*/
}

void RestrictedMasterProblem::add_ineffective_cut_constraint(Constraint<GRBConstr>* new_ineffective_cut_constraint)
{
	ineffective_cut_constraints.push_back(new_ineffective_cut_constraint);
}

void RestrictedMasterProblem::remove_from_model_ineffective_cut_constraint(Constraint<GRBConstr>* new_ineffective_cut_constraint)
{
	model->remove(*new_ineffective_cut_constraint->get_real_constraint());
}

void RestrictedMasterProblem::remove_from_cut_constraints_ineffective_cut_constraint(Constraint<GRBConstr>* new_ineffective_cut_constraint)
{
	vector<Constraint<GRBConstr>*>::iterator cut = cut_constraints.begin();
	while (*cut != new_ineffective_cut_constraint && cut != cut_constraints.end())
		cut++;
	if (cut != cut_constraints.end())
		/*new_ineffective_cut_constraint = **/cut_constraints.erase(cut);
}

void RestrictedMasterProblem::remove_branching_constraints()
{
	while (this->branching_constraints.begin() != this->branching_constraints.end()) {
		this->model->remove(*this->branching_constraints.front()->get_real_constraint());
		this->branching_constraints.erase(this->branching_constraints.begin());
	}
}

void RestrictedMasterProblem::remove_cut_constraints() {
	while (this->cut_constraints.begin() != this->cut_constraints.end()) {
		this->model->remove(*this->cut_constraints.front()->get_real_constraint());
		this->cut_constraints.erase(this->cut_constraints.begin());
	}
}

void RestrictedMasterProblem::remove_cut_constraints_from_model() {
	for (int i = 0; i < this->cut_constraints.size(); i++) {
		this->model->remove(*this->cut_constraints.at(i)->get_real_constraint());
	}		
}

void RestrictedMasterProblem::add_cut_constraints_in_model() {
	for(int i = 0; i < this->cut_constraints.size(); i++) {
		add_cut_constraint(this->cut_constraints.at(i), i, false);
	}
}

//FUNCTIONS TO ASK
//bool RestrictedMasterProblem::add_auxiliary_variables_and_induced_linking_constraints() {
//	return false;
//}

void RestrictedMasterProblem::fake_solution() {

	this->columns.clear();
	this->solution.clear();

	this->solution = { 0.604211,0.275789,0.724211,0.458947,0.395789,0.671579,0.328421,0.275789,0.671579,0.395789,0.275789,0.328421,0.724211 };

	this->columns.push_back(Exchange("cycle", 183, 3, { 424,303,77,424 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 131, 3, { 393,122,179,393 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 221, 3, { 499,268,167,499 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 178, 3, { 323,122,179,323 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 192, 3, { 115,232,77,115 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 180, 3, { 549,345,215,549 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 239, 3, { 313,467,232,313 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 199, 3, { 268,232,215,268 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 154, 3, { 313,318,263,313 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("cycle", 200, 3, { 323,424,336,323 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("chain", 220, 4, { 0,37,342,167,507 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("chain", 183, 4, { 0,54,345,263,113 }, this->instance.getVertexNumber()));
	this->columns.push_back(Exchange("chain", 184, 4, { 0,37,393,241,464 }, this->instance.getVertexNumber()));

}