/*
 * qap.c
 *
 *  Created on: 12/12/2016
 *      Author: pedro
 *
 *		http://anjos.mgi.polymtl.ca/qaplib//
 *
 *		Quadratic assignment problem:
 *      There are a set of n facilities and a set of n locations. For each pair of locations,
 *      a distance is specified and for each pair of facilities a weight or flow is specified
 *      (e.g., the amount of supplies transported between the two facilities). The problem is
 *      to assign all facilities to different locations with the goal of minimizing the sum of
 *      the distances multiplied by the corresponding flows.
 */

#include "qap.h"

#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "../config.h"
#include "../constraints/all_different.h"
#include "../constraints/element_var.h"
#include "../constraints/linear_var.h"
#include "../constraints/minimize.h"
#include "../split.h"
#include "../variables.h"

unsigned int p2_flow[] = { 0, 1, 3, 0 };
unsigned int p2_dist[] = { 0, 4, 2, 0 };

unsigned int p2c_flow[] = { 0, 1, 0, 0 };
unsigned int p2c_dist[] = { 3, 5, 7, 9 };

unsigned int p3_flow[] = { 0, 1, 4, 1, 0, 2, 4, 2, 0 };
unsigned int p3_dist[] = { 0, 2, 4, 2, 0, 1, 4, 1, 0 };

unsigned int p3b_flow[] = { 0, 1, 1, 1, 0, 1, 1, 1, 0 };
unsigned int p3b_dist[] = { 0, 1, 1, 1, 0, 1, 1, 1, 0 };

unsigned int p3c_flow[] = { 0, 0, 1, 0, 0, 0, 0, 0, 0 };
unsigned int p3c_dist[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 };

unsigned int p5_flow[] = { 0, 1, 2, 1, 0, 1, 0, 2, 1, 0, 2, 2, 0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0 };
unsigned int p5_dist[] = { 0, 5, 1, 3, 2, 5, 0, 2, 1, 4, 1, 2, 0, 2, 1, 3, 1, 2, 0, 2, 2, 4, 1, 2, 0 };

unsigned int p10_flow[] = { 0, 1, 2, 1, 0, 0, 1, 2, 1, 0, 1, 0, 2, 1, 0, 0, 1, 2, 0, 1, 2, 2, 0, 1, 0, 0, 1, 0, 2, 2, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0,
		0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 2, 2, 0, 1, 0, 0, 1, 0, 2, 2, 1, 0, 2, 1, 0, 0, 1, 2, 0, 1, 0, 1, 2, 1, 0, 0, 1, 2, 1, 0 };
unsigned int p10_dist[] = { 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0,
		1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0 };

unsigned int choi15_flow[] = {
  16,  4,  6,  8,  7,  4,  6,  8,  7,  4,  6,  8,  7,  4,  6,
   4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,
   6,  1,  4,  2,  6,  1,  9,  2,  6,  1,  9,  2,  6,  1,  9,
   8,  3,  2, 10,  8,  3,  2,  5,  8,  3,  2,  5,  8,  3,  2,
   7,  4,  6,  8, 16,  4,  6,  8,  7,  4,  6,  8,  7,  4,  6,
   4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,
   6,  1,  9,  2,  6,  1,  4,  2,  6,  1,  9,  2,  6,  1,  9,
   8,  3,  2,  5,  8,  3,  2, 10,  8,  3,  2,  5,  8,  3,  2,
   7,  4,  6,  8,  7,  4,  6,  8, 16,  4,  6,  8,  7,  4,  6,
   4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,
   6,  1,  9,  2,  6,  1,  9,  2,  6,  1,  4,  2,  6,  1,  9,
   8,  3,  2,  5,  8,  3,  2,  5,  8,  3,  2, 10,  8,  3,  2,
   7,  4,  6,  8,  7,  4,  6,  8,  7,  4,  6,  8, 16,  4,  6,
   4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,  3,  4,  6,  1,
   6,  1,  9,  2,  6,  1,  9,  2,  6,  1,  9,  2,  6,  1,  4
};

unsigned int choi15_dist[] = {
  1, 1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2,
  1, 2, 1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3,
  2, 1, 3, 1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4,
  3, 2, 1, 4, 1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5,
  4, 3, 2, 1, 5, 1, 2, 3, 4, 5, 6, 7, 8, 7, 6,
  5, 4, 3, 2, 1, 6, 1, 2, 3, 4, 5, 6, 7, 8, 7,
  6, 5, 4, 3, 2, 1, 7, 1, 2, 3, 4, 5, 6, 7, 8,
  7, 6, 5, 4, 3, 2, 1, 8, 1, 2, 3, 4, 5, 6, 7,
  8, 7, 6, 5, 4, 3, 2, 1, 7, 1, 2, 3, 4, 5, 6,
  7, 8, 7, 6, 5, 4, 3, 2, 1, 6, 1, 2, 3, 4, 5,
  6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 5, 1, 2, 3, 4,
  5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 4, 1, 2, 3,
  4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 3, 1, 2,
  3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 2, 1,
  2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 1
};

unsigned int esc16i_flow[] = {
  0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  2, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  1, 2, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 1, 2, 0, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 1, 3, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 1, 2, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 1, 2, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 1, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};

unsigned int esc16i_dist[] = {
  0, 0, 0, 1, 0, 1, 1, 2, 0, 1, 1, 2, 1, 2, 2, 3,
  0, 0, 1, 0, 1, 0, 2, 1, 1, 0, 2, 1, 2, 1, 3, 2,
  0, 1, 0, 0, 1, 2, 0, 1, 1, 2, 0, 1, 2, 3, 1, 2,
  1, 0, 0, 0, 2, 1, 1, 0, 2, 1, 1, 0, 3, 2, 2, 1,
  0, 1, 1, 2, 0, 0, 0, 1, 1, 2, 2, 3, 0, 1, 1, 2,
  1, 0, 2, 1, 0, 0, 1, 0, 2, 1, 3, 2, 1, 0, 2, 1,
  1, 2, 0, 1, 0, 1, 0, 0, 2, 3, 1, 2, 1, 2, 0, 1,
  2, 1, 1, 0, 1, 0, 0, 0, 3, 2, 2, 1, 2, 1, 1, 0,
  0, 1, 1, 2, 1, 2, 2, 3, 0, 0, 0, 1, 0, 1, 1, 2,
  1, 0, 2, 1, 2, 1, 3, 2, 0, 0, 1, 0, 1, 0, 2, 1,
  1, 2, 0, 1, 2, 3, 1, 2, 0, 1, 0, 0, 1, 2, 0, 1,
  2, 1, 1, 0, 3, 2, 2, 1, 1, 0, 0, 0, 2, 1, 1, 0,
  1, 2, 2, 3, 0, 1, 1, 2, 0, 1, 1, 2, 0, 0, 0, 1,
  2, 1, 3, 2, 1, 0, 2, 1, 1, 0, 2, 1, 0, 0, 1, 0,
  2, 3, 1, 2, 1, 2, 0, 1, 1, 2, 0, 1, 0, 1, 0, 0,
  3, 2, 2, 1, 2, 1, 1, 0, 2, 1, 1, 0, 1, 0, 0, 0
};

unsigned int esc16e_flow[] = {
  0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 1, 1, 0, 2, 3, 2, 0, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 1, 1, 2, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 1, 1, 3, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0,
  1, 0, 0, 1, 1, 2, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};

unsigned int esc16e_dist[] = {
  0, 0, 0, 1, 0, 1, 1, 2, 0, 1, 1, 2, 1, 2, 2, 3,
  0, 0, 1, 0, 1, 0, 2, 1, 1, 0, 2, 1, 2, 1, 3, 2,
  0, 1, 0, 0, 1, 2, 0, 1, 1, 2, 0, 1, 2, 3, 1, 2,
  1, 0, 0, 0, 2, 1, 1, 0, 2, 1, 1, 0, 3, 2, 2, 1,
  0, 1, 1, 2, 0, 0, 0, 1, 1, 2, 2, 3, 0, 1, 1, 2,
  1, 0, 2, 1, 0, 0, 1, 0, 2, 1, 3, 2, 1, 0, 2, 1,
  1, 2, 0, 1, 0, 1, 0, 0, 2, 3, 1, 2, 1, 2, 0, 1,
  2, 1, 1, 0, 1, 0, 0, 0, 3, 2, 2, 1, 2, 1, 1, 0,
  0, 1, 1, 2, 1, 2, 2, 3, 0, 0, 0, 1, 0, 1, 1, 2,
  1, 0, 2, 1, 2, 1, 3, 2, 0, 0, 1, 0, 1, 0, 2, 1,
  1, 2, 0, 1, 2, 3, 1, 2, 0, 1, 0, 0, 1, 2, 0, 1,
  2, 1, 1, 0, 3, 2, 2, 1, 1, 0, 0, 0, 2, 1, 1, 0,
  1, 2, 2, 3, 0, 1, 1, 2, 0, 1, 1, 2, 0, 0, 0, 1,
  2, 1, 3, 2, 1, 0, 2, 1, 1, 0, 2, 1, 0, 0, 1, 0,
  2, 3, 1, 2, 1, 2, 0, 1, 1, 2, 0, 1, 0, 1, 0, 0,
  3, 2, 2, 1, 2, 1, 1, 0, 2, 1, 1, 0, 1, 0, 0, 0
};

unsigned int esc16g_flow[] = {
		0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 1, 1, 2, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 1, 0, 1, 1, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 1, 1, 0, 1, 2, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0,
		1, 2, 1, 1, 0, 3, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
		1, 1, 0, 2, 3, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 1, 2, 3, 1, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
		0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};

unsigned int esc16g_dist[] = {
		0, 0, 0, 1, 0, 1, 1, 2, 0, 1, 1, 2, 1, 2, 2, 3,
		0, 0, 1, 0, 1, 0, 2, 1, 1, 0, 2, 1, 2, 1, 3, 2,
		0, 1, 0, 0, 1, 2, 0, 1, 1, 2, 0, 1, 2, 3, 1, 2,
		1, 0, 0, 0, 2, 1, 1, 0, 2, 1, 1, 0, 3, 2, 2, 1,
		0, 1, 1, 2, 0, 0, 0, 1, 1, 2, 2, 3, 0, 1, 1, 2,
		1, 0, 2, 1, 0, 0, 1, 0, 2, 1, 3, 2, 1, 0, 2, 1,
		1, 2, 0, 1, 0, 1, 0, 0, 2, 3, 1, 2, 1, 2, 0, 1,
		2, 1, 1, 0, 1, 0, 0, 0, 3, 2, 2, 1, 2, 1, 1, 0,
		0, 1, 1, 2, 1, 2, 2, 3, 0, 0, 0, 1, 0, 1, 1, 2,
		1, 0, 2, 1, 2, 1, 3, 2, 0, 0, 1, 0, 1, 0, 2, 1,
		1, 2, 0, 1, 2, 3, 1, 2, 0, 1, 0, 0, 1, 2, 0, 1,
		2, 1, 1, 0, 3, 2, 2, 1, 1, 0, 0, 0, 2, 1, 1, 0,
		1, 2, 2, 3, 0, 1, 1, 2, 0, 1, 1, 2, 0, 0, 0, 1,
		2, 1, 3, 2, 1, 0, 2, 1, 1, 0, 2, 1, 0, 0, 1, 0,
		2, 3, 1, 2, 1, 2, 0, 1, 1, 2, 0, 1, 0, 1, 0, 0,
		3, 2, 2, 1, 2, 1, 1, 0, 2, 1, 1, 0, 1, 0, 0, 0
};

/*
 * Solve one of the above instances of the Quadratic assignment problem
 */
void run_qap(char* csp_inst_name) {
	int n;
	int* p_flow;
	unsigned int* p_dist;
	unsigned long result;
	unsigned int d_max;
	unsigned int cost;
	unsigned int* n_pos;
	unsigned int* new_dist;
	unsigned int vs[3];
	bool symmetrical = true;
	bool symmetrical_aux;
	int i, j;

	if (!strcmp(csp_inst_name, "p2")) {
		n = 2;
		p_flow = (int*) p2_flow;
		p_dist = p2_dist;
	} else if (!strcmp(csp_inst_name, "p2c")) {
		n = 2;
		p_flow = (int*) p2c_flow;
		p_dist = p2c_dist;
	} else if (!strcmp(csp_inst_name, "p3")) {
		n = 3;
		p_flow = (int*) p3_flow;
		p_dist = p3_dist;
	} else if (!strcmp(csp_inst_name, "p3b")) {
		n = 3;
		p_flow = (int*) p3b_flow;
		p_dist = p3b_dist;
	} else if (!strcmp(csp_inst_name, "p3c")) {
		n = 3;
		p_flow = (int*) p3c_flow;
		p_dist = p3c_dist;
	} else if (!strcmp(csp_inst_name, "p5")) {
		n = 5;
		p_flow = (int*) p5_flow;
		p_dist = p5_dist;
	} else if (!strcmp(csp_inst_name, "p10")) {
		n = 10;
		p_flow = (int*) p10_flow;
		p_dist = p10_dist;
	} else if (!strcmp(csp_inst_name, "choi15")) {
		n = 15;
		p_flow = (int*) choi15_flow;
		p_dist = choi15_dist;
	} else if (!strcmp(csp_inst_name, "esc16i")) {
		n = 16;
		p_flow = (int*) esc16i_flow;
		p_dist = esc16i_dist;
	} else if (!strcmp(csp_inst_name, "esc16e")) {
		n = 16;
		p_flow = (int*) esc16e_flow;
		p_dist = esc16e_dist;
	} else if (!strcmp(csp_inst_name, "esc16g")) {
		n = 16;
		p_flow = (int*) esc16g_flow;
		p_dist = esc16g_dist;
	} else {
		printf("\n%s QAP problem is not available. Please select an instance from \"csps/qap.c\" file.\n", csp_inst_name);
		return;
	}

	int cs[] = { n, 1, 1 };

	unsigned int* pos = malloc((unsigned int)n * sizeof(unsigned int));
	unsigned int* dist = malloc((unsigned int)n * (unsigned int)n * sizeof(unsigned int));

	// see if the problem is symmetric
	for (i = 0; i < n; ++i) {
		for (j = 0; j < i; ++j) {
			if (p_flow[i * n + j] != p_flow[j * n + i] || p_dist[i * n + j] != p_dist[j * n + i]) {
				symmetrical = false;
				break;
			}
		}
	}

	// see if the cost associated with the main diagonal is 0
	symmetrical_aux = symmetrical;

	if (symmetrical) {
		for (i = 0; i < n; ++i) {
			if (p_flow[i * n + i] != 0) {
				break;
			}
		}
		if (i < n) {
			for (i = 0; i < n; ++i) {
				if (p_dist[i * n + i] != 0) {
					break;
				}
			}
		}
		if (i < n) {
			symmetrical_aux = false;
		}
	}

	d_max = 0;
	for (i = 0; i < n * n; ++i) {
		d_max += (unsigned int)p_flow[i] * p_dist[i];
	}
	d_max++;
	if (d_max > 255) {
		d_max = 255;
	}

	for (i = 0; i < n; ++i) {
		pos[i] = v_new_range(0, (unsigned int)n - 1, true);
	}
	c_all_different(pos, (unsigned int)n);

	cost = v_new_range(0, d_max, false);

	for (i = 0; i < n * n; ++i) {
		dist[i] = v_new_val(p_dist[i]);
	}

	if (!symmetrical_aux) {
		new_dist = malloc((unsigned long)n * (unsigned long)n * sizeof(unsigned int));
		n_pos = malloc((unsigned long)n * (unsigned long)n * sizeof(unsigned int));

		for (i = 0; i < n * n; ++i) {
			new_dist[i] = v_new_range(0, d_max, false);
		}

		for (i = 0; i < n * n; ++i) {
			n_pos[i] = v_new_range(1, (unsigned int)n * (unsigned int)n, false);
		}

		vs[2] = v_new_val(1);	// element_var works with Y indexs from 1 to ... Because of Flatzinc specification
		for (i = 0; i < n; ++i) {
			vs[0] = pos[i];

			for (j = 0; j < n; ++j) {
				vs[1] = pos[j];

				c_linear_var(cs, vs, 3, n_pos[i * n + j]);
				c_element_var(dist, (unsigned int)n * (unsigned int)n, n_pos[(unsigned int)i * (unsigned int)n + (unsigned int)j], new_dist[(unsigned int)i * (unsigned int)n + (unsigned int)j]);
			}
		}
		c_linear_var(p_flow, new_dist, (unsigned int)n * (unsigned int)n, cost);

	} else {
		new_dist = malloc((unsigned long)n * (unsigned long)n * sizeof(unsigned int));
		n_pos = malloc((unsigned long)n * (unsigned long)n * sizeof(unsigned int));

		for (i = 1; i < n; ++i) {
			for (j = 0; j < i; ++j) {
				new_dist[i * n + j] = new_dist[j * n + i] = v_new_range(0, d_max, false);
				n_pos[(unsigned int)i * (unsigned int)n + (unsigned int)j] = n_pos[(unsigned int)j * (unsigned int)n + (unsigned int)i] = v_new_range(1, (unsigned int)n * (unsigned int)n, false);
			}
		}

		unsigned int val = 0;
		unsigned int v_id = v_new_val(val);
		for (i = 0; i < n; ++i) {
			new_dist[i * n + i] = v_id;
			n_pos[i * n + i] = v_id;
		}

		vs[2] = v_new_val(1);	// element_var works with Y indexs from 1 to ... Because of Flatzinc specification
		for (i = 1; i < n; ++i) {
			vs[0] = pos[i];

			for (j = 0; j < i; ++j) {
				vs[1] = pos[j];

				c_linear_var(cs, vs, 3, n_pos[i * n + j]);
				c_element_var(dist, (unsigned int)n * (unsigned int)n, n_pos[(unsigned int)i * (unsigned int)n + (unsigned int)j], new_dist[(unsigned int)i * (unsigned int)n + (unsigned int)j]);
			}
		}
		c_linear_var(p_flow, new_dist, (unsigned int)n * (unsigned int)n, cost);
	}

	// for optimization
	c_minimize(cost);

	if (OPTIMIZING) {
			printf("\nOptimizing %s Quadratic Assignment Problem with %d values.\n", csp_inst_name, n);
	} else if (FINDING_ONE_SOLUTION) {
		printf("\nFinding one solution for %s Quadratic Assignment Problem with %d values.\n", csp_inst_name, n);
	} else {
		printf("\nCounting all the solutions for %s Quadratic Assignment Problem with %d values.\n", csp_inst_name, n);
	}

	// Solve the CSP
	result = solve_CSP();

	if (OPTIMIZING && result == 1) {
		printf("Best solution:\n");
		vs_print_single_val(pos, (unsigned int)n, 1);
		printf("\n");
		v_print_cost(0);

	} else if (FINDING_ONE_SOLUTION && result == 1) {
		printf("Solution:\n");
		vs_print_single_val(pos, (unsigned int)n, 1);
		printf("\n");

	} else {
		printf("%lu solution(s) found\n", result);
	}

	free(pos);
	free(dist);
	free(n_pos);
	free(new_dist);
}