ex4_enum.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <assert.h>
#include <stdbool.h>
#include <time.h>

#include "mshell.h"  // or by allocate?
#include "sid.h"
#include "allocate.h"

#define BIP_SID  0x42495078
#define BIP_COLS 32
#define BIP_ROWS 32

struct binary_program
{
   SID
   int    rows;                   
   int    cols;                   
   double a[BIP_ROWS][BIP_COLS];  
   double b[BIP_ROWS];            
};

typedef struct binary_program BIP;

#define MAX_LINE_LEN   512  

#define MAX_COEF_VAL   1e20 

#define GET_SEC(a, b)  ((b - a) / (double)CLOCKS_PER_SEC) 

typedef enum { READ_COLS, READ_ROWS, READ_COEF } LINE_MODE;

bool bip_is_valid(const BIP* bip)
{
   if (NULL == bip || !SID_ok(bip, BIP_SID))
      return false;
   
   if (bip->rows < 1 || bip->rows > BIP_ROWS)
      return false;
   if (bip->cols < 1 || bip->cols > BIP_COLS)
      return false;

   for(int r = 0; r < bip->rows; r++)
   {
      for(int c = 0; c < bip->cols; c++)
         if (bip->a[r][c] < -MAX_COEF_VAL || bip->a[r][c] > MAX_COEF_VAL)
            return false;

      if (bip->b[r] < -MAX_COEF_VAL || bip->b[r] > MAX_COEF_VAL)
         return false;
   }
   return true;
}

void bip_free(BIP* bip)
{
   assert(bip_is_valid(bip));

   SID_del(bip);

   mem_check(bip);
   
   deallocate(bip);
}

BIP* bip_read(const char* filename)
{
   assert(NULL != filename);
   assert(0    <  strlen(filename));

#ifdef USE_SCANF
   char input[] = "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf";
#endif
   
   char      buf[MAX_LINE_LEN];
   char*     s;
   FILE*     fp;
   BIP*      bip   = allocate(1, sizeof(*bip));
   int       lines = 0;
   int       rows  = 0;
   int       cols;
   LINE_MODE mode  = READ_COLS;
   
   if (NULL == (fp = fopen(filename, "r")))
   {
      fprintf(stderr, "Can't open file %s\n", filename);
      goto read_error;
   }
   printf("Reading %s\n", filename);

   while(NULL != (s = fgets(buf, sizeof(buf), fp)))
   {
      char* t = strpbrk(s, "#\n\r"); 

      lines++;
      
      if (NULL != t) /* else line is not terminated or too long */
         *t = '\0';  /* clip comment or newline */
      
/* Skip over leading space
       */
      while(isspace(*s))
         s++;

      /* Skip over empty lines
       */
      if (!*s)  /* <=> (*s == '\0') */
         continue;

#ifdef USE_SCANF
      double v[33];
      int    nums;      
      
      nums = sscanf(s, input,
         &v[ 0], &v[ 1], &v[ 2], &v[ 3], &v[ 4], &v[ 5], &v[ 6], &v[ 7], 
         &v[ 8], &v[ 9], &v[10], &v[11], &v[12], &v[13], &v[14], &v[15],
         &v[16], &v[17], &v[18], &v[19], &v[20], &v[21], &v[22], &v[23],
         &v[24], &v[25], &v[26], &v[27], &v[28], &v[29], &v[30], &v[31]);

      if (nums < 1)
      {
         fprintf(stderr, "Error line %d: Expected number, got \"%s\"\n", lines, s);
         goto read_error;
      }

      /* do processing here:
       * mode tells which kind of line is next.
       */
      switch(mode)
      {
      case READ_COLS :
         bip->cols = (int)v[0];
         
         if (bip->cols < 1 || bip->cols > BIP_COLS)
         {
fprintf(stderr, "Error line %d: 1 <= cols=%d <= %d\n", lines, bip->cols, BIP_COLS);
goto read_error;
         }
         /* Insert <= into the input string
          */
         strncpy(input + bip->cols * 4, "<= ", 3);
         
         mode = READ_ROWS;
         break;
      case READ_ROWS :
         bip->rows = (int)v[0];

         if (bip->rows < 1 || bip->rows > BIP_ROWS)
         {
fprintf(stderr, "Error line %d: 1 <= rows <= %d\n", lines, BIP_ROWS);
            goto read_error;
         }
         mode = READ_COEF;
         break;
      case READ_COEF :
         if (rows >= bip->rows)
         {
            fprintf(stderr, "Error line=%d too many rows (expetced %d)\n",
               lines, bip->rows);
            goto read_error;
         }
         if (nums != bip->cols + 1)
         {
            fprintf(stderr, "Error line=%d Expected %d cols, got %d\n",
               lines, bip->cols, nums);
            goto read_error;
         }
         for(cols = 0; cols < bip->cols; cols++)
            bip->a[rows][cols] = v[cols];

         bip->b[rows] = v[cols];
         
         rows++;
         break;
      default :
         abort();
      }
#else         
      /* do processing here:
       * mode tells which kind of line is next.
       */
      switch(mode)
      {
      case READ_COLS :
         bip->cols = atoi(s);
         
         if (bip->cols < 1 || bip->cols > BIP_COLS)
         {
            fprintf(stderr, "Error line %d: 1 <= cols=%d <= %d\n", lines, bip->cols, BIP_COLS);
            goto read_error;
         }
         mode = READ_ROWS;
         break;
      case READ_ROWS :
         bip->rows = atoi(s);

         if (bip->rows < 1 || bip->rows > BIP_ROWS)
         {
            fprintf(stderr, "Error line %d: 1 <= rows <= %d\n", lines, BIP_ROWS);
            goto read_error;
         }
         mode = READ_COEF;
         break;
      case READ_COEF :
         if (rows >= bip->rows)
         {
            fprintf(stderr, "Error line=%d too many rows (expetced %d)\n",
               lines, bip->rows);
            goto read_error;
         }
#ifdef USE_STRTOK  // not  thread safe, rarely a problem in reading
         const char* delim = " \t";
         cols = 0;
         
         for(s = strtok(s, delim); s != NULL; s = strtok(NULL, delim)) /* not thread safe */
         {
            if (cols == bip->cols && strcmp(s, "<="))
               bip->b[rows] = atof(s);
            else if (cols < bip->cols)
            {
               bip->a[rows][cols] = atof(s);
               cols++;
            }
         }
#else
         for(cols = 0; cols < bip->cols; cols++)
         {
            bip->a[rows][cols] = (double)strtol(s, &t, 10);

            if (s == t)
            {
               fprintf(stderr, "Error line=%d Expected %d cols, got %d\n",
                  lines, bip->cols, cols);
               goto read_error;
            }
            s = t;
         }
         while(isspace(*s))
            s++;

         if (strncmp(s, "<=", 2))
         {
            fprintf(stderr, "Error line=%d Expected <= , got \"%s\"\n",
               lines, s);
            goto read_error;
         }
         s += 2;
         bip->b[rows] = (double)strtol(s, &t, 10);

         if (s == t)
         {
            fprintf(stderr, "Error line=%d Expected rhs\n", lines);
            goto read_error;
         }
#endif // !USE_STRTOK
         rows++;
         break;
      default :
         abort();
      }
#endif // !USE_SCANF
   }
   fclose(fp);

   if (rows != bip->rows)
   {
      fprintf(stderr, "Error line=%d Expetced %d rows, got %d\n",
         lines, bip->rows, rows);
      goto read_error;
   }

   SID_set(bip, BIP_SID);

   assert(bip_is_valid(bip));

   printf("Read %d rows, %d cols\n", bip->rows, bip->cols);
   
   return bip;

 read_error:
      deallocate(bip);
      return NULL;   
}

void bip_print(FILE* fp, const BIP* bip)
{
   assert(NULL != fp);
   assert(bip_is_valid(bip));

   for(int r = 0; r < bip->rows; r++)
   {
      for(int c = 0; c < bip->cols; c++)
         fprintf(fp, "%5g ", bip->a[r][c]);
      fprintf(fp, "<= %5g\n", bip->b[r]);
   }
}

void print_solu(
   FILE* fp,         
   int   vars,       
   const double* x)  
{
   assert(NULL != fp);
   assert(vars >= 0);
   assert(vars <  BIP_COLS);
   assert(NULL != x);

   for(int r = 0; r < vars; r++)
      fprintf(fp, "%2g ", x[r]);
   fputc('\n', fp);
}

bool bip_is_feasible(const BIP* bip, const double* x)
{
   assert(bip_is_valid(bip));
   
   int r;
   
   for(r = 0; r < bip->rows; r++)
   {
      double lhs = 0.0;
      
      for(int c = 0; c < bip->cols; c++)
      {
         assert(x[c] == 0.0 || x[c] == 1.0);
         
         lhs += bip->a[r][c] * x[c];
      }
      if (lhs > bip->b[r])
         break;
   }
   return r == bip->rows;
}

int bip_enumerate(const BIP* bip, bool with_output)
{
   assert(bip_is_valid(bip));
   
   clock_t start = clock();
   int     solus = 0;
   int     count = 0;
   
   /* generate all permutation
    */
   for(unsigned long long bitvec = 0; bitvec < (1uL << bip->cols); bitvec++)
   {      
      unsigned long long mask = 1;  // long long because long is 32 bit on 32 bit computers
      double             x[BIP_COLS];

      assert(sizeof(bitvec) * 8 > BIP_COLS); //lint !e506
      
      /* construct solution vector
       */
      for(int j = 0; j < bip->cols; j++, mask <<= 1)
         x[j] = (bitvec & mask) ? 1.0 : 0.0; 

      if (bip_is_feasible(bip, &x[0])) //lint !e772 
      {
         if (with_output)
            print_solu(stdout, bip->rows, &x[0]);
         
         solus++;
      }
      count++;
   }
   assert((unsigned)count == (1u << bip->cols));
   
   double elapsed = GET_SEC(start, clock());
   
   printf("Checked %d vectors in %.3f s = %.3f kvecs/s\n",
      count, elapsed, count / elapsed / 1000.0);

   return solus;
}

int main(int argc, const char** argv)
{
   if (argc < 2)
   {
      fprintf(stderr, "usage: %s filename", argv[0]);
      return EXIT_FAILURE;
   }
   BIP* bip = bip_read(argv[1]);

   if (NULL == bip)
      return EXIT_FAILURE;

   bip_print(stdout, bip);
   
   printf("%d solutions\n", bip_enumerate(bip, true));

   bip_free(bip);

   mem_maximum(stdout);
   
   return EXIT_SUCCESS;
}