#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#include <time.h>

struct Set
{
  int n;    /* number of elements in the set */
  int min;  /* minimum distance between elements of this set */
  int best; /* best minimum distance for all sets up to this set */
};

/* Set sizes */
struct Set sets[] = {
  { 19 },
  { 19 },
  { 18 },
  { 16 },
  { 16 },
  { 14 },
  { 11 },
  { 11 },
  { 11 },
  { 9 },
  { 9 },
  { 7 },
  { 7 },
  { 7 },
  { 6 },
  { 5 },
  { 4 },
  { 3 },
  { 2 },
  { 2 },
};
#define Singles 6 /* number of one-element sets */

#define RT_BEGIN 10.        /* starting reciprocal "temperature" */
#define RT_END 1000.        /* ending reciprocal "temperature" */
#define RT_STEP 3.0e-5      /* speed of "temperature" change */
#define RT_FINE_STEP 1.0e-5 /* speed of "temperature" change on final step*/
#define IMPROVE_LIMIT 0.5/*portion of array values with harder qality function*/
#define SRAND 1             /* randomize */
#define SKIP_GROUPS 0       /* perform simulated annealing only once */

#define M (sizeof(sets) / sizeof(struct Set))
#define SinglesSet (-1)
#define INF INT_MAX

#define SWAP(a, b) { struct A t = a; a = b; b = t; }

struct A
{
  int set; /* set number */
  int l;   /* distance to the nearest element of the same set to the left */
  int r;   /* distance to the nearest element of the same set to the right */
};

struct Align
{
  struct A* a;    /* array */
  double* l;      /* pre-computed table of logarithms */
  int* m;         /* histograms of distances */
  int n;          /* number of elements in the array */
  int nn;         /* number of elements in largest sets */
  int bestMin;    /* optimal min.distance for largest sets */
  int startVar;   /* number of largest sets */
  int startWatch; /* start of not-yet optimized portion of sets */
};

/* Quality function */
static double eval(double* l, int n, int x)
{
  ++x;

  if (x > 0)
  {
    return (x >= n)? 0: l[x];
  }
  else
  {
    return -(1 << (1 - x));
  }
}

/* Calculates how much quality may be improved after swapping elements
   left and left+1 */
static double compare(struct Align* a, int left)
{
  double result = 0.;
  int right = left + 1;

  if (a->a[left].set != SinglesSet)
  {
    int best = sets[a->a[left].set].best;

    result +=
      eval(a->l, a->n, a->a[left].l + 1 - best) -
      eval(a->l, a->n, a->a[left].l - best) +
      eval(a->l, a->n, a->a[left].r - 1 - best) -
      eval(a->l, a->n, a->a[left].r - best);
  }

  if (a->a[right].set != SinglesSet)
  {
    int best = sets[a->a[right].set].best;

    result +=
      eval(a->l, a->n, a->a[right].l - 1 - best) -
      eval(a->l, a->n, a->a[right].l - best) +
      eval(a->l, a->n, a->a[right].r + 1 - best) -
      eval(a->l, a->n, a->a[right].r - best);
  }

  return result;
}

/* Harden quality function for not-yet-optimized portion of sets if possible */
static void improve(struct Align* a)
{
  int i;
  for (i = 0; i < a->startWatch; ++i)
  {
    if (sets[i].min < sets[i].best)
      return;
  }

  for (i = a->startWatch; i < M; ++i)
  {
    if (sets[i].min <= sets[i].best)
      return;
  }

  for (i = a->startWatch; i < M; ++i)
  {
    ++sets[i].best;
  }
}

/* Update distances for one of the swapped elements,
  update minimal distance for correspondent set */
static void
update(struct Align* a, int* d, int theSet, int pos, int diff, int lr)
{
  int better = 0;
  --a->m[a->n * theSet + *d];

  if (sets[theSet].min == *d &&
      a->m[a->n * theSet + *d] == 0)
  {
    sets[theSet].min += diff;
    better = (diff == 1)? 1: 0;
  }

  *d += diff;
  pos += *d * diff * lr * -1;
  ++a->m[a->n * theSet + *d];

  if (*d == sets[theSet].min - 1)
  {
    sets[theSet].min = *d;
  }

  if (diff * lr == 1) /* left */
  {
    a->a[pos].r += diff;
  }
  else /* right */
  {
    a->a[pos].l += diff;
  }

  if (better)
  {
    improve(a);
  }
}

/* Swap elements left and left+1 */
static void swap(struct Align* a, int left)
{
  int right = left + 1;
  SWAP(a->a[left], a->a[right])

  if (a->a[left].l < a->n)
    update(a, &(a->a[left].l), a->a[left].set, left, -1, -1);

  if (a->a[right].r < a->n)
    update(a, &(a->a[right].r), a->a[right].set, right, -1, 1);

  if (a->a[left].r < a->n)
    update(a, &(a->a[left].r), a->a[left].set, left, 1, -1);

  if (a->a[right].l < a->n)
    update(a, &(a->a[right].l), a->a[right].set, right, 1, 1);
}

/* Main loop of simulated annealing */
static void iterate(struct Align* a, double step)
{
  int i;
  double rt = RT_BEGIN;

  while (rt < RT_END)
  {
    /*printf("\nrt=%g\n", rt);*/

    for (i = 0; i != a->n; ++i)
    {
      double d;
      int left;

      do {
	left  = rand() % (a->n - 1);
      } while (a->a[left].set == a->a[left+1].set);

      d = compare(a, left);

      if (d >= 0. || exp(d * rt) * RAND_MAX > rand())
      {
	swap(a, left);
      }
    }

    rt *= 1.0 + step;
  }
}

/* Perform simulated annealing for equally-sized groups of sets,
  which allows to incrementally tighten quality function for these groups */
static void groups(struct Align* a)
{
  for (a->startWatch = a->startVar; a->startWatch < M;)
  {
    int oldN = sets[a->startWatch].n;
    iterate(a, RT_STEP);

#if SKIP_GROUPS
    return;
#endif

    for (; a->startWatch < M && sets[a->startWatch].n == oldN; ++a->startWatch)
    {
      a->nn += sets[a->startWatch].n;

      if (a->nn > a->n * IMPROVE_LIMIT)
      {
	goto exit;
      }
    }
  }

 exit:
  iterate(a, RT_FINE_STEP);
}

/* Allocate memory, initialize everything */
static void setup()
{
  int i;
  int theSet = 0;
  int theSetElems = sets[theSet].n;
  struct Align a;

  a.a = 0;
  a.l = 0;
  a.m = 0;
  a.n = Singles;
  a.startVar = 0;

  for (i = 0; i < M; ++i)
  {
    a.n += sets[i].n;

    if (!a.startVar && sets[i].n != sets[0].n)
    {
      a.startVar = i;
    }

    if (!a.startVar)
    {
      a.nn += sets[i].n;
    }
  }

  a.a = calloc(a.n, sizeof(struct A));
  a.l = calloc(a.n, sizeof(double));
  a.m = calloc(a.n, M * sizeof(int));

  for (i = 1; i < a.n; ++i)
  {
    a.l[i] = log(i);
  }

  a.bestMin = (a.n - a.startVar) / (sets[0].n - 1);
  i = 0;
  while (1)
  {
    a.a[i].set = theSet;
    a.a[i].l = (theSetElems == sets[theSet].n)? INF: 1;
    --theSetElems;
    a.a[i].r = theSetElems? 1: INF;
    ++i;

    if (!theSetElems)
    {
      sets[theSet].min = 1;
      sets[theSet].best = a.bestMin;
      a.m[a.n * theSet + 1] = sets[theSet].n - 1;
      ++theSet;

      if (theSet == M)
	break;

      theSetElems = sets[theSet].n;
    }
  }
  for (; i < a.n; ++i)
  {
    a.a[i].set = SinglesSet;
    a.a[i].l = INF;
    a.a[i].r = INF;
  }

  groups(&a);

  /*for (i = 0; i < a.n; ++i)
  {
    printf("%d ", a.a[i].set);
  }*/

  free(a.m);
  free(a.l);
  free(a.a);
}

int main()
{
#if SRAND
  srand(time(0));
#endif

  setup();

  /* Print set size, minimal distance, limit for minimal distance */
  int i;
  for (i = 0; i < M; ++i)
  {
    printf("\n%d\t%d", sets[i].n, sets[i].min);

    if (sets[i].best != sets[i].min)
    {
      printf("\t%d", sets[i].best);
    }
  }
  puts("\n");

  return 0;
}