/*
  Quadratic Assignment Problem

  computes min( \Sum_{l,c} distance(perm(l), perm(c)) * flow(l, c) )

  (with --swap, swaps the flow and distance matrixes, working with the
  dual problem)
*/

#include <stdio.h>
#include <stdlib.h>
#include <fcntl.h>
#include <errno.h>
#include <ctype.h>
#include <string.h>
#include <alloca.h>

#ifdef SPLITGO_MPI
#  include <mpi.h>
#endif

#include "fdc_int.h"

// try to save on variables
static fd_int fd_const[MAX_VALUE + 1];

#define FD_CONST(n) (fd_const[n] ? fd_const[n] : (fd_const[n] = fd_new(n, n)))

//#define FD_CONST(n) fd_new(n, n)

#define NMAX 128

typedef struct {
  int size;
  int *flow, *dist;
  char *name;
  int opt;
  int *sol;
} qap;

// args: problem name, size, optimal cost, solution (permutation)
#define QAP(p, s, ...) \
  qap p = { s, p ## _flow, p ## _dist, #p, ## __VA_ARGS__ }

#include "qap.h"	// problems

#define Pflow(a,b) (P->flow[(a) * N + (b)])
#define Pdist(a,b) (P->dist[(a) * N + (b)])
#define Psol(a)    (P->sol[a])

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

QAP(p2, 2);

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

QAP(p2c, 2);

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

QAP(p3, 3);

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

QAP(p3b, 3);

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

QAP(p3c, 3);

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
};

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
};

QAP(p5, 5);

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
};

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
};

int p10_sol[] = { 5, 1, 2, 3, 4, 8, 9, 6, 10, 7 }; // ?
		// 4 1 2 3 5 8 10 6 9 7
		// 1 2 3 4 5 9 10 6 8 7

QAP(p10, 10, 58 /* ? */);

/*
  sunfire: 4 - 6:38 (0.10); 2 - 13:16 (0.13); 1 - 26:26 (0.13) (symmetric)
	       3.98		 1.99

  diaz:    8 - 1:44 (0.02); 4 - 3:26 (0.02); 2 - 6:41 (0.02); 1 - 13:17 (0.02)
	       7.68		3.86		 1.99		     (symmetric)

  ism:     4 - 4:48 (0.03); 2 - 9:37 (0.03); 1 - 19:10 (0.03) (symmetric)
	       4.00		1.99

  bicho:   1:11:22; solution: 4 - 0.58; 1 - 0.26 (symmetric)

  qap.pl:  12:49 (0.07)
*/

static qap *findqap(char *s)
{
  static qap *qaps[] = {
    &p2, &p2c, &p3, &p3b, &p3c, &p5, &p10,
    &esc16a, &esc16d, &esc16e, &esc16f, &esc16g, &esc16h, &esc16i, &esc16j,
    &esc32a, &esc32e, &esc64a, &esc128,
    &choi16, &choi15,
    0
  };
  int i = 0;

  while (qaps[i] && strcmp(qaps[i]->name, s))
    ++i;

  return qaps[i];
}

static qap *readqap(char *s)
{
  FILE *f;
  qap *p;
  int i;
  int _;	// nuke gcc warnings

  if ((f = fopen(s, "r")) == NULL)
    return 0;

  p = malloc(sizeof(*p));

  p->name = s;

  // read the problem size
  _ = fscanf(f, "%d", &p->size);

  // worry about the suitability problem later

  // skip what else may be on the same line
  while (getc(f) != '\n')
    ;

  p->opt = -1; p->sol = NULL;

  // allocate memory for the matrixes
  p->flow = malloc(p->size * p->size * sizeof(*p->flow));
  p->dist = malloc(p->size * p->size * sizeof(*p->dist));

  // and read them
  for (i = 0; i < p->size * p->size; ++i)
    _ = fscanf(f, "%d", p->flow + i);

  for (i = 0; i < p->size * p->size; ++i)
    _ = fscanf(f, "%d", p->dist + i);

  return p;
}

int main(int argc, char *argv[])
{
  int solutions = 0, one_solution = 1;
  qap *P = 0;
  int N;
  fd_int perm[NMAX][NMAX], sums[NMAX], prods[2][NMAX * NMAX], cost;
  fd_int ndist[NMAX][NMAX], pos[NMAX];
  fd_int idist[NMAX][NMAX], distt[NMAX][NMAX];
  fd_int tv[NMAX];
  int ti[NMAX];
  int l, c, i, j;
  int symmetrical = 1;
  int symmetrical0;  // symmetric problem and at least one of the diagonals is 0
  int generic = 0;   // whether to specialise for symmetric problems
  int swap = 0;      // swap flow and distance matrixes
  int all_diff = 0;  // use an all-different constraint (redundant)
  int use_label = 0; // restrict labelling to the permutation variables
  int fixed_cost = 0;
  int cost_min = MIN_VALUE, cost_max = MAX_VALUE;

  fd_init(&argc, &argv);

  for (j = 1; j < argc; ++j)
    if (!strcmp(argv[j], "--all"))
      one_solution = 0;
    else if (!strcmp(argv[j], "--gen"))		// ignore symmetries
      generic = 1;
    else if (!strcmp(argv[j], "--swap"))	// swap distance & flow matrixes
      swap = 1;
    else if (!strcmp(argv[j], "--diff"))	// use all-different constraint
      all_diff = 1;
    else if (!strcmp(argv[j], "--label"))	// restricted labelling
      use_label = 1;
    else if (!strncmp(argv[j], "--fix", strlen("--fix")))
      {
	fixed_cost = 1;

	if (argv[j][strlen("--fix")] == '=')
	  cost_min = cost_max = atoi(argv[j] + strlen("--fix") + 1);
      }
    else if (!strncmp(argv[j], "--cost=", strlen("--cost=")))
      {
	char *a = argv[j] + strlen("--cost=");

	cost_min = atoi(a);

	if (a = strchr(a, '-'))
	  cost_max = atoi(a + 1);
      }
    else if (isdigit(*argv[j]))
      break;
    else if (!P && ((P = findqap(argv[j])) || (P = readqap(argv[j]))))
      ;
    else
      {
	_fd_error("%s: invalid argument `%s'\n", argv[0], argv[j]);

	return 2;
      }

  if (!P)
    P = findqap("esc16i");

  N = P->size;

  if (N > NMAX)
    _fd_fatal("qap: problem too big (increase NMAX)");
#ifdef QAP_MODEL_B
  if (N * N > DOMAIN_BITS)
    _fd_fatal("qap: problem too big (model B) (increase DOMAIN_BITS)");
#endif

  if (swap)
    {
      int *p = P->flow; P->flow = P->dist; P->dist = p;
    }

  if (P->name)
    _fd_debug("%s\n", P->name);

  {
    int diag_cost = 0, tb = 0;

    for (i = 0; i < N; ++i)
      diag_cost += Pflow(i, i) * Pdist(i, i);

    for (l = 0; l < N; ++l)
      for (c = 0; c < N; ++c)
	tb += Pflow(l, c) * Pdist(l, c);

    _fd_debug("tight bound %d (diagonal %d)\n", tb, diag_cost);
  }

  // see if the problem is symmetric
  for (l = 0; l < N; ++l)
    for (c = 0; c < l; ++c)
      if (Pflow(l, c) != Pflow(c, l) || Pdist(l, c) != Pdist(c, l))
	{
	  symmetrical = 0;

	  break;
	}

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

  if (symmetrical)
    {
      for (i = 0; i < N; ++i)
	if (Pflow(i, i) != 0)
	  break;

      if (i < N)
	for (i = 0; i < N; ++i)
	  if (Pdist(i, i) != 0)
	    break;

      if (i < N)
	symmetrical0 = 0;
    }

  if (symmetrical)
    if (symmetrical0)
      _fd_debug("symmetrical problem with 0 diagonal\n");
    else
      _fd_debug("symmetrical problem\n");
  else
    _fd_debug("asymmetrical problem (%d, %d)\n", l, c);

  if (generic)
    symmetrical0 = 0;

  // location of the facilities (permutation)
  // (used to serve only for visualisation, but may be a more efective
  // way of generating the permutations)
#if 01
  for (i = 0; j < argc; ++i, ++j)
    {
      int lb, ub;
      char *s;

      lb = ub = atoi(argv[j]) - 1;

      if (s = strchr(argv[j], '-'))
	ub = atoi(s + 1) - 1;

      pos[i] = fd_new(lb, ub);
    }

  for (; i < N; ++i)
    pos[i] = fd_new(0, N - 1);
#else
  for (i = 0; i < N; ++i)
    pos[i] = fd_new(Psol(i) - 1, Psol(i) - 1);
#endif

  if (use_label)
    fd_label(pos, N);

  // cost variable
  cost = fd_new(cost_min, cost_max);

#ifndef QAP_MODEL_B
#ifdef SPLITGO_MPI
  {
    int rank;

    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    if (rank == 0)
      printf("model A\n");
  }
#else
  printf("model A\n");
#endif

  if (all_diff)
    fd_all_different(pos, N);

  // permutation matrix
  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      perm[l][c] = fd_new(0, 1);

  // there can only be one 1 per row ...
  for (l = 0; l < N; ++l)
    fd_exactly(perm[l], N, 1, 1);

  // ... and per column
  for (c = 0; c < N; ++c)
    {
      for (l = 0; l < N; ++l)
	tv[l] = perm[l][c];

      fd_exactly(tv, N, 1, 1);
    }

  // the position of the one within row i corresponds to the
  // permutated location of facility i
  for (i = 0; i < N; ++i)
    fd_element(perm[i], N, pos[i], 1);

  // permutated distances = perm x dist x perm^T
  // calculate idist = perm x dist
  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      idist[l][c] = fd_new(0, MAX_VALUE);

  // the *transpose* of the distance matrix
  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      distt[l][c] = FD_CONST(Pdist(c, l));

  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      fd_knapsack2(perm[l], distt[c], N, idist[l][c]);

  // calculate ndist = idist x perm^T
  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      ndist[l][c] = fd_new(0, MAX_VALUE);

  for (l = 0; l < N; ++l)
    for (c = 0; c < N; ++c)
      fd_knapsack2(idist[l], perm[c], N, ndist[l][c]);

  // costs

  if (!symmetrical0)
    {
#if 0
      for (i = 0; i < N; ++i)
	sums[i] = fd_new(0, MAX_VALUE);

      for (l = 0; l < N; ++l)
	for (c = 0; c < N; ++c)
	  distt[l][c] = FD_CONST(Pflow(l, c));

      for (l = 0; l < N; ++l)
	fd_knapsack2(ndist[l], distt[l], N, sums[l]);

      fd_sum2(sums, N, cost);
#elif 0
      i = 0;

      for (l = 0; l < N; ++l)
	for (c = 0; c < N; ++c)
	  {
	    prods[0][i] = ndist[l][c];
	    prods[1][i] = FD_CONST(Pflow(l, c));

	    ++i;
	  }

      fd_knapsack2(prods[0], prods[1], i, cost);
#else
      i = 0;

      for (l = 0; l < N; ++l)
	for (c = 0; c < N; ++c)
	  {
	    prods[0][i] = ndist[l][c];

	    ++i;
	  }

      fd_poly_eq(P->flow, prods[0], N * N, cost);
#endif
    }
  else
    {
#if 0
      for (i = 0; i < N - 1; ++i)
	sums[i] = fd_new(0, MAX_VALUE);

      for (l = 1; l < N; ++l)
	for (c = 0; c < l; ++c)
	  distt[l][c] = FD_CONST(Pflow(l, c));

      for (l = 1; l < N; ++l)
	fd_knapsack2(ndist[l], distt[l], l, sums[l - 1]);

      fd_sum2(sums, N - 1, cost);
#elif 01
      i = 0;

      for (l = 1; l < N; ++l)
	for (c = 0; c < l; ++c)
	  {
	    prods[0][i] = ndist[l][c];
	    prods[1][i] = FD_CONST(Pflow(l, c));

	    ++i;
	  }

      fd_knapsack2(prods[0], prods[1], i, cost);
#else
      fd_int ps[NMAX * NMAX];

      i = 0;

      for (l = 1; l < N; ++l)
	for (c = 0; c < l; ++c)
	  {
	    prods[0][i] = ndist[l][c];
	    prods[1][i] = FD_CONST(Pflow(l, c));

	    ++i;
	  }

      for (j = 0; j < i; ++j)
	{
	  ps[j] = fd_new(0, MAX_VALUE);

	  fd_var_eq_times(ps[j], prods[0][j], prods[1][j]);
	}

      fd_sum2(ps, i, cost);
#endif
    }
#else /* !QAP_MODEL_B */
  fd_int dist[NMAX * NMAX];
  fd_int *new_dist;		// permutated distance matrix
  fd_int *npos;			// new_dist[npos[(l,c)]] = dist[(l,c)]
  fd_int vs[2];			// for computing p' = l' * N + c'
  int cs[] = { 0, 1 };		// 		   "

#ifdef SPLITGO_MPI
  {
    int rank;

    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    if (rank == 0)
      printf("model B\n");
  }
#else
  printf("model B\n");
#endif

  fd_all_different(pos, N);		// required in this model

  cs[0] = N;

  // distances
  for (i = 0; i < N * N; ++i)
    dist[i] = FD_CONST(P->dist[i]);

  // costs

  if (!symmetrical0)
    {
      new_dist = alloca(N * N * sizeof(*new_dist));
      npos = alloca(N * N * sizeof(*npos));

      for (i = 0; i < N * N; ++i)
	new_dist[i] = fd_new(0, MAX_VALUE);

      for (i = 0; i < N * N; ++i)
	npos[i] = fd_new(0, N * N - 1);

      for (l = 0; l < N; ++l)
	{
	  vs[0] = pos[l];

	  for (c = 0; c < N; ++c)
	    {
	      vs[1] = pos[c];

	      fd_poly_eq(cs, vs, 2, npos[l * N + c]);

	      fd_element_var(dist, N * N, npos[l * N + c], new_dist[l * N + c]);
	    }
	}

      fd_poly_eq(P->flow, new_dist, N * N, cost);
    }
  else
    {
      new_dist = alloca(N * N * sizeof(*new_dist));
      npos = alloca(N * N * sizeof(*npos));

      for (l = 1; l < N; ++l)
	for (c = 0; c < l; ++c) 
	{
	  new_dist[l * N + c] = new_dist[c * N + l] = fd_new(0, MAX_VALUE);
	  npos[l * N + c] = npos[c * N + l] = fd_new(0, N * N - 1);
	}

      for (l = 0; l < N; ++l)
	{
	  new_dist[l * N + l] = FD_CONST(0);
	  npos[l * N + l] = FD_CONST(0);
	}

      for (l = 1; l < N; ++l)
	{
	  vs[0] = pos[l];

	  for (c = 0; c < l; ++c)
	    {
	      vs[1] = pos[c];

	      fd_poly_eq(cs, vs, 2, npos[l * N + c]);

	      fd_element_var(dist, N * N, npos[l * N + c], new_dist[l * N + c]);
	    }
	}

      fd_poly_eq(P->flow, new_dist, N * N, cost);
    }  
#endif /* !QAP_MODEL_B */

  if (!fixed_cost)
    fd_min(cost);

  while (fd_solve() == FD_OK)
    {
      printf("solution %d:\n", ++solutions);

      if (P->name)
	printf("%s\n", P->name);

#ifndef QAP_MODEL_B
      printf("cost = %d\n", (symmetrical0 ? 2 : 1) * fd_var_value(cost));
#else
      printf("cost = %d\n", fd_var_value(cost));
#endif

      putchar('(');
      for (i = 0; i < N; ++i)
	{
	  printf("%d", fd_var_value(pos[i]) + 1);

	  if (i != N - 1)
	    putchar(' ');
	}
      printf(")\n");

#if !(defined(LOCAL_SEARCH) || defined(DISTRIBUTED_SOLVER))
      if (one_solution)
#endif
	break;
    }

  if (solutions)
    printf("%d solutions found\n", solutions);
  else
    {
      if (P->name)
	printf("%s\n", P->name);

      printf("inconsistent CSP\n");
    }

  fd_end();

  return !solutions;
}