hungarian.c

匈牙利算法求解最优分配问题

                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: ");
  if (prob->mode == HUNGARIAN_MIN)
  {
	  for(i=0;i<prob->m;i++)
	  {
		  printf("  [ ");
		  for(j=0;j<prob->n;j++)
		  {
			  printf("%4d ", prob->maxutil - prob->r[i][j]);
		  }
		  puts(" ]");
	  }

  }
  else
  { 
	  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);
}