C CODE

#ifndef HUNGARIAN_H
#define HUNGARIAN_H

#ifdef __cplusplus
extern "C" {
#endif

#include <sys/types.h> // for size_t
#include <string.h>
#include <stdio.h>
#include <assert.h>
#include <stdlib.h>
#define USAGE "Usage: ./test [-m <m>] [-n <n>] [-f <fname>]"

#include <math.h>
#include <limits.h> 
#include <sys/param.h>

#ifndef bzero
  #define bzero(b,len) (memset((b), '\0', (len)), (void) 0)
#endif

/* are we maximizing or minimizing? */
#define HUNGARIAN_MIN 0
#define HUNGARIAN_MAX 1


typedef struct
{
  int* i;
  int* j;
  int k;
} hungarian_sequence_t;


typedef struct
{
  size_t m,n;  // problem dimensions
  int** r;     // the rating (utility) matrix
  int** q;     // the Q matrix
  int* u;      // the U vector
  int* v;      // the V vector
  int* ess_rows; // list of essential rows
  int* ess_cols; // list of essential columns
  hungarian_sequence_t seq; // sequence of i's and j's
  int row_total, col_total; // row and column totals
  int* a;  // assignment vector
  int maxutil;  // maximum utility
  int mode; // are we maximizing or minimizing?
} hungarian_t;

/*
 * initialize the given object as an mXn problem.  allocates storage, which
 * should be freed with hungarian_fini().
 */
void hungarian_init(hungarian_t* prob, int* r, size_t m, size_t n, int mode);

/*
 * frees storage associated with the given problem object.  you must have
 * called hungarian_init() first.
 */
void hungarian_fini(hungarian_t* prob);

/*
 * solve the given problem.  runs the Hungarian Method on the rating matrix
 * to produce optimal assignment, which is stored in the vector prob->a.
 * you must have called hungarian_init() first.
 */
void hungarian_solve(hungarian_t* prob);

/*
 * prints out the resultant assignment in a 0-1 matrix form.  also computes
 * and prints out the benefit from the assignment.  you must have called
 * hungarian_solve() first.
 */
void hungarian_print_assignment(hungarian_t* prob);

/*
 * prints out the rating matrix for the given problem.  you must have called
 * hungarian_solve() first.
 */
void hungarian_print_rating(hungarian_t* prob);

/*
 * check whether an assigment is feasible.  returns 1 if the assigment is
 * feasible, 0 otherwise.  you must have called hungarian_solve() first.
 */
int hungarian_check_feasibility(hungarian_t* prob);

/*
 * computes and returns the benefit from the assignment.  you must have
 * called hungarian_solve() first.
 */
int hungarian_benefit(hungarian_t* prob);

/*
 * makes and returns a pointer to an mXn rating matrix with values uniformly 
 * distributed between 1 and MAXUTIL
 */
int* make_random_r(size_t m, size_t n);
int* make_r_from_ORlib(char* fname, int* m, int* n);

#ifdef __cplusplus
}
#endif

#endif

typedef enum
{
  HUNGARIAN_ERROR,
  HUNGARIAN_ROUTINE_ONE,
  HUNGARIAN_ROUTINE_TWO,
  HUNGARIAN_DONE
} hungarian_code_t;

// the Q matrix is filled with instances of this type
typedef enum
{
  HUNGARIAN_ZERO,
  HUNGARIAN_ONE,
  HUNGARIAN_STAR
} hungarian_q_t;


// some internal routines
void hungarian_make_cover(hungarian_t* prob);
void hungarian_build_q(hungarian_t* prob);
void hungarian_add_stars(hungarian_t* prob);
hungarian_code_t hungarian_routine_one(hungarian_t* prob);
hungarian_code_t hungarian_routine_two(hungarian_t* prob);

/*
 * initialize the given object as an mXn problem.  allocates storage, which
 * should be freed with hungarian_fini()
 */
void hungarian_init(hungarian_t* prob, int* r, size_t m, size_t n, int mode)
{
  int i,j;
  // did we get a good object pointer?
  assert(prob);

  // we can't work with matrices that have more rows than columns.
  //
  // TODO: automatically transpose such matrices
  assert(m<=n);

  // init the struct
  prob->m = m;
  prob->n = n;
  assert(prob->r = (int**)calloc(m,sizeof(int*)));
  assert(prob->q = (int**)calloc(m,sizeof(int*)));
  prob->mode = mode;
  prob->maxutil = 0;
  for(i=0;i<m;i++)
  {
    assert(prob->r[i] = (int*)calloc(n,sizeof(int)));
    assert(prob->q[i] = (int*)calloc(n,sizeof(int)));
    for(j=0;j<n;j++)
    {
      prob->r[i][j] = r[i*n+j];
      if(prob->r[i][j] > prob->maxutil)
        prob->maxutil = prob->r[i][j];
    }
  }
  // if we're going to minimize, rather than maximize, we need to subtract
  // each utility from the maximum.  this operation will be reversed before
  // computing the benefit.
  if(mode == HUNGARIAN_MIN)
  {
    for(i=0;i<m;i++)
    {
      for(j=0;j<n;j++)
        prob->r[i][j] = prob->maxutil - prob->r[i][j];
    }
  }

  assert(prob->a = (int*)calloc(m,sizeof(int)));
  assert(prob->u = (int*)calloc(m,sizeof(int)));
  assert(prob->v = (int*)calloc(n,sizeof(int)));
  assert(prob->seq.i = (int*)calloc(m*n,sizeof(int)));
  assert(prob->seq.j = (int*)calloc(m*n,sizeof(int)));
  assert(prob->ess_rows = (int*)calloc(m,sizeof(int)));
  assert(prob->ess_cols = (int*)calloc(n,sizeof(int)));
}

/*
 * frees storage associated with the given problem object.  you must have
 * called hungarian_init() first.
 */
void hungarian_fini(hungarian_t* prob)
{
  int i;
  assert(prob);
  for(i=0;i<prob->m;i++)
  {
    free(prob->q[i]);
    free(prob->r[i]);
  }
  free(prob->q);
  free(prob->r);
  free(prob->a);
  free(prob->u);
  free(prob->v);
  free(prob->seq.i);
  free(prob->seq.j);
  free(prob->ess_rows);
  free(prob->ess_cols);
}

/*
 * make an initial cover
 */
void hungarian_make_cover(hungarian_t* prob)
{
  int i,j;
  int* row_max; 
  int* col_max; 

  assert(prob);
  assert(row_max = (int*)calloc(prob->m,sizeof(int)));
  assert(col_max = (int*)calloc(prob->n,sizeof(int)));

  prob->row_total = prob->col_total = 0;

  // find the max in each row (col) and sum them
  for(i=0;i<prob->m;i++)
  {
    for(j=0;j<prob->n;j++)
    {
      if(prob->r[i][j] > row_max[i])
        row_max[i] = prob->r[i][j];
    }
    prob->row_total += row_max[i];
  }

  for(j=0;j<prob->n;j++)
  {
    for(i=0;i<prob->m;i++)
    {
      if(prob->r[i][j] > col_max[j])
        col_max[j] = prob->r[i][j];
    }
    prob->col_total += col_max[j];
  }

  // make u and v into an initial cover, based on row and col totals
  if(prob->row_total <= prob->col_total)
    memcpy(prob->u,row_max,sizeof(int)*prob->m);
  else
    memcpy(prob->v,col_max,sizeof(int)*prob->n);

  free(row_max);
  free(col_max);
}

void hungarian_build_q(hungarian_t* prob)
{
  int i,j;
  assert(prob);

  for(i=0;i<prob->m;i++)
  {
    for(j=0;j<prob->n;j++)
    {
      if((prob->u[i] + prob->v[j]) == prob->r[i][j])
      {
        if(prob->q[i][j] == HUNGARIAN_ZERO)
          prob->q[i][j] = HUNGARIAN_ONE;
      }
      else
        prob->q[i][j] = HUNGARIAN_ZERO;
    }
  }
}

void hungarian_add_stars(hungarian_t* prob)
{
  int i,j,k;

  assert(prob);

  if(prob->row_total <= prob->col_total)
  {
    for(i=0;i<prob->m;i++)
    {
      for(j=0;j<prob->n;j++)
      {
        if(prob->q[i][j] == HUNGARIAN_ONE)
        {
          for(k=0;k<prob->m;k++)
          {
            if((k!=i) && (prob->q[k][j] == HUNGARIAN_STAR))
              break;
          }
          if(k==prob->m)
          {
            prob->q[i][j] = HUNGARIAN_STAR;
            break;
          }
        }
      }
    }
  }
  else
  {
    for(j=0;j<prob->n;j++)
    {
      for(i=0;i<prob->m;i++)
      {
        if(prob->q[i][j] == HUNGARIAN_ONE)
        {
          for(k=0;k<prob->n;k++)
          {
            if((k!=j) && (prob->q[i][k] == HUNGARIAN_STAR))
              break;
          }
          if(k==prob->n)
          {
            prob->q[i][j] = HUNGARIAN_STAR;
            break;
          }
        }
      }
    }
  }
}

/*
 * Kuhn's Routine I
 */
hungarian_code_t hungarian_routine_one(hungarian_t* prob)
{
  int i,j,k;
  char jumpcase1,jumpprelim;

  assert(prob);
  prob->seq.k=0;
  for(i=0;i<prob->m;i++)
    prob->ess_rows[i]=0;

  j=0;
  // look for a 1* in each column
  while(j<prob->n)
  {
    i=0;
    while(i<prob->m)
    {
      if(prob->q[i][j] == HUNGARIAN_STAR)
      {
        // found a 1*; next column
        break;
      }
      i++;
    }
    if(i<prob->m)
    {
      // found a 1*; next column
      j++;
    }
    else
    {
      // didn't find a 1*; column j is eligible; search it for a 1
      i=0;
      while(i<prob->m)
      {
        if(prob->q[i][j] == HUNGARIAN_ONE)
        {
          // found a 1 (i,j); start recording
          prob->seq.i[prob->seq.k] = i;
          prob->seq.j[prob->seq.k] = j;
          prob->seq.k++;
          // CASE 1
          jumpprelim=0;
          while(!jumpprelim)
          {
            // look in row ik for a 1*
            j=0;
            jumpcase1=0;
            while(j<prob->n && !jumpcase1)
            {
              if(prob->q[i][j] == HUNGARIAN_STAR)
              {
                // CASE 2
                i=0;
                while(!jumpcase1)
                {
                  // found a 1* in (ik,jk); search col jk for a 1
                  while(i<prob->m)
                  {
                    if(prob->q[i][j] == HUNGARIAN_ONE)
                    {
                      // test i for distinctness
                      for(k=0;k<prob->seq.k;k++)
                      {
                        if(prob->seq.i[k] == i)
                          break;
                      }
                      if(k<prob->seq.k)
                      {
                        // i is not distinct
                        i++;
                        continue;
                      }
                      else
                      {
                        // i is distinct; record and back to Case 1
                        prob->seq.i[prob->seq.k] = prob->seq.i[prob->seq.k-1];
                        prob->seq.j[prob->seq.k] = j;
                        prob->seq.k++;
                        prob->seq.i[prob->seq.k] = i;
                        prob->seq.j[prob->seq.k] = j;
                        prob->seq.k++;

                        jumpcase1=1;
                        break;
                      }
                    }
                    else
                      i++;
                  }
                  if(i>=prob->m)
                  {
                    // didn't find a 1 in col jk; row ik is essential
                    prob->ess_rows[prob->seq.i[prob->seq.k-1]] = 1;
                    
                    // delete last two elts of sequence
                    j=prob->seq.j[prob->seq.k-1];
                    i=prob->seq.i[prob->seq.k-1]+1;

                    if(prob->seq.k > 1)
                    {
                      // back to Case 2
                      prob->seq.k-=2;
                      continue;
                    }
                    else
                    {
                      // k==1; back to prelim search for 1 in (i1+1,j0)
                      prob->seq.k--;
                      jumpcase1=jumpprelim=1;
                      break;
                    }
                  }
                }
              }
              else
                j++;
            }
            if(j>=prob->n)
            {
              // didn't find a 1* in row ik; toggle and Alternative Ia
              for(;prob->seq.k > 0;prob->seq.k--)
              {
                if(prob->q[prob->seq.i[prob->seq.k-1]][prob->seq.j[prob->seq.k-1]]
                   == HUNGARIAN_ONE)
                {
                  prob->q[prob->seq.i[prob->seq.k-1]][prob->seq.j[prob->seq.k-1]] 
                          = HUNGARIAN_STAR;
                }
                else
                {
                  prob->q[prob->seq.i[prob->seq.k-1]][prob->seq.j[prob->seq.k-1]] 
                          = HUNGARIAN_ONE;
                }
              }
              // Alternative Ia
              return(HUNGARIAN_ROUTINE_ONE);
            }
          }
        }
        else
          i++;
      }
      // didn't find a 1 in col j; next col
      j++;
    }
  }
  // out of cols; Alternative Ib
  // determine ess cols
  for(j=0;j<prob->n;j++)
  {
    prob->ess_cols[j]=0;
    for(i=0;i<prob->m;i++)
    {
      if(prob->q[i][j] == HUNGARIAN_STAR && !prob->ess_rows[i])
      {
        prob->ess_cols[j] = 1;
        break;
      }
    }
  }
  return(HUNGARIAN_ROUTINE_TWO);
}


/*
 * Kuhn's Routine II
 */
hungarian_code_t hungarian_routine_two(hungarian_t* prob)
{
  int i,j,d,m,dtmp;
  int oldsum,newsum;

  assert(prob);

  oldsum = 0;
  for(i=0;i<prob->m;i++)
    oldsum += prob->u[i];
  for(j=0;j<prob->n;j++)
    oldsum += prob->v[j];

  d=0;
  for(i=0;i<prob->m;i++)
  {
    if(prob->ess_rows[i])
      continue;
    for(j=0;j<prob->n;j++)
    {
      if(prob->ess_cols[j])
        continue;
      dtmp = prob->u[i] + prob->v[j] - prob->r[i][j];
      if(dtmp<0)
      {
        printf("SUPERMOO: %d + %d < %d\n", 
               prob->u[i], prob->v[j], prob->r[i][j]);
      }
      if(!d || ((dtmp>0) && dtmp < d))
        d = dtmp;
    }
  }

  if(d<0)
    printf("MOO: %d < 0\n", d);

  //if(d<=0)
  if(!d)
    return(HUNGARIAN_DONE);
  else
  {
    // is there some u[i]==0?
    for(i=0;i<prob->m;i++)
    {
      if(!prob->u[i])
        break;
    }
    if(i < prob->m)
    {
      // CASE 2; some u[i] == 0; compute m as the min of d and inessential v[j]
      m = d;
      for(j=0;j<prob->n;j++)
      {
        if((!prob->ess_cols[j]) && (prob->v[j] < m))
          m = prob->v[j];
      }

      // adjust the cover
      for(i=0;i<prob->m;i++)
      {
        if(prob->ess_rows[i])
          prob->u[i] += m;
      }
      for(j=0;j<prob->n;j++)
      {
        if(!prob->ess_cols[j])
          prob->v[j] -= m;
      }
    }
    else
    {
      // CASE 1; all u[i] > 0; compute m as the min of d and inessential u[j]
      m = d;
      for(i=0;i<prob->m;i++)
      {
        if((!prob->ess_rows[i]) && (prob->u[i] < m))
          m = prob->u[i];
      }

      // adjust the cover
      for(i=0;i<prob->m;i++)
      {
        if(!prob->ess_rows[i])
          prob->u[i] -=m;
      }
      for(j=0;j<prob->n;j++)
      {
        if(prob->ess_cols[j])
          prob->v[j] += m;
      }
    }
  }

  newsum = 0;

  for(i=0;i<prob->m;i++)
  {
    for(j=0;j<prob->n;j++)
    {
      if(prob->u[i]+prob->v[j]<prob->r[i][j])
      {
        printf("SUPERMOO (%d,%d): %d + %d < %d\n", i,j,
               prob->u[i], prob->v[j], prob->r[i][j]);
      }
    }
  }

  // Alternative IIa; build new q and return to routine I
  hungarian_build_q(prob);
  return(HUNGARIAN_ROUTINE_ONE);
}


/*
 * solve the given problem.  runs the Hungarian Method on the rating matrix
 * to produce optimal assignment, which is stored in the vector prob->a.
 * you must have called hungarian_init() first.
 */
void 
hungarian_solve(hungarian_t* prob)
{
  hungarian_code_t state = HUNGARIAN_ROUTINE_ONE;
  int i,j;

  assert(prob);

  bzero(prob->u,sizeof(int)*prob->m);
  bzero(prob->v,sizeof(int)*prob->n);
  prob->seq.k=0;
  bzero(prob->ess_rows,sizeof(int)*prob->m);
  bzero(prob->ess_cols,sizeof(int)*prob->n);
  for(i=0;i<prob->m;i++)
    bzero(prob->q[i],sizeof(int)*prob->n);

  hungarian_make_cover(prob);
  hungarian_build_q(prob);
  hungarian_add_stars(prob);

  while(state != HUNGARIAN_DONE)
  {
    if(state == HUNGARIAN_ROUTINE_ONE)
    {
      state = hungarian_routine_one(prob);
    }
    else
    {
      state = hungarian_routine_two(prob);
    }
  }

  // fill in the assignment vector
  for(i=0;i<prob->m;i++)
  {
    for(j=0;j<prob->n;j++)
    {
      if(prob->q[i][j] == HUNGARIAN_STAR)
      {
        prob->a[i] = j;
        break;
      }
    }
  }
}

/*
 * prints out the current state of the Q matrix.  also computes and prints
 * out the benefit from the current assignment.  mostly useful for
 * debugging.
 */
void hungarian_print_stars(hungarian_t* prob)
{
  int i,j;
  int benefit=0;
  assert(prob);
  puts("\nQ: ");
  for(i=0;i<prob->m;i++)
  {
    printf("  [ ");
    for(j=0;j<prob->n;j++)
    {
      if(prob->q[i][j] == HUNGARIAN_ZERO)
        printf("%4d ", 0);
      else if(prob->q[i][j] == HUNGARIAN_ONE)
        printf("%4d ", 1);
      else
      {
        printf("%3d%s ", 1,"*");
        if(prob->mode == HUNGARIAN_MIN)
          benefit += prob->maxutil - prob->r[i][j];
        else
          benefit += prob->r[i][j];
      }
    }
    puts(" ]");
  }
  printf("\nBenefit: %d\n\n", benefit);
}

/*
 * prints out the resultant assignment in a 0-1 matrix form.  you must have 
 * called hungarian_solve() first.
 */
void 
hungarian_print_assignment(hungarian_t* prob)
{
  int i,j;
  assert(prob);
  puts("\nA:");
  for(i=0;i<prob->m;i++)
  {
    printf("  [ ");
    for(j=0;j<prob->n;j++)
      printf("%4d ", (j==prob->a[i]) ? 1 : 0);
    printf(" ]\n");
  }
}

/*
 * prints out the rating matrix for the given problem.  you must have called
 * hungarian_solve() first.
 */
void 
hungarian_print_rating(hungarian_t* prob)
{
  int i,j;

  puts("\nR: ");
  for(i=0;i<prob->m;i++)
  {
    printf("  [ ");
    for(j=0;j<prob->n;j++)
    {
      printf("%4d ", prob->r[i][j]);
    }
    puts(" ]");
  }
}

/*
 * check whether an assigment is feasible.  returns 1 if the assigment is
 * feasible, 0 otherwise.  you must have called hungarian_solve() first.
 */
int
hungarian_check_feasibility(hungarian_t* prob)
{
  int i,j;
  char assigned;

  // check for over/under assigned rows
  for(i=0;i<prob->m;i++)
  {
    assigned=0;
    for(j=0;j<prob->n;j++)
    {
      if(prob->q[i][j] == HUNGARIAN_STAR)
      {
        if(assigned)
          return(0);
        else
          assigned=1;
      }
    }
    if((prob->m <= prob->n) && !assigned)
      return(0);
  }
  // check for over/under assigned cols
  for(j=0;j<prob->n;j++)
  {
    assigned=0;
    for(i=0;i<prob->m;i++)
    {
      if(prob->q[i][j] == HUNGARIAN_STAR)
      {
        if(assigned)
          return(0);
        else
          assigned=1;
      }
    }
    if((prob->n <= prob->m) && !assigned)
      return(0);
  }
  return(1);
}

/*
 * computes and returns the benefit from the assignment.  you must have
 * called hungarian_solve() first.
 */
int 
hungarian_benefit(hungarian_t* prob)
{
  int i;
  int benefit=0;
  assert(prob);
  for(i=0;i<prob->m;i++)
  {
    if(prob->mode == HUNGARIAN_MIN)
      benefit += prob->maxutil - prob->r[i][prob->a[i]];
    else
      benefit += prob->r[i][prob->a[i]];
  }

  return(benefit);
}

#ifndef max
  #define max(a, b) ((a) > (b) ? (a) : (b))
#endif

// problem dimensions
int m,n;

// input data file
char input_fname[PATH_MAX];

void parse_args(int argc, char** argv);

int
main(int argc, char** argv)
{
  hungarian_t prob;

  /*
  int r[4][4] =  {{  363, 626, 376,  46  },
                  {  802, 993, 784, 953  },
                  {  673,  23, 823, 207  },
                  {  380, 451, 520, 646  }};
                  */
  int* r;
  
  m = n = 4;

  parse_args(argc,argv);

  if(!strlen(input_fname))
    r = make_random_r(m,n);
  else 
    r = make_r_from_ORlib(input_fname,&m,&n);

  if(!r)
  {
    puts("got NULL input");
    exit(-1);
  }

  hungarian_init(&prob,(int*)r,m,n,HUNGARIAN_MIN);
  hungarian_print_rating(&prob);
  hungarian_solve(&prob);
  hungarian_print_assignment(&prob);

  printf("\nfeasible? %s\n", 
         (hungarian_check_feasibility(&prob)) ? "yes" : "no");
  printf("benefit: %d\n\n", hungarian_benefit(&prob));

  hungarian_fini(&prob);
  free(r);
  return(0);
}

void
parse_args(int argc, char** argv)
{
  int i;
  // parse command-line args
  for(i=1; i<argc; i++)
  {
    if(!strcmp(argv[i],"-m"))
    {
      if(++i < argc)
      {
        m = atoi(argv[i]);
      }
      else
      {
        puts(USAGE);
        exit(-1);
      }
    }
    else if(!strcmp(argv[i],"-n"))
    {
      if(++i < argc)
      {
        n = atoi(argv[i]);
      }
      else
      {
        puts(USAGE);
        exit(-1);
      }
    }
    else if(!strcmp(argv[i],"-f"))
    {
      if(++i < argc)
      {
        strncpy(input_fname,argv[i],sizeof(input_fname));
      }
      else
      {
        puts(USAGE);
        exit(-1);
      }
    }
    else
    {
      puts(USAGE);
      exit(-1);
    }
  }
}