bip_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 <limits.h>
#include <math.h>
#include <errno.h>
#include <float.h>
#include <fenv.h>

#include "mshell.h"  // or by allocate?
#include "splitline.h"
#include "bip_enum.h"

#if defined(USE_INT)
typedef int Numb;
#elif defined(USE_LONGLONG) 
typedef long long Numb;
#elif defined(USE_FLOAT) 
typedef float Numb;
#elif defined(USE_DOUBLE)
typedef double Numb;
#elif defined(USE_LONGDOUBLE)
typedef long double Numb;
#else
#error "No number type defined"
#endif

#define BIP_MAX_COLS 32
#define BIP_MAX_ROWS 32

typedef enum { READ_ERROR, READ_COLS, READ_ROWS, READ_COEF } LINE_MODE;

typedef enum { LE = 0, GE = 1, EQ = 2 } Sense;

struct binary_program
{
   int    rows;                   
   int    max_rows;               
   int    cols;                   
   Numb   a[BIP_MAX_ROWS][BIP_MAX_COLS];  
   Numb   b[BIP_MAX_ROWS];            
   Sense  sense[BIP_MAX_ROWS];        
};


#define MAX_LINE_LEN   512  

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

static Numb max_coef_val;
static Numb min_coef_val;
static bool int_coef;

static bool bip_is_valid(const BIP* bip)
{
   if (NULL == bip)
      return false;

   if (max_coef_val <=  min_coef_val) // bip_init not called.
      return false;
   
   if (bip->rows < 1 || bip->rows > BIP_MAX_ROWS)
      return false;
   if (bip->cols < 1 || bip->cols > BIP_MAX_COLS)
      return false;

   for(int r = 0; r < bip->rows; r++)
   {
      for(int c = 0; c < bip->cols; c++)
         if (bip->a[r][c] < min_coef_val || bip->a[r][c] > max_coef_val) // can only happen with floating point Numb
            return false;

      if (bip->b[r] < min_coef_val || bip->b[r] > max_coef_val) // can only happen with floating point Numb
         return false;
   }
   return true;
}

void bip_init()
{
#if defined(USE_INT)

   max_coef_val = INT_MAX;
   min_coef_val = INT_MIN;
   int_coef     = true;
   
#elif defined(USE_LONGLONG) 

   max_coef_val = LLONG_MAX;
   min_coef_val = LLONG_MIN;
   int_coef     = true;

#elif defined(USE_FLOAT) 

   max_coef_val =  powf(10, (float)FLT_DIG);
   min_coef_val = -powf(10, (float)FLT_DIG);
   int_coef     =  false;

#elif defined(USE_DOUBLE)

   max_coef_val =  pow(10, (double)DBL_DIG);
   min_coef_val = -pow(10, (double)DBL_DIG);
   int_coef     =  false;

#elif defined(USE_LONGDOUBLE)

   max_coef_val =  powl(10, (long double)LDBL_DIG);
   min_coef_val = -powl(10, (long double)LDBL_DIG);
   int_coef     =  false;

#else
#error "No number type defined"
#endif   

   //   printf("Min/Max coefficient values [%.0Lf,%.0Lf]\n",
   //      (long double)min_coef_val, (long double)max_coef_val);
}

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

   mem_check(bip);
   
   free(bip);
}

static bool bip_can_overflow(const BIP* bip)
{
   assert(bip_is_valid(bip));

   for(int r = 0; r < bip->rows; r++)
   {
      Numb row_max = 0;
      Numb row_min = 0;
      
      for(int c = 0; c < bip->cols; c++)
      {
         Numb val = bip->a[r][c];
         
         if (val > 0)
         {
            if (row_max < max_coef_val - val)
               row_max += val;
            else
               return fprintf(stderr, "Error: row %d numerical overflow\n", r), true;
         }
         else if (val < 0)
         {
            if (row_min > min_coef_val - val)
               row_min += val;
            else
               return fprintf(stderr, "Error: row %d numerical negative overflow\n", r), true;
         }
         else
         {
            assert(val == 0);
         }
      }
   }
   return false;
}


static LINE_MODE process_line(LINE_MODE mode, const LFS* lfs, const int lines, BIP* bip)
{
   /* If line is not empty
    */
   if (lfs_used_fields(lfs) > 0)
   {
      long double val;
      char*       t;
      
      /* do processing here:
       * mode tells which kind of line is next.
       */
      switch(mode)
      {
      case READ_COLS :
         if (lfs_used_fields(lfs) != 1)
            return fprintf(stderr, "Error line %d: Got %d fields, expected 1\n",
               lines, lfs_used_fields(lfs)),
               READ_ERROR;

         val = strtold(lfs_get_field(lfs, 0), &t);

         if (errno == ERANGE || isnan(val) || isinf(val) || rintl(val) != val || val < 1.0 || val > (long double)BIP_MAX_COLS || *t != '\0')
            return fprintf(stderr, "Error line=%d Number of cols %s=%Lf out of range or not integer %s\n",
               lines, lfs_get_field(lfs, 0), val, errno ? strerror(errno) : ""),
               READ_ERROR;

         bip->cols = (int)val;
         mode      = READ_ROWS;
         break;
      case READ_ROWS :
         if (lfs_used_fields(lfs) != 1)
            return fprintf(stderr, "Error line %d: Got %d fields, expected 1\n",
               lines, lfs_used_fields(lfs)),
               READ_ERROR;

         val = strtold(lfs_get_field(lfs, 0), &t);

         if (errno == ERANGE || isnan(val) || isinf(val) || rintl(val) != val || val < 1.0 || val > (long double)BIP_MAX_ROWS || *t != '\0')
            return fprintf(stderr, "Error line=%d Number of rows %s=%Lf out of range or not integer %s\n",
               lines, lfs_get_field(lfs, 0), val, errno ? strerror(errno) : ""),
               READ_ERROR;

         bip->max_rows = (int)val;
         mode          = READ_COEF;
         break;
      case READ_COEF :
         if (bip->rows >= bip->max_rows)
            return fprintf(stderr, "Error line=%d Expetced %d rows, got more\n",
               lines, bip->max_rows),
               READ_ERROR;

         if (lfs_used_fields(lfs) != bip->cols + 2)
            return fprintf(stderr, "Error line %d: Got %d fields, expected %d\n",
               lines, lfs_used_fields(lfs), bip->cols + 2),
               READ_ERROR;

         for(int col = 0; col < bip->cols + 2; col++)
         {
            if (col != bip->cols)
            {
               val = strtold(lfs_get_field(lfs, col), &t);

               if (errno == ERANGE || isnan(val) || isinf(val)
                  || (int_coef && rintl(val) != val)
                  || val > (long double)max_coef_val
                  || val < (long double)min_coef_val
                  || *t != '\0')
               return fprintf(stderr, "Error line=%d Number %s=%Lf out of range %s: %s\n",
                  lines, lfs_get_field(lfs, col), val, int_coef ? "or not integer" : "", errno ? strerror(errno) : ""),
                  READ_ERROR;

               if (col < bip->cols)
                  bip->a[bip->rows][col] = (Numb)val;
               else
                  bip->b[bip->rows] = (Numb)val;
            }
            else
            {
               const char* sense = lfs_get_field(lfs, bip->cols);
         
               if (!strcmp(sense, "<="))
                  bip->sense[bip->rows] = LE;
               else if (!strcmp(sense, ">="))
                  bip->sense[bip->rows] = GE;
               else if (!strcmp(sense, "=") || !strcmp(sense, "=="))
                  bip->sense[bip->rows] = EQ;
               else
                  return fprintf(stderr, "Error line=%d Expected <=, >=, ==, got \"%s\"\n",
                     lines, sense),
                     READ_ERROR;
            }
         }
         bip->rows++;
         break;
      default :
         abort();
      }
   }
   return mode;
}


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

   char      buf[MAX_LINE_LEN];
   char*     s;
   FILE*     fp;
   BIP*      bip   = calloc(1, sizeof(*bip));
   LFS*      lfs   = NULL;
   int       lines = 0;
   LINE_MODE mode  = READ_COLS;
   
   if (NULL == (fp = fopen(filename, "r")))
   {
      fprintf(stderr, "Can't open file %s\n", filename);
      free(bip);
      return NULL;
   }
   printf("Reading %s\n", filename);

   while(mode != READ_ERROR && NULL != (s = fgets(buf, sizeof(buf), fp)))
   {
      lines++;
      
      lfs  = lfs_split_line(lfs, s, "#");
      mode = process_line(mode, lfs, lines, bip);
   }
   fclose(fp);

   if (NULL != lfs)
      lfs_free(lfs);
   
   if (READ_ERROR == mode)
   {
      free(bip);

      return NULL;
   }
   if (bip->cols == 0 || bip->max_rows == 0 || bip->rows < bip->max_rows)
   {
      fprintf(stderr, "Error: unexpected EOF\n");
      free(bip);
      return NULL;
   }
   
   assert(bip->rows == bip->max_rows);

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

   if (bip_can_overflow(bip))
   {
      bip_free(bip);

      return NULL;
   }   
   return bip;
}

void bip_print(const BIP* bip, FILE* fp)
{
   const char* sense[3] = { "<=", ">=", "==" };
      
   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, "%Lf ", (long double)bip->a[r][c]);
      fprintf(fp, "%s %Lf\n", sense[bip->sense[r]], (long double)bip->b[r]);
   }
}

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

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

static 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];

         assert(fetestexcept(FE_ALL_EXCEPT & ~FE_INEXACT) == 0);
      }
      
      switch(bip->sense[r])
      {
      case LE :
         if (lhs <= bip->b[r])
            continue;
         break;
      case GE :
         if (lhs >= bip->b[r])
            continue;
         break;
      case EQ :
         if (lhs == bip->b[r])
            continue;
         break;
      default :
         abort();
      }
      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_MAX_COLS]; // having this always double is questionable

      assert(sizeof(bitvec) * 8 > BIP_MAX_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->cols, &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;
}