/*
    This program reads a file pbn.dat (which contains
    a description of a paint-by-numbers puzzle) and 
    then constructs a file pbn.lp (a CPLEX ".lp" file
    that describes an integer program that can be used
    to solve it).

    It was written by Robert Bosch in September 2000
    to accompany the "Painting by Numbers" installment 
    of his Mindsharpener column, which is scheduled 
    to appear in issue 65 of Optima, the newsletter
    of the Mathematical Programming Society.  It was
    modified in January 2001 to make the code's variable 
    names consistent with the notation used in the
    article.  To access the electronic version of 
    Optima, point your browser to:

      http://www.ise.ufl.edu/~optima/

    Please send comments or questions to: 
	
      Robert Bosch
      Dept. of Mathematics
      Oberlin College
      Oberlin, Ohio  44074

      bobb@cs.oberlin.edu
      http://www.oberlin.edu/~math/people/faculty/bosch.html
 */

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

#define MAX_M 30   /* maximum # of rows     */
#define MAX_N 30   /* maximum # of columns  */
#define MAX_K 11   /* maximum # of clusters */

main()
{

int m, n;
int i, ip, j, jp, t;
int sum;

int kr[MAX_M+1];           /* kr[i] = # clusters in row i             */
int sr[MAX_M+1][MAX_K+1];  /* sr[i][t] = size of row i's tth cluster  */
int er[MAX_M+1][MAX_K+1];  /* er[i][t] = smallest j such that row i's */
                           /*  tth cluster can be placed with its     */
                           /*  leftmost pixel occupying column j      */
int lr[MAX_M+1][MAX_K+1];  /* lr[i][t] = largest j such that row i's  */
                           /*  tth cluster can be placed with its     */
                           /*  leftmost pixel occupying column j      */

int kc[MAX_N+1];           /* kc[j] = # clusters in column j          */
int sc[MAX_N+1][MAX_K+1];  /* sc[j][t] = size of col j's tth cluster  */
int ec[MAX_N+1][MAX_K+1];  /* ec[j][t] = smallest i such that column  */
                           /*  j's tth cluster can be placed with its */
                           /*  uppermost pixel occupying row i        */
int lc[MAX_N+1][MAX_K+1];  /* lc[j][t] = largest i such that column   */
                           /*  j's tth cluster can be placed with its */
                           /*  uppermost pixel occupying row i        */

int fix[MAX_M+1][MAX_N+1]; /* fix[i][j] = 1 if Zij should be set = 1  */

char trash_string[25];

FILE *input, *lp_file;

/*
    Read the file describing the paint-by-numbers puzzle. 
 */

if ((input = fopen("pbn.dat", "r")) == NULL)
  {
    printf("Couldn't open pbn.dat!\n");
    exit(0);
  }

fscanf(input, "%s", trash_string);
fscanf(input, "%d", &m);
fscanf(input, "%s", trash_string);
fscanf(input, "%d", &n);

sum = 0;

for (i = 1; i <= m; i++)
  {
    fscanf(input, "%s", trash_string);  /*  row_2:                 */
    fscanf(input, "%s", trash_string);  /*  number_of_clusters: 4  */
    fscanf(input, "%d", &kr[i]);
    if (kr[i] > 0)
      {
        fscanf(input, "%s", trash_string);  /*  size(s): 1 2 1 1   */
        for (t = 1; t <= kr[i]; t++)
          {
            fscanf(input, "%d", &sr[i][t]);
            sum += sr[i][t];
          }
      }
  }

for (j = 1; j <= n; j++)
  {
    fscanf(input, "%s", trash_string);  /*  col_2:                 */
    fscanf(input, "%s", trash_string);  /*  number_of_clusters: 3  */
    fscanf(input, "%d", &kc[j]);
    if (kc[j] > 0)
      {
        fscanf(input, "%s", trash_string);  /*  size(s): 1 5 2     */
        for (t = 1; t <= kc[j]; t++)
          fscanf(input, "%d", &sc[j][t]);
      }
  }

fclose(input);

/*
    For each cluster, compute the "earliest" and "latest" values.
 */

for (i = 1; i <= m; i++)
  if (kr[i] > 0)
    {
      er[i][1] = 1;
      for (t = 2; t <= kr[i]; t++)
        er[i][t] = er[i][t-1] + sr[i][t-1] + 1;
    }

for (i = 1; i <= m; i++)
  if (kr[i] > 0)
    {
      lr[i][kr[i]] = n + 1 - sr[i][kr[i]];
      for (t = kr[i]-1; t >= 1; t--)
        lr[i][t] = lr[i][t+1] - sr[i][t] - 1;
    }

for (j = 1; j <= n; j++)
  if (kc[j] > 0)
    {
      ec[j][1] = 1;
      for (t = 2; t <= kc[j]; t++)
        ec[j][t] = ec[j][t-1] + sc[j][t-1] + 1;
    }

for (j = 1; j <= n; j++)
  if (kc[j] > 0)
    {
      lc[j][kc[j]] = m + 1 - sc[j][kc[j]];
      for (t = kc[j]-1; t >= 1; t--)
        lc[j][t] = lc[j][t+1] - sc[j][t] - 1;
    }

/*
    Write the lp file.
 */

if ((lp_file = fopen("pbn.lp", "w")) == NULL)
  {
    printf("Couldn't open pbn.lp!\n");
    exit(0);
  }

/*
    There's no need for an objective function.
    I decided to include one anyway.
 */

fprintf(lp_file, "minimize\n");

for (i = 1; i <= m; i++)
  for (j = 1; j <= n; j++)
    fprintf(lp_file, " + Z%d,%d\n", i, j);

fprintf(lp_file, "\n");

fprintf(lp_file, "subject to\n");

/*
    Row i's tth cluster must appear in row i exactly once.
 */

for (i = 1; i <= m; i++)
  if (kr[i] > 0)
    for (t = 1; t <= kr[i]; t++)
      {
        for (j = er[i][t]; j <= lr[i][t]; j++)
          fprintf(lp_file, " + Y%d,%d,%d", i, t, j);
        fprintf(lp_file, " = 1\n");
      }
 
fprintf(lp_file, "\n");

/*
    Row i's (t+1)th cluster must be placed to the right
    of its tth cluster.
 */

for (i = 1; i <= m; i++)
  if (kr[i] > 1)
    for (t = 1; t <= kr[i]-1; t++)
      for (j = er[i][t]; j <= lr[i][t]; j++)
        {
          fprintf(lp_file, " + Y%d,%d,%d", i, t, j);
          for (jp = j + sr[i][t] + 1; jp <= lr[i][t+1]; jp++)
            fprintf(lp_file, " - Y%d,%d,%d", i, t+1, jp);
          fprintf(lp_file, " <= 0\n");
        }

fprintf(lp_file, "\n");

/*
    Column j's tth cluster must appear in column j exactly once.
 */

for (j = 1; j <= n; j++)
  if (kc[j] > 0)
    for (t = 1; t <= kc[j]; t++)
      {
        for (i = ec[j][t]; i <= lc[j][t]; i++)
          fprintf(lp_file, " + X%d,%d,%d", j, t, i);
        fprintf(lp_file, " = 1\n");
      }

fprintf(lp_file, "\n");
 
/*
    Column j's (t+1)th cluster must be placed below its
    tth cluster.
 */

for (j = 1; j <= n; j++)
  if (kc[j] > 1)
    for (t = 1; t <= kc[j]-1; t++)
      for (i = ec[j][t]; i <= lc[j][t]; i++)
        {
          fprintf(lp_file, " + X%d,%d,%d", j, t, i);
          for (ip = i + sc[j][t] + 1; ip <= lc[j][t+1]; ip++)
            fprintf(lp_file, " - X%d,%d,%d", j, t+1, ip);
          fprintf(lp_file, " <= 0\n");
        }

fprintf(lp_file, "\n");

/*
    The remaining constraints are "double coverage" constraints.
 */

for (i = 1; i <= m; i++)
  for (j = 1; j <= n; j++)
    if ((kr[i] > 0) && (kc[j] > 0))
      {
        fprintf(lp_file, " + Z%d,%d", i, j);
        for (t = 1; t <= kr[i]; t++)
          for (jp = j - sr[i][t] + 1; jp <= j; jp++)
            if ((jp >= er[i][t]) && (jp <= lr[i][t]))
              fprintf(lp_file, " - Y%d,%d,%d", i, t, jp);
        fprintf(lp_file, " <= 0\n");

        fprintf(lp_file, " + Z%d,%d", i, j);
        for (t = 1; t <= kc[j]; t++)
          for (ip = i - sc[j][t] + 1; ip <= i; ip++)
            if ((ip >= ec[j][t]) && (ip <= lc[j][t]))
              fprintf(lp_file, " - X%d,%d,%d", j, t, ip);
        fprintf(lp_file, " <= 0\n");
      }

fprintf(lp_file, "\n");
 
for (i = 1; i <= m; i++)
  for (j = 1; j <= n; j++)
    for (t = 1; t <= kr[i]; t++)
      for (jp = j - sr[i][t] + 1; jp <= j; jp++)
        if ((jp >= er[i][t]) && (jp <= lr[i][t]))
          {
            fprintf(lp_file, " + Z%d,%d", i, j);
            fprintf(lp_file, " - Y%d,%d,%d >= 0\n", i, t, jp);
          }

fprintf(lp_file, "\n");

for (i = 1; i <= m; i++) 
  for (j = 1; j <= n; j++)
    for (t = 1; t <= kc[j]; t++)
      for (ip = ec[j][t]; ip <= lc[j][t]; ip++)
        if ((ip >= i - sc[j][t] + 1) && (ip <= i))
          {
            fprintf(lp_file, " + Z%d,%d", i, j);
            fprintf(lp_file, " - X%d,%d,%d >= 0\n", j, t, ip);
          }

fprintf(lp_file, "\n");

/*
    Determine which Zij can be set = 1.
 */

for (i = 1; i <= m; i++)
  for (j = 1; j <= n; j++)
    fix[i][j] = -1;

for (i = 1; i <= m; i++)
  for (t = 1; t <= kr[i]; t++)
    for (j = lr[i][t]; j <= er[i][t] + sr[i][t] - 1; j++)
      fix[i][j] = 1;

for (j = 1; j <= n; j++)
  for (t = 1; t <= kc[j]; t++)
    for (i = lc[j][t]; i <= ec[j][t] + sc[j][t] - 1; i++)
      fix[i][j] = 1;

fprintf(lp_file, "bounds\n");

for (i = 1; i <= m; i++)
  for (j = 1; j <= n; j++)
    if (fix[i][j] == 1)
      fprintf(lp_file, " Z%d,%d = 1\n", i, j);
    else
      fprintf(lp_file, " 0 <= Z%d,%d <= 1\n", i, j);

fprintf(lp_file, "\n");
 
for (i = 1; i <= m; i++)
  if (kr[i] > 0)
    for (t = 1; t <= kr[i]; t++)
      for (j = er[i][t]; j <= lr[i][t]; j++)
        fprintf(lp_file, " 0 <= Y%d,%d,%d <= 1\n", i, t, j);

fprintf(lp_file, "\n");
 
for (j = 1; j <= n; j++)
  if (kc[j] > 0)
    for (t = 1; t <= kc[j]; t++)
      for (i = ec[j][t]; i <= lc[j][t]; i++)
        fprintf(lp_file, " 0 <= X%d,%d,%d <= 1\n", j, t, i);

fprintf(lp_file, "\n");

fprintf(lp_file, "integers\n");

for (i = 1; i <= m; i++)
  for (j = 1; j <= n; j++)
    fprintf(lp_file, " Z%d,%d\n", i, j);

fprintf(lp_file, "\n");
 
for (i = 1; i <= m; i++)
  if (kr[i] > 0)
    for (t = 1; t <= kr[i]; t++)
      for (j = er[i][t]; j <= lr[i][t]; j++)
        fprintf(lp_file, " Y%d,%d,%d\n", i, t, j);

fprintf(lp_file, "\n");
 
for (j = 1; j <= n; j++)
  if (kc[j] > 0)
    for (t = 1; t <= kc[j]; t++)
      for (i = ec[j][t]; i <= lc[j][t]; i++)
        fprintf(lp_file, " X%d,%d,%d\n", j, t, i);

fprintf(lp_file, "\n");

fprintf(lp_file, "end\n");

fclose(lp_file);

}