/* Balanced Incomplete Block Design (BIBD) (à la Gecode) */

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

#include "fdc_int.h"

#define MAX_V MAX_VARIABLES
#define MAX_B MAX_VARIABLES

int main(int argc, char *argv[])
{
  int V = 7;	  // no. of distinct objects
  int B = 7;	  // no. of blocks (groups)
  int K = 3;	  // objects per block
  int R = 3;	  // blocks per object (no. of blocks an object appears in)
  int lambda = 1; // blocks per 2-objects (no. of blocks a pair of
		  // objects appears in)

  fd_int *vs;
#define vs(r, c) vs[(r) * B + (c)]
  fd_int *ts, *us;
  int v, b, v0, v1;
  int solutions = 0, one_solution = 1;
  int i;

  fd_init(&argc, &argv);

  for (i = 1; i < argc; ++i)
    if (!strcmp(argv[i], "--all"))
      one_solution = 0;
#ifndef GECODE_STYLE
    else if (argc - i < 3)
      {	       
	_fd_debug("must give V, B, and K (or nothing)\n");
	return 1;
      }
    else
      {
	V = atoi(argv[i++]);
	B = atoi(argv[i++]);
	K = atoi(argv[i]);
      }

  R = K * B / V;	// must be an integer

  if (R * V != K * B)
    {
      _fd_debug("%d * %d / %d must be an integer\n", K, B, V);
      return 1;
    }

  lambda = R * (K - 1) / (V - 1);	// idem

  if (lambda * (V - 1) != R * (K - 1))
    {
      _fd_debug("%d * (%d - 1) / (%d - 1) must be an integer\n", R, K, V);
      return 1;
    }
#else /* GECODE_STYLE */
    else if (argc - i < 3)
      {	       
	_fd_debug("must give V, K, and lambda (or nothing)\n");
	return 1;
      }
    else
      {
	V = atoi(argv[i++]);
	K = atoi(argv[i++]);
	lambda = atoi(argv[i]);
      }

  B = lambda * V * (V - 1) / K / (K - 1);

  if (B * K * (K - 1) != lambda * V * (V - 1))
    {
      _fd_debug("warning: %d * %d * (%d - 1) / %d / (%d - 1) must be an integer\n",
		lambda, V, V, K, K);
//      return 1;
    }

  R = K * B / V;	// must be an integer

  if (R * V != K * B)
    {
      _fd_debug("warning: %d * %d / %d must be an integer\n", K, B, V);
//      return 1;
    }
#endif /* GECODE_STYLE */

  _fd_debug("V = %d, B = %d, K = %d, R = %d, lambda = %d\n", V, B, K, R, lambda);

  vs = calloc(V * B, sizeof(fd_int));
  ts = calloc(V > B ? V : B, sizeof(fd_int));
  us = calloc(B, sizeof(fd_int));

  for (v = 0; v < V; ++v)
    for (b = 0; b < B; ++b)
      vs(v, b) = fd_new(0, 1);

  // an object appears in R blocks
  for (v = 0; v < V; ++v)
    fd_sum(vs + B * v, B, R);

  // K objects per block
  for (b = 0; b < B; ++b)
    {
      for (v = 0; v < V; ++v)
	ts[v] = vs(v, b);

      fd_sum(ts, V, K);
    }

  // the scalar (or dot) product between any two rows must be lambda
#if 0
  for (v0 = 0; v0 < V; ++v0)
    for (v1 = v0 + 1; v1 < V; ++v1)
      {
	for (b = 0; b < B; ++b)
	  {
	    ts[b] = fd_new(0, 1);

	    fd_var_eq_times(ts[b], vs(v0, b), vs(v1, b));
	  }

	fd_sum(ts, B, lambda);
      }
#else
  for (v0 = 0; v0 < V; ++v0)
    for (v1 = v0 + 1; v1 < V; ++v1)
      {
	for (b = 0; b < B; ++b)
	  {
	    ts[b] = vs(v0, b);
	    us[b] = vs(v1, b);
	  }

	fd_sum_prod(ts, us, B, lambda);
      }
#endif

#if 01
  // try to eliminate solutions which result from exchanging rows
  // and/or columns
  // XXX: quite incomplete

  // constrain the first row to be of the form 1 ... 1 0 ... 0
  for (b = 1; b < B; ++b)
    fd_ge(vs(0, b - 1), vs(0, b));

  // constrain the second row to be of the form 1 ... 1 0 ... 0 1 ... 1 0 ... 0
  // XXX: could constrain vs[1,R-1] < vs[1,R]
  for (b = 1; b < R; ++b)
    fd_ge(vs(1, b - 1), vs(1, b));
  for (b = R + 1; b < B; ++b)
    fd_ge(vs(1, b - 1), vs(1, b));

  // constrain the first column to be of the form 1 ... 1 0 ... 0
  for (v = 1; v < V; ++v)
    fd_ge(vs(v - 1, 0), vs(v, 0));

  // constrain the second column to be of the form 1 ... 1 0 ... 0 1 ... 1 0 ... 0
  // XXX: could constrain vs[K-1,1] < vs[K,1]
  for (v = 1; v < K; ++v)
    fd_ge(vs(v - 1, 1), vs(v, 1));
  for (v = K + 1; v < V; ++v)
    fd_ge(vs(v - 1, 1), vs(v, 1));
#endif

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

      for (v = 0; v < V; ++v)
	{
	  for (b = 0; b < B; ++b)
	    {
	      fd_print(vs(v, b));
	      putchar(' ');
	    }
	  putchar('\n');
	}

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

  if (solutions)
    printf("%d solutions found\n", solutions);
  else
    printf("inconsistent CSP\n");

  fd_end();

  return !solutions;
}