#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <iostream>
#include <gurobi_c++.h>

#ifdef VISUALC
#include <windows.h>
#else
#include <unistd.h>
#include <sys/times.h>
#endif

using namespace std;

/* CONSTANTS */

#define EPSILON 0.00001
#define VIOLATION_TRIINEQ 0.1
#define VIOLATION_PART2  0.1
#define VIOLATION_COVERINEQ  0.1
#define VIOLATION_EXHAUSTIVE_PERCENTAGERHS 0.01
#define MAXTIME 3600.0
#define BIGVALUE 9999.9
#define MINTENURE 5
#define MAXINFEASTENURE 40
#define MAXFEASTENURE 20
#define MAXCOVERTEST 7

/* ONLY FOR RESOLUTION OF A REAL-WORLD INSTANCE (CHOOSE AT MOST ONLY ONE) */

//#define INST_2014
//#define INST_2015
//#define INST_2015_PROPOSAL

/* GENERAL FEATURES (COMMENT IT FOR TURNING OFF) */

#define BETA_EQUALITY
#define USE_TABU_SEARCH
//#define PUREBYB
#define WITH_WEIGHTS
#define ENABLE_CUTS
#define PART2_CUTS
#define WEIGHTED_COVER_CUTS
//#define WEIGHTED_COVER_CUTS_EXHAUSTIVE
//#define WEIGHTED_UPPERBOUND_CUTS_EXHAUSTIVE
//#define WEIGHTED_LOWERBOUND_CUTS_EXHAUSTIVE
#define PERCENTAGE_TRIINEQ_IN_FORM 20
#define PERCENTAGE_TRIINEQ_IN_CUTS 70
//#define HISTOGRAM
//#define ADD_ARTIFICIAL_NODES

#if PERCENTAGE_TRIINEQ_IN_FORM < 100
#define TRIINEQ_LAZYS
#endif

#ifdef HISTOGRAM
int hist[10];
#endif

/* GLOBAL VARIABLES */

double start_t; /* initial timestamp */
int n, k; /* number of vertices and clusters */
float **d; /* distance matrix */
float *w; /* weights of vertices */
float WL, WU; /* weight bounds */
int FL, FU; /* size bounds */
int *best; /* best assignment so far */
float best_cost; /* cost of the best assignment */

struct triangular_ineq {
	int i, j, r;
	int sign_ij, sign_ir, sign_jr;
	float dist;
} *t_ineq;
int t_ineq_count, t_ineq_count_in_formulation, t_ineq_count_during_cuts; /* triangular combinations: "all", "inside formulation", "for cuts" */

struct order_elements_by_val {
	int v;
	float val;
};

/* FUNCTIONS */

/*
 * Tclock - obtain a timestamp (in seconds)
 */
double ECOclock() {
#ifdef VISUALC
	/* measure standard wall-clock: use it on Windows
	* (Linux clock() function overflows every 72 minutes) */
	return ((double)clock()) / (double)CLOCKS_PER_SEC;
#else
	/* measure user-time: use it on single-threaded system in Linux (is more accurate) */
	struct tms buf;
	times(&buf);
	return ((double)buf.tms_utime) / (double)sysconf(_SC_CLK_TCK);
#endif
}

/*
 * set_color - change font color
 */
void set_color(int color)
{
#ifdef VISUALC
	SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), color);
#endif
}

/*
 * bye - finalize execution and write a string on screen
 */
void bye(char *string)
{
	set_color(12);
	cout << string << endl;
	set_color(7);
	exit(1);
}

/* auxiliary function to sort triangular inequalities */
int cmpfunc_t_ineq(const void *a, const void *b)
{
	struct triangular_ineq *ta = (struct triangular_ineq*)a;
	struct triangular_ineq *tb = (struct triangular_ineq*)b;
	if (ta->dist < tb->dist) return -1;
	return 1;
}

/* auxiliary function to sort nodes according to some value from highest to lowest */
int cmpfunc_t_order(const void *a, const void *b)
{
	struct order_elements_by_val *oa = (struct order_elements_by_val*)a;
	struct order_elements_by_val *ob = (struct order_elements_by_val*)b;

	if (oa->val > ob->val) return -1;
	return 1;
}

/*
 * check - check if a solution is a k-partition and compute obj. function
 * (same idea as the one given for lazy constraints, see it)
 * input: matrix x[i][j] = true if and only if i and j belong to the same cluster, for all i < j
 * output: objective value
 */
float check(bool **x)
{
	float objval = 0.0;
	/* ask for memory */
	bool *labeled = new bool[n];
	bool *cluster = new bool[n];
	int *scanned_vertices = new int[n];
	for (int i = 0; i < n; i++) labeled[i] = false;
	bool unlabeled_vertices_exist;
	int vert;
	int count_scanned_vertices = 0;
	do  {
		unlabeled_vertices_exist = false;
		/* search for unlabeled vertices */
		for (int i = 0; i < n; i++) {
			if (labeled[i] == false) {
				vert = i;
				unlabeled_vertices_exist = true;
				scanned_vertices[count_scanned_vertices++] = i;
				labeled[i] = true;
				break;
			}
		}
		if (unlabeled_vertices_exist) {
			/* now, "vert" was unlabeled - scan a clique around "vert" and labels all of its vertices */
			for (int i = 0; i < n; i++) cluster[i] = false;
			cluster[vert] = true;
			for (int j = 0; j < n; j++) {
				if (vert < j) {
					if (x[vert][j] == true) {
						if (labeled[j] == true) bye("Internal error, some triangle inequality is not working!");
						labeled[j] = true;
						cluster[j] = true;
					}
				}
				if (j < vert) {
					if (x[j][vert] == true) {
						if (labeled[j] == true) bye("Internal error, some triangle inequality is not working!");
						labeled[j] = true;
						cluster[j] = true;
					}
				}
			}
			/* show the cluster and computes the value */
			cout << "Cluster " << count_scanned_vertices << ":";
			for (int i = 0; i < n; i++) if (cluster[i] == true) cout << " " << i;
			float akku = 0.0;
			for (int i = 0; i < n - 1; i++) {
				for (int j = i+1; j < n; j++) {
					if (cluster[i] == true && cluster[j] == true) akku += d[i][j];
				}
			}
#ifndef WITH_WEIGHTS
			cout << "  (dist = " << akku << ")" << endl;
#else
			float akku_weight = 0.0;
			for (int i = 0; i < n; i++) if (cluster[i] == true) akku_weight += w[i];
			cout << "  (dist = " << akku << ")  (weight = " << akku_weight << ")";
			if (akku_weight < WL || akku_weight > WU) cout << " <-- OUTSIDE OF [WL, WU] !!!";
			cout << endl;
#endif
			objval += akku;
		}
	} while (unlabeled_vertices_exist);

	if (count_scanned_vertices == k) cout << "Solution is a " << k << "-partition :)" << endl;
	else cout << "Solution is not a " << k << "-partition :(" << endl;

	/* free memory */
	delete[] scanned_vertices;
	delete[] cluster;
	delete[] labeled;
	return objval;
}

/* generates numbers between 0 and N-1 with uniform distribution */
int Random(int N)
{
	return rand() % N;
}

/*
 * tabusearch - use a tabu search heuristic
 * note: "best" must have an initial (not necessarily feasible) solution!
 */
void tabusearch()
{
#if defined(INST_2014) | defined(INST_2015)
	bye("Not implemented!");
#endif
#ifndef WITH_WEIGHTS
	bye("Works only with weights!");
#endif
	int ndivk = n / k;

	/* ask for memory */
	int *current = new int[n];
	int *sizecluster = new int[k]; /* elements of each cluster */
	long double *wcluster = new long double[k]; /* weight of each cluster */
	int **currentv = new int*[k]; /* list of lists with vertices of each cluster */
	for (int j = 0; j < k; j++)	currentv[j] = new int[ndivk + 1];
	int **tenure = new int*[n]; /* life of a (vertex,cluster) in the tabu list */
	for (int i = 0; i < n; i++) {
		tenure[i] = new int[k];
		for (int j = 0; j < k; j++) tenure[i][j] = 0;
	}

	/* set current solution */
	for (int i = 0; i < n; i++) current[i] = best[i];

	/* set initial weights per cluster */
	for (int j = 0; j < k; j++) wcluster[j] = 0.0;
	for (int i = 0; i < n; i++) {
		int j = current[i];
		wcluster[j] += w[i];
	}

	/* compute initial cost ("long double" to diminish numerical errors) */
	long double current_cost = 0.0;
	for (int i = 0; i < n - 1; i++) for (int j = i + 1; j < n; j++) if (current[i] == current[j]) current_cost += d[i][j];
	for (int j = 0; j < k; j++) {
		long double wj = wcluster[j];
		if (wj < WL) current_cost += BIGVALUE;
		if (wj > WU) current_cost += BIGVALUE;
	}
	cout << "Initial solution = " << current_cost << endl;
	best_cost = current_cost;
	set_color(7);

	int maxiter = (int)exp(0.26571*(double)n - 0.052978);
	cout << "Number of iterations estimated = " << maxiter << endl;

	/* make iterations until a limit is reached */
	int iter = 0;
	start_t = ECOclock();
	while (iter < maxiter) {
		/* create solution from current one */
		for (int j = 0; j < k; j++) sizecluster[j] = 0;
		for (int i = 0; i < n; i++) {
			int j = current[i];
			currentv[j][sizecluster[j]++] = i;
		}

		/* show current sol
		cout << "Solution in iteration " << iter << " with cost " << current_cost << ":" << endl;
		for (int j = 0; j < k; j++) {
		int p = 0;
		cout << "  Cluster " << j+1 << ":";
		for (int p = 0; p < sizecluster[j]; p++) cout << " " << currentv[j][p];
		cout << "   (weight = " << wcluster[j] << ")" << endl;
		} */

		/* explore solutions where a vertex of a "bigcluster" is sent to a "lowcluster" */
		long double min_cost = BIGVALUE*BIGVALUE;
		int min_v1 = -1;
		int min_v2 = -1;
		int min_j1 = -1;
		int min_j2 = -1;
		long double min_wj1 = -1.0;
		long double min_wj2 = -1.0;
		int min_movement = -1;
		for (int j1 = 0; j1 < k; j1++) {
			if (sizecluster[j1] == ndivk + 1) {
				for (int p1 = 0; p1 < sizecluster[j1]; p1++) {
					int v1 = currentv[j1][p1];
					for (int j2 = 0; j2 < k; j2++) {
						if (sizecluster[j2] == ndivk) {
							if (tenure[v1][j2] == 0) { /* is (v1,j2) tabu ? */
								/* compute new cost */
								long double new_cost = current_cost;
								for (int pp1 = 0; pp1 < sizecluster[j1]; pp1++) {
									if (pp1 != p1) new_cost -= d[v1][currentv[j1][pp1]];
								}
								for (int pp2 = 0; pp2 < sizecluster[j2]; pp2++) {
									new_cost += d[v1][currentv[j2][pp2]];
								}
								/* compute new weights */
								long double wj1 = wcluster[j1];
								long double wj2 = wcluster[j2];
								long double new_wj1 = wj1 - w[v1];
								long double new_wj2 = wj2 + w[v1];
								/* we also penalize solutions where weights do not hold */
								if (wj1 < WL) new_cost -= BIGVALUE;
								if (wj1 > WU) new_cost -= BIGVALUE;
								if (wj2 < WL) new_cost -= BIGVALUE;
								if (wj2 > WU) new_cost -= BIGVALUE;
								if (new_wj1 < WL) new_cost += BIGVALUE;
								if (new_wj1 > WU) new_cost += BIGVALUE;
								if (new_wj2 < WL) new_cost += BIGVALUE;
								if (new_wj2 > WU) new_cost += BIGVALUE;
								if (new_cost < min_cost) {
									/* if the cost is better, update! */
									min_cost = new_cost;
									min_v1 = v1;
									min_j1 = j1;
									min_j2 = j2;
									min_wj1 = new_wj1;
									min_wj2 = new_wj2;
									min_movement = 1;
								}
							}
						}
					}
				}
			}
		}

		/* explore solutions where two vertices of clusters of the same size are interchanged */
		for (int j1 = 0; j1 < k - 1; j1++) {
			for (int j2 = j1 + 1; j2 < k; j2++) {
				if (sizecluster[j1] == sizecluster[j2]) {
					for (int p1 = 0; p1 < sizecluster[j1]; p1++) {
						int v1 = currentv[j1][p1];
						if (tenure[v1][j2] == 0) { /* is (v1,j2) tabu ? */
							for (int p2 = 0; p2 < sizecluster[j2]; p2++) {
								int v2 = currentv[j2][p2];
								if (tenure[v2][j1] == 0) { /* is (v2,j1) tabu ? */
									/* compute new cost */
									long double new_cost = current_cost;
									for (int pp1 = 0; pp1 < sizecluster[j1]; pp1++) {
										if (pp1 != p1) new_cost += d[v2][currentv[j1][pp1]] - d[v1][currentv[j1][pp1]];
									}
									for (int pp2 = 0; pp2 < sizecluster[j2]; pp2++) {
										if (pp2 != p2) new_cost += d[v1][currentv[j2][pp2]] - d[v2][currentv[j2][pp2]];
									}
									/* compute new weights */
									long double wj1 = wcluster[j1];
									long double wj2 = wcluster[j2];
									long double new_wj1 = wj1 - w[v1] + w[v2];
									long double new_wj2 = wj2 + w[v1] - w[v2];
									/* we also penalize solutions where weights do not hold */
									if (wj1 < WL) new_cost -= BIGVALUE;
									if (wj1 > WU) new_cost -= BIGVALUE;
									if (wj2 < WL) new_cost -= BIGVALUE;
									if (wj2 > WU) new_cost -= BIGVALUE;
									if (new_wj1 < WL) new_cost += BIGVALUE;
									if (new_wj1 > WU) new_cost += BIGVALUE;
									if (new_wj2 < WL) new_cost += BIGVALUE;
									if (new_wj2 > WU) new_cost += BIGVALUE;
									if (new_cost < min_cost) {
										/* if the cost is better, update! */
										min_cost = new_cost;
										min_v1 = v1;
										min_v2 = v2;
										min_j1 = j1;
										min_j2 = j2;
										min_wj1 = new_wj1;
										min_wj2 = new_wj2;
										min_movement = 2;
									}
								}
							}
						}
					}
				}
			}
		}

		/* update current solution */
		if (min_movement == -1) bye("Tenure too high!");
		switch (min_movement) {
		case 1:
			//cout << "Send vertex " << min_v1 << " to cluster " << min_j2 + 1 << ", new cost = " << min_cost << endl;
			current[min_v1] = min_j2;
			current_cost = min_cost;
			wcluster[min_j1] = min_wj1;
			wcluster[min_j2] = min_wj2;
			break;
		case 2:
			//cout << "Exchange vertices " << min_v1 << " and " << min_v2 << " from clusters " << min_j1 + 1 << " and " << min_j2 + 1 << ", new cost = " << min_cost << endl;
			current[min_v1] = min_j2;
			current[min_v2] = min_j1;
			current_cost = min_cost;
			wcluster[min_j1] = min_wj1;
			wcluster[min_j2] = min_wj2;
			break;
		default:
			bye("Tenure too high!");
		}

		/* update best solution if necessary (and if it is better in at least 0.001%) */
		if (min_cost < best_cost * 0.99999) {
			if (best_cost >= BIGVALUE && min_cost < BIGVALUE) cout << "Feasibility reached :)" << endl;
			cout << "New best solution with cost = " << min_cost << " at iteration " << iter << "." << endl;
			best_cost = min_cost;
			for (int i = 0; i < n; i++) best[i] = current[i];
		}

		/* next iteration, mark (v1,j1) as tabu and decrease life of elements in tabu list */
		for (int i = 0; i < n; i++) {
			for (int j = 0; j < k; j++) {
				if (tenure[i][j] > 0) tenure[i][j]--;
			}
		}
		if (min_cost < BIGVALUE) {
			/* when solutions are feasible --> use a low tenure */
			tenure[min_v1][min_j1] = Random(MAXFEASTENURE - MINTENURE + 1) + MINTENURE;
		}
		else {
			/* when solutions are infeasible --> use a high tenure */
			tenure[min_v1][min_j1] = Random(MAXINFEASTENURE - MINTENURE + 1) + MINTENURE;
		}
		iter++;

		/* show time elapsed */
		if (iter % 10000 == 0) {
			cout << " in iteration " << iter << ":  best = " << best_cost << ",  time = " << ECOclock() - start_t << endl;
		}
	}
	cout << "Time elapsed in optimization: " << ECOclock() - start_t << " sec." << endl;

	for (int i = 0; i < n; i++) delete[] tenure[i];
	delete[] tenure;
	for (int j = 0; j < k; j++)	delete[] currentv[j];
	delete[] currentv;
	delete[] sizecluster;
	delete[] current;
}

/*
 * XIJ - return index of variable x(i,j)
 */
int XIJ(int i, int j)
{
	if (i < 0 || i >= n || j < 0 || j >= n || i == j) { cout << "i = " << i << ", j = " << j << endl;  bye("Out of bounds in x(i,j)"); }
	if (i > j) return j + i * n;
	return i + j * n;
}

#ifdef TRIINEQ_LAZYS
/* this block of code is for generating triangular inequalities "on the fly" as lazy constraints (and also other cuts) */
class triineq_lazys : public GRBCallback
{
public:
	GRBVar *vars;
	int lazys_added, lazys_called, cut_called;
	int tricut_added, part2_added;
	int nodecnt;
	float *x;
	order_elements_by_val *order;
	int *Wset, *Tset;
	bool *forbidden;
	bool *Tcharacteristic;

#ifdef WITH_WEIGHTS
	int cover_added, newcuts_added;
#ifdef WEIGHTED_COVER_CUTS_EXHAUSTIVE
	int cover2_added, cover3_added, cover4_added;
#endif
#ifdef WEIGHTED_LOWERBOUND_CUTS_EXHAUSTIVE
	int lower4_added, lower5_added;
#endif
#ifdef WEIGHTED_UPPERBOUND_CUTS_EXHAUSTIVE
	int upper4_added, upper5_added, upper6_added;
#endif
#endif
	triineq_lazys(GRBVar *xvars) {
		vars = xvars;
		lazys_added = 0;
		lazys_called = 0;
		tricut_added = 0;
		part2_added = 0;
		cut_called = 0;
		nodecnt = 0;

		x = new float[n*n];
		order = new order_elements_by_val[n];
		Wset = new int[n];
		Tset = new int[n];
		forbidden = new bool[n];
		Tcharacteristic = new bool[n];

#ifdef WITH_WEIGHTS
		cover_added = 0;
		newcuts_added = 0;
#ifdef WEIGHTED_COVER_CUTS_EXHAUSTIVE
		cover2_added = 0;
		cover3_added = 0;
		cover4_added = 0;
#endif
#ifdef WEIGHTED_LOWERBOUND_CUTS_EXHAUSTIVE
		lower4_added = 0;
		lower5_added = 0;
#endif
#ifdef WEIGHTED_UPPERBOUND_CUTS_EXHAUSTIVE
		upper4_added = 0;
		upper5_added = 0;
		upper6_added = 0;
#endif
#endif
	}

	~triineq_lazys() {
		/* free memory */
		delete[] Tcharacteristic;
		delete[] forbidden;
		delete[] Tset;
		delete[] Wset;
		delete[] order;
		delete[] x;
	}

protected:
	void callback() {
		try {
			switch (where) {
			case GRB_CB_MIP: /* IN MIP (update node count) */
				nodecnt = (int)getDoubleInfo(GRB_CB_MIP_NODCNT);
				if (nodecnt % 100 == 1) {
					cout << "Lazy constraints added: " << lazys_added << endl;
#ifdef ENABLE_CUTS
					cout << "2-partition inequalities added: " << part2_added << endl;
					cout << "triangular cuts added: " << tricut_added << endl;

#ifdef WITH_WEIGHTS
					cout << "Weight-cover inequalities added: " << cover_added << endl;
					cout << "New cuts (weight, lowerbound, upperbound) added: " << newcuts_added << endl;
#endif

#ifdef HISTOGRAM
					cout << "Histogram:" << endl;
					for (int i = 0; i < 10; i++) {
						cout << "  From " << i * 10 << "% to " << (i + 1) * 10 - 1 << "% ---> " << hist[i] << " tri. ineqs." << endl;
					}
#endif

#endif
				}
				break;
			case GRB_CB_MIPSOL: /* LAZY CONSTRAINT ROUTINE */
				lazys_called++;

				/* Found an integer solution - does it satisfy all triangular inequalities? */
				for (int i = 0; i < n*n; i++) x[i] = getSolution(vars[i]);

				for (int ttcount = t_ineq_count_in_formulation; ttcount < t_ineq_count; ttcount++) {
					struct triangular_ineq *ti = &t_ineq[ttcount];
					int i = ti->i;
					int j = ti->j;
					int r = ti->r;
					int sign_ij = ti->sign_ij;
					int sign_ir = ti->sign_ir;
					int sign_jr = ti->sign_jr;
					/* compute l.h.s. of the inequality */
					float akku = 0.0;
					if (sign_ij > 0) akku += x[XIJ(i, j)]; else akku -= x[XIJ(i, j)];
					if (sign_ir > 0) akku += x[XIJ(i, r)]; else akku -= x[XIJ(i, r)];
					if (sign_jr > 0) akku += x[XIJ(j, r)]; else akku -= x[XIJ(j, r)];
					if (akku > 1.0 + EPSILON) {
						/* add the inequality */
						GRBLinExpr lazy = 0.0;
						if (sign_ij > 0) lazy += vars[XIJ(i, j)]; else lazy -= vars[XIJ(i, j)];
						if (sign_ir > 0) lazy += vars[XIJ(i, r)]; else lazy -= vars[XIJ(i, r)];
						if (sign_jr > 0) lazy += vars[XIJ(j, r)]; else lazy -= vars[XIJ(j, r)];
						addLazy(lazy <= 1);
						lazys_added++;
					}
				}
				break;
			case GRB_CB_MIPNODE: /* CUT ROUTINES */
#ifdef ENABLE_CUTS
				if (getIntInfo(GRB_CB_MIPNODE_STATUS) == GRB_OPTIMAL) {
					cut_called++;

					/* Found a non integer solution */
					for (int i = 0; i < n*n; i++) x[i] = getNodeRel(vars[i]);

					//cout << "Fractional variables = " << frac_count << endl;
					//for (int i = 0; i < n - 1; i++) {
					//	for (int j = i + 1; j < n; j++) {
					//		cout << "x(" << i << "," << j << ") = " << x[XIJ(i, j)] << ",   ";
					//	}
					//}
					//cout << endl;

					/* Routine 1: we check for violated triangular inequalities
					 * in root node: check all inequalities
					 * in the remaining nodes: check only the "best" ones */
					int tt_number = t_ineq_count_during_cuts;
					if (nodecnt == 0) tt_number = t_ineq_count;

					bool t_ineq_added = false;
					for (int ttcount = t_ineq_count_in_formulation; ttcount < tt_number; ttcount++) {
						struct triangular_ineq *ti = &t_ineq[ttcount];
						int i = ti->i;
						int j = ti->j;
						int r = ti->r;
						int sign_ij = ti->sign_ij;
						int sign_ir = ti->sign_ir;
						int sign_jr = ti->sign_jr;
						/* compute l.h.s. of the inequality */
						float akku = 0.0;
						if (sign_ij > 0) akku += x[XIJ(i, j)]; else akku -= x[XIJ(i, j)];
						if (sign_ir > 0) akku += x[XIJ(i, r)]; else akku -= x[XIJ(i, r)];
						if (sign_jr > 0) akku += x[XIJ(j, r)]; else akku -= x[XIJ(j, r)];
						if (akku > 1.0 + VIOLATION_TRIINEQ) {
							/* add the inequality */
							GRBLinExpr cut = 0.0;
							if (sign_ij > 0) cut += vars[XIJ(i, j)]; else cut -= vars[XIJ(i, j)];
							if (sign_ir > 0) cut += vars[XIJ(i, r)]; else cut -= vars[XIJ(i, r)];
							if (sign_jr > 0) cut += vars[XIJ(j, r)]; else cut -= vars[XIJ(j, r)];
							addCut(cut <= 1);
							tricut_added++;
							t_ineq_added = true;
#ifdef HISTOGRAM
							int place = (ttcount * 10) / t_ineq_count;
							if (place == 10) place = 9;
							hist[place]++;
#endif
						}
					}

#if defined(WITH_WEIGHTS) && defined(WEIGHTED_COVER_CUTS)
					/* Routine 2: we check for violated cover inequalities
					 * 1) order vertices "i" according to value w_i + sum_{j in V-i} w_j x*_ij in descending order
					 * 2) consider a set with MAXCOVERTEST-2 vertices from the ordering and 2 vertices from V
					 * 3) check for violated 3- 4- cover inequalities by picking elements from every combination of the previous set
					 */

					if (t_ineq_added == false) {
						for (int i = 0; i < n; i++) {
							float val = w[i];
							for (int j = 0; j < n; j++) {
								if (j != i) val += w[j] * x[XIJ(i, j)];
							}
							order[i].v = i;
							order[i].val = val;
						}

						qsort(order, n, sizeof(struct order_elements_by_val), cmpfunc_t_order);

						for (int ii = 0; ii < MAXCOVERTEST - 4; ii++) {
							int i = order[ii].v;
							float wi = w[i];
							for (int jj = ii + 1; jj < MAXCOVERTEST - 3; jj++) {
								int j = order[jj].v;
								/* check necessary condition */
								float wj = w[j];
								if (wi + wj > WU) continue;
								float xij = x[XIJ(i, j)];

								for (int rr = jj + 1; rr < MAXCOVERTEST - 2; rr++) {
									int r = order[rr].v;
									/* check necessary condition */
									float wr = w[r];
									if (wi + wj + wr > WU) {
										/* check 2-cover inequality {i,j,r} (it is seldom violated) */
										//if (wi + wr > WU) continue;
										//if (wj + wr > WU) continue;
										//if (x[XIJ(i, j)] + x[XIJ(i, r)] + x[XIJ(j, r)] <= 1.0 + VIOLATION_COVERINEQ) continue;
										//GRBLinExpr cut = 0.0;
										//cut += vars[XIJ(i, j)] + vars[XIJ(i, r)] + vars[XIJ(j, r)];
										//addCut(cut <= 1); /* no more than 1 edge can happen in {i,j,r} */
										//cover_added++;
										continue;
									}
									float xir = x[XIJ(i, r)];
									float xjr = x[XIJ(j, r)];

									for (int ss = rr + 1; ss < /*MAXCOVERTEST*/ n - 1; ss++) {
										int s = order[ss].v;
										float ws = w[s];
										if (wi + wj + wr + ws > WU) {
											/* check 3-cover inequality {i,j,r,s} */
											if (wi + wj + ws > WU) continue;
											if (wi + wr + ws > WU) continue;
											if (wj + wr + ws > WU) continue;

											float xis = x[XIJ(i, s)];
											float xjs = x[XIJ(j, s)];
											float xrs = x[XIJ(r, s)];
											if (xij + xir + xis + xjr + xjs + xrs <= 3.0 + VIOLATION_COVERINEQ) continue;
											GRBLinExpr cut = 0.0;
											cut += vars[XIJ(i, j)] + vars[XIJ(i, r)] + vars[XIJ(i, s)] + vars[XIJ(j, r)] + vars[XIJ(j, s)] + vars[XIJ(r, s)];
											addCut(cut <= 3); /* no more than a triangle can happen in {i,j,r,s} */
											cover_added++;
											continue;
										}
										float xis = x[XIJ(i, s)];
										float xjs = x[XIJ(j, s)];
										float xrs = x[XIJ(r, s)];

										for (int tt = ss + 1; tt < /* MAXCOVERTEST */ n; tt++) {
											int t = order[tt].v;
											float wt = w[t];
											if (wi + wj + wr + ws + wt > WU) {
												/* check 4-cover inequality {i,j,r,s,t} */
												if (wi + wj + wr + wt > WU) continue;
												if (wi + wj + ws + wt > WU) continue;
												if (wi + wr + ws + wt > WU) continue;
												if (wj + wr + ws + wt > WU) continue;

												float xit = x[XIJ(i, t)];
												float xjt = x[XIJ(j, t)];
												float xrt = x[XIJ(r, t)];
												float xst = x[XIJ(s, t)];
												if (xij + xir + xis + xit + xjr + xjs + xjt + xrs + xrt + xst <= 6.0 + VIOLATION_COVERINEQ) continue;
												GRBLinExpr cut = 0.0;
												cut += vars[XIJ(i, j)] + vars[XIJ(i, r)] + vars[XIJ(i, s)] + vars[XIJ(i, t)];
												cut += vars[XIJ(j, r)] + vars[XIJ(j, s)] + vars[XIJ(j, t)];
												cut += vars[XIJ(r, s)] + vars[XIJ(r, t)] + vars[XIJ(s, t)];
												addCut(cut <= 6); /* no more than a clique of 4 nodes (6 edges) can happen in {i,j,r,s,t} */
												cover_added++;
											}
										}
									}
								}
							}
						}
					}
#endif

#ifdef PART2_CUTS
					/* Routine 3: add 2-partition cuts (Tchai's version separation routine) */
					if (t_ineq_added == false) {
						for (int v = 0; v < n; v++) {
							/* create W set */
							int Wcard = 0;
							for (int u = 0; u < n; u++) {
								if (u != v) {
									float xuv = x[XIJ(u, v)];
									if (xuv > EPSILON && xuv < 1 - EPSILON) Wset[Wcard++] = u;
								}
							}

							if (Wcard >= 5) { /* usually |W| >= 3 is enough, but it is preferable more fractional variables (less depth in the tree) */
								for (int i = 0; i < n; i++) forbidden[i] = false;

								/* repeat this process several times */
								for (int count = 0; count < 5; count++) {
									/* choose two different non-forbidden random vertices from W */
									int i = Wset[Random(Wcard)];
									int j = Wset[Random(Wcard)];

									if (i != j && forbidden[i] == false && forbidden[j] == false) {
										/* create T = {i,j} */
										int Tcard = 0;

										Tset[Tcard++] = i;
										Tset[Tcard++] = j;

										/* find non-forbidden vertices from W - T */
										for (int rindex = 0; rindex < Wcard; rindex++) {
											int r = Wset[rindex];
											if (forbidden[r] == false) {
												bool flag = true;
												for (int tindex = 0; tindex < Tcard; tindex++) {
													if (r == Tset[tindex]) {
														flag = false;
														break;
													}
												}
												if (flag) {
													/* at this point, r is in W - T and not forbidden */
													float akku = x[XIJ(r, v)];
													for (int tindex = 0; tindex < Tcard; tindex++) {
														int t = Tset[tindex];
														akku -= x[XIJ(r, t)];
													}
													/* add r to T when x(r,v) - sum_{t in T} x(r,t) > 0 */
													if (akku > EPSILON) Tset[Tcard++] = r;
												}
											}
										}

										/* try to add the cut */
										float akku = 0;
										for (int tindex = 0; tindex < Tcard; tindex++) {
											int t = Tset[tindex];
											akku += x[XIJ(t, v)];
										}
										for (int tind1 = 0; tind1 < Tcard - 1; tind1++) {
											int t1 = Tset[tind1];
											for (int tind2 = tind1 + 1; tind2 < Tcard; tind2++) {
												int t2 = Tset[tind2];
												akku -= x[XIJ(t1, t2)];
											}
										}
										if (akku > 1.0 + VIOLATION_PART2) {
											GRBLinExpr cut = 0.0;

											for (int tindex = 0; tindex < Tcard; tindex++) {
												int t = Tset[tindex];
												cut += vars[XIJ(t, v)];
											}
											for (int tind1 = 0; tind1 < Tcard - 1; tind1++) {
												int t1 = Tset[tind1];
												for (int tind2 = tind1 + 1; tind2 < Tcard; tind2++) {
													int t2 = Tset[tind2];
													cut -= vars[XIJ(t1, t2)];
												}
											}
											addCut(cut <= 1);
											part2_added++;
											//cout << "2-partition cut added with |T| = " << Tcard << endl;
										}

										/* forbids all vertices of T */
										for (int tindex = 0; tindex < Tcard; tindex++) forbidden[Tset[tindex]] = true;
									}
								}
							}
						}
					}
#endif

#ifdef WITH_WEIGHTS
#ifdef WEIGHTED_COVER_CUTS_EXHAUSTIVE
					/* generate all cover-2 inequalities */
					//for (int i = 0; i < n - 2; i++) {
					//	float wi = w[i];
					//	for (int j = i + 1; j < n - 1; j++) {
					//		float wj = w[j];
					//		if (wi + wj <= WU) {
					//			for (int r = j + 1; r < n; r++) {
					//				float wr = w[r];
					//				if (wi + wr <= WU && wj + wr <= WU && wi + wj + wr > WU) {
					//					/* check if the inequality is violated */
					//					if (x[XIJ(i, j)] + x[XIJ(i, r)] + x[XIJ(j, r)] > 1.0 * (1.0 + VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
					//						GRBLinExpr cut = 0.0;
					//						cut += vars[XIJ(i, j)] + vars[XIJ(i, r)] + vars[XIJ(j, r)];
					//						addCut(cut <= 1); /* no more than 1 edge can happen in {i,j,r} */
					//						newcuts_added++;
					//						cover2_added++;
					//					}
					//				}
					//			}
					//		}
					//	}
					//}
					/* generate all cover-3 inequalities */
					for (int i = 0; i < n - 3; i++) {
						float wi = w[i];
						for (int j = i + 1; j < n - 2; j++) {
							float wj = w[j];
							if (wi + wj <= WU) {
								for (int r = j + 1; r < n - 1; r++) {
									float wr = w[r];
									if (wi + wj + wr <= WU) {
										for (int s = r + 1; s < n; s++) {
											float ws = w[s];
											if (wi + wj + ws <= WU && wi + wr + ws <= WU && wj + wr + ws <= WU && wi + wj + wr + ws > WU) {
												/* check if the inequality is violated */
												if (x[XIJ(i, j)] + x[XIJ(i, r)] + x[XIJ(j, r)] + x[XIJ(i, s)] + x[XIJ(j, s)] + x[XIJ(r, s)] > 3.0 * (1.0 + VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
													GRBLinExpr cut = 0.0;
													cut += vars[XIJ(i, j)] + vars[XIJ(i, r)] + vars[XIJ(j, r)];
													cut += vars[XIJ(i, s)] + vars[XIJ(j, s)] + vars[XIJ(r, s)];
													addCut(cut <= 3); /* no more than a triangle can happen in {i,j,r,s} */
													newcuts_added++;
													cover3_added++;
												}
											}
										}
									}
								}
							}
						}
					}
					/* generate all cover-4 inequalities */
					for (int i = 0; i < n - 4; i++) {
						float wi = w[i];
						for (int j = i + 1; j < n - 3; j++) {
							float wj = w[j];
							if (wi + wj <= WU) {
								for (int r = j + 1; r < n - 2; r++) {
									float wr = w[r];
									if (wi + wj + wr <= WU) {
										for (int s = r + 1; s < n - 1; s++) {
											float ws = w[s];
											if (wi + wj + wr + ws <= WU) {
												for (int t = s + 1; t < n; t++) {
													float wt = w[t];
													if (wi + wj + wr + wt <= WU && wi + wj + ws + wt <= WU && wi + wr + ws + wt <= WU
														&& wj + wr + ws + wt <= WU && wi + wj + wr + ws + wt > WU) {
														/* check if the inequality is violated */
														if (x[XIJ(i, j)] + x[XIJ(i, r)] + x[XIJ(j, r)] + x[XIJ(i, s)] + x[XIJ(j, s)] + x[XIJ(r, s)] + x[XIJ(i, t)] + x[XIJ(j, t)] + x[XIJ(r, t)] + x[XIJ(s, t)] > 6.0 * (1.0 + VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
															GRBLinExpr cut = 0.0;
															cut += vars[XIJ(i, j)] + vars[XIJ(i, r)] + vars[XIJ(j, r)];
															cut += vars[XIJ(i, s)] + vars[XIJ(j, s)] + vars[XIJ(r, s)];
															cut += vars[XIJ(i, t)] + vars[XIJ(j, t)] + vars[XIJ(r, t)] + vars[XIJ(s, t)];
															addCut(cut <= 6); /* no more than a 4-clique (6 edges) can happen in {i,j,r,s,t} */
															newcuts_added++;
															cover4_added++;
														}
													}
												}
											}
										}
									}
								}
							}
						}
					}
#endif

#ifdef WEIGHTED_LOWERBOUND_CUTS_EXHAUSTIVE
					/* explore the Weight-lowerbound inequalities exhaustively for |T| = 4, 5
					(since FL = 5 for n = 44 and k = 8, which is where we make the experiment) */
					if (FL != 5) bye("Error. FL != 5!");

					for (int i = 0; i < n; i++) {
						float wi = w[i];

						/* start with T = {i} */
						for (int j = 0; j < n; j++) Tcharacteristic[j] = false;
						Tcharacteristic[i] = true;
						double wT = wi;

						for (int t2 = 0; t2 < n - 2; t2++) {
							if (t2 == i) continue;
							/* add t2 to T */
							Tcharacteristic[t2] = true;
							wT += w[t2];

							for (int t3 = t2 + 1; t3 < n - 1; t3++) {
								if (t3 == i) continue;
								/* add t3 to T */
								Tcharacteristic[t3] = true;
								wT += w[t3];

								for (int t4 = t3 + 1; t4 < n; t4++) {
									if (t4 == i) continue;
									/* add t4 to T */
									Tcharacteristic[t4] = true;
									wT += w[t4];

									/* check the inequality */
									float r = wT - WL;
									if (r > 0) {
										float rhs = wT - wi;
										float akku = 0.0;
										for (int j = 0; j < n; j++) {
											if (i != j) {
												if (Tcharacteristic[j]) {
													/* variables from T-i */
													float coef = w[j];
													if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
												}
												else {
													/* variables from V-T */
													float coef = w[j] + r;
													if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
												}
											}
										}
										if (akku < rhs * (1.0 - VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
											/* add it as a cut */
											GRBLinExpr cut = 0.0;
											for (int j = 0; j < n; j++) {
												if (i != j) {
													if (Tcharacteristic[j]) {
														/* variables from T-i */
														float coef = w[j];
														if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
													}
													else {
														float coef = w[j] + r;
														if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
													}
												}
											}
											addCut(cut >= rhs);
											newcuts_added++;
											lower4_added++;
										}
									}

									for (int t5 = t4 + 1; t5 < n; t5++) {
										if (t5 == i) continue;
										/* add t5 to T */
										Tcharacteristic[t5] = true;
										wT += w[t5];

										/* check the inequality */
										float r = wT - WL;
										if (r > 0) {
											float rhs = wT - wi;
											float akku = 0.0;
											for (int j = 0; j < n; j++) {
												if (i != j) {
													if (Tcharacteristic[j]) {
														/* variables from T-i */
														float coef = w[j];
														if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
													}
													else {
														/* variables from V-T */
														float coef = w[j] + r;
														if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
													}
												}
											}
											if (akku < rhs * (1.0 - VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
												/* add it as a cut */
												GRBLinExpr cut = 0.0;
												for (int j = 0; j < n; j++) {
													if (i != j) {
														if (Tcharacteristic[j]) {
															/* variables from T-i */
															float coef = w[j];
															if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
														}
														else {
															float coef = w[j] + r;
															if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
														}
													}
												}
												addCut(cut >= rhs);
												newcuts_added++;
												lower5_added++;
											}
										}

										/* remove t5 from T */
										Tcharacteristic[t5] = false;
										wT -= w[t5];
									}
									/* remove t4 from T */
									Tcharacteristic[t4] = false;
									wT -= w[t4];
								}
								/* remove t3 from T */
								Tcharacteristic[t3] = false;
								wT -= w[t3];
							}
							/* remove t2 from T */
							Tcharacteristic[t2] = false;
							wT -= w[t2];
						}
					}
#endif

#ifdef WEIGHTED_UPPERBOUND_CUTS_EXHAUSTIVE
					/* explore the Weight-upperbound inequalities exhaustively for |T| = 4, 5, 6
					(since F_L-1 = 4 and FU = 6 for n = 44 and k = 8, which is where we make the experiment) */
					if (FL != 5 || FU != 6) bye("Error. FL != 5 or FU != 6!");

					for (int i = 0; i < n; i++) {
						float wi = w[i];

						/* start with T = {i} */
						for (int j = 0; j < n; j++) Tcharacteristic[j] = false;
						Tcharacteristic[i] = true;
						double wT = wi;

						for (int t2 = 0; t2 < n - 2; t2++) {
							if (t2 == i) continue;
							/* add t2 to T */
							Tcharacteristic[t2] = true;
							wT += w[t2];

							for (int t3 = t2 + 1; t3 < n - 1; t3++) {
								if (t3 == i) continue;
								/* add t3 to T */
								Tcharacteristic[t3] = true;
								wT += w[t3];

								for (int t4 = t3 + 1; t4 < n; t4++) {
									if (t4 == i) continue;
									/* add t4 to T */
									Tcharacteristic[t4] = true;
									wT += w[t4];

									/* check the inequality */
									float r = WU - wT;
									if (r > 0) {
										float rhs = wT - wi;
										float akku = 0.0;
										for (int j = 0; j < n; j++) {
											if (i != j) {
												if (Tcharacteristic[j]) {
													/* variables from T-i */
													float coef = w[j];
													if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
												}
												else {
													float excess = w[j] - r;
													/* variables from S */
													if (excess > EPSILON) akku += excess * x[XIJ(i, j)];
												}
											}
										}
										if (akku > rhs * (1.0 + VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
											/* add it as a cut */
											GRBLinExpr cut = 0.0;
											for (int j = 0; j < n; j++) {
												if (i != j) {
													if (Tcharacteristic[j]) {
														/* variables from T-i */
														float coef = w[j];
														if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
													}
													else {
														float excess = w[j] - r;
														/* variables from S */
														if (excess > EPSILON) cut += excess * vars[XIJ(i, j)];
													}
												}
											}
											addCut(cut <= rhs);
											newcuts_added++;
											upper4_added++;
										}
									}

									for (int t5 = t4 + 1; t5 < n; t5++) {
										if (t5 == i) continue;
										/* add t5 to T */
										Tcharacteristic[t5] = true;
										wT += w[t5];

										/* check the inequality */
										float r = WU - wT;
										if (r > 0) {
											float rhs = wT - wi;
											float akku = 0.0;
											for (int j = 0; j < n; j++) {
												if (i != j) {
													if (Tcharacteristic[j]) {
														/* variables from T-i */
														float coef = w[j];
														if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
													}
													else {
														float excess = w[j] - r;
														/* variables from S */
														if (excess > EPSILON) akku += excess * x[XIJ(i, j)];
													}
												}
											}
											if (akku > rhs * (1.0 + VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
												/* add it as a cut */
												GRBLinExpr cut = 0.0;
												for (int j = 0; j < n; j++) {
													if (i != j) {
														if (Tcharacteristic[j]) {
															/* variables from T-i */
															float coef = w[j];
															if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
														}
														else {
															float excess = w[j] - r;
															/* variables from S */
															if (excess > EPSILON) cut += excess * vars[XIJ(i, j)];
														}
													}
												}
												addCut(cut <= rhs);
												newcuts_added++;
												upper5_added++;
											}
										}

										for (int t6 = t5 + 1; t6 < n; t6++) {
											if (t6 == i) continue;
											/* add t6 to T */
											Tcharacteristic[t6] = true;
											wT += w[t6];

											/* check the inequality */
											float r = WU - wT;
											if (r > 0) {
												float rhs = wT - wi;
												float akku = 0.0;
												for (int j = 0; j < n; j++) {
													if (i != j) {
														if (Tcharacteristic[j]) {
															/* variables from T-i */
															float coef = w[j];
															if (coef > EPSILON) akku += coef * x[XIJ(i, j)];
														}
														else {
															float excess = w[j] - r;
															/* variables from S */
															if (excess > EPSILON) akku += excess * x[XIJ(i, j)];
														}
													}
												}
												if (akku > rhs * (1.0 + VIOLATION_EXHAUSTIVE_PERCENTAGERHS)) {
													/* add it as a cut */
													GRBLinExpr cut = 0.0;
													for (int j = 0; j < n; j++) {
														if (i != j) {
															if (Tcharacteristic[j]) {
																/* variables from T-i */
																float coef = w[j];
																if (coef > EPSILON) cut += coef * vars[XIJ(i, j)];
															}
															else {
																float excess = w[j] - r;
																/* variables from S */
																if (excess > EPSILON) cut += excess * vars[XIJ(i, j)];
															}
														}
													}
													addCut(cut <= rhs);
													newcuts_added++;
													upper6_added++;
												}
											}

											/* remove t6 from T */
											Tcharacteristic[t6] = false;
											wT -= w[t6];
										}
										/* remove t5 from T */
										Tcharacteristic[t5] = false;
										wT -= w[t5];
									}
									/* remove t4 from T */
									Tcharacteristic[t4] = false;
									wT -= w[t4];
								}
								/* remove t3 from T */
								Tcharacteristic[t3] = false;
								wT -= w[t3];
							}
							/* remove t2 from T */
							Tcharacteristic[t2] = false;
							wT -= w[t2];
						}
					}
#endif
#endif

				}
				break;
#endif
			default: break;
			}
		}
		catch (GRBException e) {
			cout << "Error number: " << e.getErrorCode() << endl;
			cout << e.getMessage() << endl;
		}
		catch (...) {
			cout << "Error during callback" << endl;
		}
	}
};
#endif

/*
 * optimize - create model and optimize it
 */
void optimize()
{
	float objval;

	try {
		GRBEnv env = GRBEnv();
		GRBModel model = GRBModel(env);
		set_color(7);

		/* set some parameters */
		model.getEnv().resetParams();
		model.getEnv().set(GRB_DoubleParam_NodefileStart, 0.5);
		model.getEnv().set(GRB_DoubleParam_TimeLimit, MAXTIME);
		model.getEnv().set(GRB_DoubleParam_MIPGap, 0.0);
		model.getEnv().set(GRB_DoubleParam_MIPGapAbs, 0.0);
		model.getEnv().set(GRB_DoubleParam_IntFeasTol, EPSILON);
		model.getEnv().set(GRB_IntParam_Threads, 1);
		model.getEnv().set(GRB_IntParam_Seed, 1);

#ifdef PUREBYB
		model.getEnv().set(GRB_IntParam_Presolve, GRB_PRESOLVE_CONSERVATIVE);
		model.getEnv().set(GRB_IntParam_Cuts, GRB_CUTS_OFF);
		model.getEnv().set(GRB_IntParam_RINS, 0);
		model.getEnv().set(GRB_DoubleParam_Heuristics, 0.0);
		model.getEnv().set(GRB_IntParam_PumpPasses, 0);
#endif

		/* create binary variables x(i,j)
		 * this is somewhat ugly: those variables such that i >= j are also created (for the sake of simplicity)
		 * but are set to zero, and are never further accessed
		 * also, variables whose distance is at least BIGVALUE
		 * also, variables such that w_i + w_j > W_U (weighted version)
		*/
		GRBVar *vars = new GRBVar[n*n];
		char str[128];
		int index = 0;
		int nonzero = 0;
#ifdef WITH_WEIGHTS
		int c1_count = 0;
#endif

		for (int j = 0; j < n; j++) {
			for (int i = 0; i < n; i++) {
				sprintf(str, "x(%d,%d)", i, j);
				float ub = 1.0;
				if (i >= j) ub = 0.0;
				if (d[i][j] >= BIGVALUE - EPSILON) ub = 0.0;
#ifdef WITH_WEIGHTS
				if (i < j && w[i] + w[j] > WU) {
					ub = 0.0;
					c1_count++;
				}
#endif
				vars[index++] = model.addVar(0.0, ub, 0.0, GRB_BINARY, str);
				if (ub > EPSILON) nonzero++;
			}
		}
		set_color(11);
		cout << "Number of total variables added: " << index << " (non-zero variables: " << nonzero << ")" << endl;
#ifdef WITH_WEIGHTS
		cout << "Number of 1-cover equalities (x_ij = 0): " << c1_count << endl;
#endif
		model.update();

		/* generate objective function */
		GRBLinExpr fobj = 0.0;
		for (int i = 0; i < n - 1; i++) for (int j = i + 1; j < n; j++) fobj += d[i][j] * vars[XIJ(i, j)];
		model.setObjective(fobj, GRB_MINIMIZE);

#ifdef WITH_WEIGHTS
		/* generate weights constraints */
		int w_count = 0;
		for (int i = 0; i < n; i++) {
			GRBLinExpr restr = 0.0;
			for (int j = 0; j < n; j++) if (i != j) restr += w[j] * vars[XIJ(i, j)];
			model.addConstr(restr >= WL - w[i]); /* sum is between WL-wi and WU-wi */
			model.addConstr(restr <= WU - w[i]);
			w_count += 2;
		}
		cout << "Number of weight constraints added: " << w_count << endl;
#endif

#if defined(INST_2014) || defined(INST_2015)
		/* generate "amazonas <= 3" and "cabezas <= 1" constraint:
		 * For "amazonas <= 3" we do the following:
		 *  Since at most 3 cities of {13,14,15,16,19} can belong to the same cluster
		 *  we impose that LESS than 6 edges (a complete of 4 vertices) can be chosen from the
		 *  edges in those cities (which are 10 edges in total).
		 *  Since it is impossible to choose 5 edges from the 10 available, we use "<= 4"
		 * For "cabezas <= 1" we simply set to zero the edges {0,1}, {0,2}, {1,2} */
		const int amazonas[] = { 13, 14, 15, 16, 19 };
		GRBLinExpr amrestr = 0.0;
		for (int ai = 0; ai < 4; ai++) {
			for (int aj = ai + 1; aj < 5; aj++) {
				int i = amazonas[ai];
				int j = amazonas[aj];
				amrestr += vars[XIJ(i, j)];
			}
		}
		model.addConstr(amrestr <= 4);
		model.addConstr(vars[XIJ(0, 1)] == 0);
		model.addConstr(vars[XIJ(1, 2)] == 0);
		model.addConstr(vars[XIJ(0, 2)] == 0);
		cout << "Number of amazonas constraints added: 4" << endl;
#endif

		/* generate equitable constraints */
		int eq_count = 0;
		int r = n % k;
		int alpha = n / k;
		for (int i = 0; i < n; i++) {
			GRBLinExpr restr = 0.0;
			for (int j = 0; j < n; j++) if (i != j) restr += vars[XIJ(i, j)];
			if (r == 0) {
				/* sum is equal to n/k - 1 */
				model.addConstr(restr == alpha - 1);
				FL = alpha;
				FU = alpha;
				eq_count++;
			}
			else {
				model.addConstr(restr >= alpha - 1); /* sum is between floor(n/k)-1 and ceil(n/k)-1 */
				model.addConstr(restr <= alpha);
				FL = alpha;
				FU = alpha + 1;
				eq_count += 2;
			}
		}
		cout << "Number of equitable constraints added: " << eq_count << endl;

		/* generate triangular inequalities:
		 * we sort by distances and we put first those inequalities whose distances are small */
		t_ineq = new struct triangular_ineq[n * (n - 1) * (n - 2) / 2];
		int ttcount = 0;
		for (int i = 0; i < n - 2; i++) {
			for (int j = i + 1; j < n - 1; j++) {
				for (int r = j + 1; r < n; r++) {
					/* first combination */
					if (d[i][j] < BIGVALUE && d[i][r] < BIGVALUE) {
						t_ineq[ttcount].i = i; t_ineq[ttcount].sign_ij = +1;
						t_ineq[ttcount].j = j; t_ineq[ttcount].sign_ir = +1;
						t_ineq[ttcount].r = r; t_ineq[ttcount].sign_jr = -1;
						t_ineq[ttcount].dist = d[i][j] + d[i][r];
						ttcount++;
					}
					/* second combination */
					if (d[i][j] < BIGVALUE && d[j][r] < BIGVALUE) {
						t_ineq[ttcount].i = i; t_ineq[ttcount].sign_ij = +1;
						t_ineq[ttcount].j = j; t_ineq[ttcount].sign_ir = -1;
						t_ineq[ttcount].r = r; t_ineq[ttcount].sign_jr = +1;
						t_ineq[ttcount].dist = d[i][j] + d[j][r];
						ttcount++;
					}
					/* third combination */
					if (d[i][r] < BIGVALUE && d[j][r] < BIGVALUE) {
						t_ineq[ttcount].i = i; t_ineq[ttcount].sign_ij = -1;
						t_ineq[ttcount].j = j; t_ineq[ttcount].sign_ir = +1;
						t_ineq[ttcount].r = r; t_ineq[ttcount].sign_jr = +1;
						t_ineq[ttcount].dist = d[i][r] + d[j][r];
						ttcount++;
					}
				}
			}
		}
		t_ineq_count = ttcount;

#ifdef TRIINEQ_LAZYS
		qsort(t_ineq, t_ineq_count, sizeof(struct triangular_ineq), cmpfunc_t_ineq);
#endif

		t_ineq_count_in_formulation = (t_ineq_count * PERCENTAGE_TRIINEQ_IN_FORM) / 100;
		t_ineq_count_during_cuts = (t_ineq_count * PERCENTAGE_TRIINEQ_IN_CUTS) / 100;
		for (ttcount = 0; ttcount < t_ineq_count_in_formulation; ttcount++) {
			int i = t_ineq[ttcount].i;
			int j = t_ineq[ttcount].j;
			int r = t_ineq[ttcount].r;
			GRBLinExpr restr = 0.0;
			if (t_ineq[ttcount].sign_ij > 0) restr += vars[XIJ(i, j)]; else restr -= vars[XIJ(i, j)];
			if (t_ineq[ttcount].sign_ir > 0) restr += vars[XIJ(i, r)]; else restr -= vars[XIJ(i, r)];
			if (t_ineq[ttcount].sign_jr > 0) restr += vars[XIJ(j, r)]; else restr -= vars[XIJ(j, r)];
			model.addConstr(restr <= 1);
		}
		cout << "Number of triangular constraints added: " << t_ineq_count_in_formulation << "  (" << PERCENTAGE_TRIINEQ_IN_FORM << "% of " << t_ineq_count << ")" << endl;

#ifdef BETA_EQUALITY
		/* generate equality: sum of x_ij = r e(alpha+1) + (k - r) e(alpha), only if r != 0
		 * they are needed ONLY if n does not divide k; they also replace clique inequalities */
		if (r > 0) {
			int beta = alpha * (r * (alpha + 1) + (k - r) * (alpha - 1)) / 2;
			cout << "Equality: Beta = " << beta << endl;
			GRBLinExpr restr = 0.0;
			for (int i = 0; i < n - 1; i++) for (int j = i + 1; j < n; j++) restr += vars[XIJ(i, j)];
			model.addConstr(restr == beta);
		}
#else
		cout << "Beta equality not added. Caution: partition of size k may not be guaranteed!" << endl;
#endif

		start_t = ECOclock();
#ifdef USE_TABU_SEARCH
		/* execute tabu search with a trivial initial solution */
		best = new int[n];
		int cc = 0;
		for (int i = 0; i < n; i++) {
			best[i] = cc++;
			if (cc == k) cc = 0;
		}
		tabusearch();

		/* inject solution */
		bool **xtabu = new bool*[n];
		for (int i = 0; i < n; i++) {
			xtabu[i] = new bool[n];
			for (int j = 0; j < n; j++) {
				xtabu[i][j] = false;
				if (i != j)	if (best[i] == best[j]) xtabu[i][j] = true;
			}
		}
		objval = check(xtabu);
		cout << "Inject solution with objval = " << objval << endl;
		for (int i = 0; i < n - 1; i++) {
			for (int j = i + 1; j < n; j++) vars[XIJ(i, j)].set(GRB_DoubleAttr_Start, xtabu[i][j] ? 1 : 0);
		}
		for (int i = 0; i < n; i++) delete[] xtabu[i];
		delete[] xtabu;
		delete[] best;

		/* turn off GuRoBi heuristics and focus on proving optimality */
		model.getEnv().set(GRB_DoubleParam_Heuristics, 0.0);
		model.getEnv().set(GRB_IntParam_PumpPasses, 0);
		model.getEnv().set(GRB_IntParam_RINS, 0);
#endif

		model.update();
		//model.write("kpartmodel.lp");
		//cout << "MIP saved to kpartmodel.lp." << endl;


#ifdef TRIINEQ_LAZYS
		cout << "Triangular ineq. are added as lazy constraints." << endl;
#ifdef ENABLE_CUTS
		cout << "Cuts enabled." << endl;
		model.getEnv().set(GRB_IntParam_PreCrush, 1);
		/* disable some GuRoBi cuts that are seldom used */
		model.getEnv().set(GRB_IntParam_CliqueCuts, 0);
		model.getEnv().set(GRB_IntParam_FlowCoverCuts, 0);
		model.getEnv().set(GRB_IntParam_FlowPathCuts, 0);
		model.getEnv().set(GRB_IntParam_ImpliedCuts, 0);
		model.getEnv().set(GRB_IntParam_MIPSepCuts, 0);
		model.getEnv().set(GRB_IntParam_NetworkCuts, 0);
		model.getEnv().set(GRB_IntParam_ModKCuts, 0);

#ifdef HISTOGRAM
		for (int i = 0; i < 10; i++) hist[i] = 0;
#endif
#endif
		model.getEnv().set(GRB_IntParam_LazyConstraints, 1);
		triineq_lazys cb = triineq_lazys(vars);
		model.setCallback(&cb);
#endif

		set_color(7);
		cout << "Starting optimization..." << endl;
		model.optimize();
		cout << "Time elapsed in optimization: " << ECOclock() - start_t << " sec." << endl;

#ifdef TRIINEQ_LAZYS
		cout << "Lazy constraints added: " << cb.lazys_added << "  (lazy callbacks = " << cb.lazys_called << ")." << endl;
#ifdef ENABLE_CUTS
		cout << "Added cuts --> 2-part: " << cb.part2_added << ", triang.: " << cb.tricut_added << "  (cut callbacks = " << cb.cut_called << ")." << endl;
#ifdef WITH_WEIGHTS		
		cout << "(Weight) cuts --> cover: " << cb.cover_added << ", new cuts: " << cb.newcuts_added << endl;
#ifdef WEIGHTED_COVER_CUTS_EXHAUSTIVE
		cout << "Weight-Cover as \"new cuts\":  2 --> " << cb.cover2_added << ", 3 --> " << cb.cover3_added << ", 4 --> " << cb.cover4_added << endl;
#endif
#ifdef WEIGHTED_LOWERBOUND_CUTS_EXHAUSTIVE
		cout << "Weight-Lowerbound as \"new cuts\":  4 --> " << cb.lower4_added << ", 5 --> " << cb.lower5_added << endl;
#endif
#ifdef WEIGHTED_UPPERBOUND_CUTS_EXHAUSTIVE
		cout << "Weight-Upperbound as \"new cuts\":  4 --> " << cb.upper4_added << ", 5 --> " << cb.upper5_added << ", 6 --> " << cb.upper6_added << endl;
#endif
#endif
#endif
#endif

		int status = model.get(GRB_IntAttr_Status);
		if (status != GRB_OPTIMAL) {
			switch (status) {
			case GRB_INF_OR_UNBD:
			case GRB_INFEASIBLE: bye("Infeasible :(");
			case GRB_TIME_LIMIT: bye("Time limit reached :(");
			default: bye("Unexpected error :(");
			}
		}

		set_color(10);
		cout << "Optimal solution found :)" << endl;
		set_color(7);

		/* check solution */
		bool **xsol = new bool*[n];
		for (int i = 0; i < n; i++) {
			xsol[i] = new bool[n];
			for (int j = 0; j < n; j++) xsol[i][j] = false;
		}
		for (int i = 0; i < n - 1; i++) {
			for (int j = i + 1; j < n; j++) if (vars[XIJ(i, j)].get(GRB_DoubleAttr_X) > 0.2) xsol[i][j] = true;
		}
		float objval = check(xsol);
		cout << "Optimal value = " << objval << "   (opt*4 = " << objval * 4 << ")" << endl;

		for (int i = 0; i < n; i++) delete[] xsol[i];
		delete[] xsol;

		return;
	}
	catch (GRBException e) {
		cout << "Error code = " << e.getErrorCode() << endl;
		cout << e.getMessage() << endl;
	}
	catch (...) {
		cout << "Exception during optimization" << endl;
	}
	bye("Chau!");
}

/*
 * main - Funcion principal
 */
int main(int argc, char **argv)
{
	set_color(15);
	cout << "KPART - Solver of weighted balanced k-clique partitioning problem." << endl;
	cout << "Made by Daniel Severin." << endl;
	cout << "Partially supported by 15MathAmSud06 project and CONICET." << endl;

	/* check if there are enough arguments */
	if (argc < 2) {
		cout << "Usage: kpart file.txt" << endl;
		cout << "where: file.txt = file with the instance" << endl;
		exit(1);
	}

	/* read number of points and partitions */
	char *filename = argv[1];
	FILE *stream = fopen(filename, "rt");
	if (!stream) {
		perror("Error openinig file");
		bye("Bye!");
	}
	fscanf(stream, "%d", &n); /* first line: number of vertices */
	fscanf(stream, "%d", &k); /* second line : number of clusters */
	if (n <= 3 || n >= 10000 || k <= 1 || k >= 5000) bye("Out of bounds!");

	/* ask for memory */
	w = new float[n];
	d = new float*[n];
	for (int i = 0; i < n; i++) d[i] = new float[n];

	/* read weights and distance matrix */
	float max_dist = -BIGVALUE;
	for (int i = 0; i < n; i++) {
		fscanf(stream, "%f", &w[i]);
		for (int j = 0; j < n; j++) {
			float dist;
			fscanf(stream, "%f", &dist);
			d[i][j] = dist;
			if (dist > max_dist) max_dist = dist;
		}
	}
	fclose(stream);

	set_color(14);
	cout << "We have n = " << n << " and k = " << k << "." << endl;
	cout << "Maximum distance found = " << max_dist << endl;

#ifdef WITH_WEIGHTS
	/* compute WL and WU */
	float mu = 0.0;
	for (int i = 0; i < n; i++) mu += w[i];
	mu /= (float)n;
	float sigma = 0.0;
	for (int i = 0; i < n; i++) sigma += (w[i] - mu) * (w[i] - mu);
	sigma = sqrt(sigma / (float)(n-1));
	//Lo que yo creo:
	//WL = ((float)n / (float)k) * (mu - sigma);
	//WU = ((float)n / (float)k) * (mu + sigma);
	WL = ((float)n / (float)k) * mu - sigma;
	WU = ((float)n / (float)k) * mu + sigma;

#ifdef INST_2014
	WL = 2.50717; WU = 4.82616;
	sigma = 1.159495;
#endif

#ifdef INST_2015
	WL = 2.37586; WU = 4.64914;
	sigma = 1.13664;
#endif

#ifdef INST_2015_PROPOSAL
	WL = 2.08412; WU = 4.14835;
	sigma = 1.032115;

	cout << "Forbidding pairs:";
	for (int i = 0; i < n; i += 2) {
		d[i][i + 1] = BIGVALUE;
		d[i + 1][i] = BIGVALUE;
		cout << " (" << i << "," << i + 1 << ")";
	}
	cout << endl;
#endif

	cout << "Weights: mu = " << mu << ", sigma = " << sigma << ", WL = " << WL << ", WU = " << WU << "." << endl;
	/* now we add artificial nodes (only for "weighted" version) */
#ifdef ADD_ARTIFICIAL_NODES
	int rest = n % k;
	if (rest != 0) {
		/* define a new instance with additional vertices */
		int nn = n + (k - rest);
		cout << "Creating new instance with " << nn << " vertices." << endl;
		float *ww = new float[nn];
		float **dd = new float*[nn];
		for (int i = 0; i < nn; i++) dd[i] = new float[nn];

		/* copy data from the old instance */
		for (int i = 0; i < nn; i++) {
			if (i < n) ww[i] = w[i];
			else ww[i] = 0.0;
			for (int j = 0; j < nn; j++) {
				if (i < n && j < n) dd[i][j] = d[i][j];
				else {
					if (i >= n && j >= n) dd[i][j] = BIGVALUE; /* <-- two artificial nodes can not be in the same cluster */
					else dd[i][j] = 0.0;
				}
			}
		}

		/* swap instances */
		for (int i = 0; i < n; i++) delete[] d[i];
		delete[] d;
		delete[] w;
		n = nn;
		w = ww;
		d = dd;
	}
#endif
#endif

	/* compute cardinals of S and T for 2-partition ineq. to be facet defining (Labbe et al.) */
	//int r = n % k;
	//if (r == 0) bye("k must not divide n");
	//int alpha = n / k;
	//cout << "2-part. ineq: r = " << r << ", F_L = " << alpha << ", F_U = " << alpha + 1 << endl;
	//int kval = n / alpha; /* Case 1: kval = floor(n/F_L)*/
	//int tcase = 1;
	//if (n % alpha != 1) {
	//	kval = 1 + n / (alpha + 1); /* Case 2: kval = ceil(n/F_U) */
	//	tcase = 2;
	//	if (n % (alpha + 1) != alpha) bye("No case!");
	//}
	//cout << "   kval = " << kval << ", case = " << tcase << endl;
	//for (int s = 1; s <= n - 2; s++) {
	//	for (int t = s + 1; t <= n - 1; t++) {
	//		int delta = t - s;
	//		if (delta < kval) {
	//			if ((tcase == 1 && s <= delta * ((alpha - 1) / 2) + (kval - delta - 1) * (alpha / 2) + (alpha - 2) / 2) ||
	//				(tcase == 2 && s <= (alpha - 1) / 2 + (kval - delta - 1) * (alpha + 1) / 2 + (delta - 1) * (alpha / 2) + (alpha - 3) / 2)) {
	//				cout << "   |S| = " << s << ",  |T| = " << t << endl;
	//			}
	//		}
	//	}
	//}

	/* optimize it */
	optimize();
	set_color(7);

	/* free memory */
	for (int i = 0; i < n; i++) delete[] d[i];
	delete[] d;
	delete[] w;
	return 0;
}