1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using HeuristicLab.Common;


23  using HeuristicLab.Core;


24  using HeuristicLab.Data;


25  using HeuristicLab.Encodings.PermutationEncoding;


26  using HeuristicLab.Operators;


27  using HeuristicLab.Parameters;


28  using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;


29 


30  namespace HeuristicLab.Problems.LinearAssignment {


31  [Item("LinearAssignmentProblemSolver", "Uses the hungarian algorithm to solve linear assignment problems.")]


32  [StorableClass]


33  public sealed class LinearAssignmentProblemSolver : SingleSuccessorOperator {


34  private const int UNASSIGNED = 1;


35 


36  public ILookupParameter<DoubleMatrix> CostsParameter {


37  get { return (ILookupParameter<DoubleMatrix>)Parameters["Costs"]; }


38  }


39  public ILookupParameter<Permutation> AssignmentParameter {


40  get { return (ILookupParameter<Permutation>)Parameters["Assignment"]; }


41  }


42  public ILookupParameter<DoubleValue> QualityParameter {


43  get { return (ILookupParameter<DoubleValue>)Parameters["Quality"]; }


44  }


45 


46  [StorableConstructor]


47  private LinearAssignmentProblemSolver(bool deserializing) : base(deserializing) { }


48  private LinearAssignmentProblemSolver(LinearAssignmentProblemSolver original, Cloner cloner) : base(original, cloner) { }


49  public LinearAssignmentProblemSolver()


50  : base() {


51  Parameters.Add(new LookupParameter<DoubleMatrix>("Costs", LinearAssignmentProblem.CostsDescription));


52  Parameters.Add(new LookupParameter<Permutation>("Assignment", "The assignment solution to create."));


53  Parameters.Add(new LookupParameter<DoubleValue>("Quality", "The quality value of the solution."));


54  }


55 


56  public override IDeepCloneable Clone(Cloner cloner) {


57  return new LinearAssignmentProblemSolver(this, cloner);


58  }


59 


60  public override IOperation Apply() {


61  var costs = CostsParameter.ActualValue;


62  double quality;


63  var solution = Solve(costs, out quality);


64 


65  AssignmentParameter.ActualValue = new Permutation(PermutationTypes.Absolute, solution);


66  QualityParameter.ActualValue = new DoubleValue(quality);


67 


68  return base.Apply();


69  }


70 


71  /// <summary>


72  /// Uses the Hungarian algorithm to solve the linear assignment problem (LAP).


73  /// The LAP is defined as minimize f(p) = Sum(i = 1..N, c_{i, p(i)}) for a permutation p and an NxN cost matrix.


74  ///


75  /// The runtime complexity of the algorithm is O(n^3). The algorithm is deterministic and terminates


76  /// returning one of the optimal solutions and the corresponding quality.


77  /// </summary>


78  /// <remarks>


79  /// The algorithm is written similar to the fortran implementation given in http://www.seas.upenn.edu/qaplib/code.d/qapglb.f


80  /// </remarks>


81  /// <param name="costs">An NxN costs matrix.</param>


82  /// <param name="quality">The quality value of the optimal solution.</param>


83  /// <returns>The optimal solution.</returns>


84  public static int[] Solve(DoubleMatrix costs, out double quality) {


85  int length = costs.Rows;


86  // solve the linear assignment problem f(p) = Sum(i = 1..p, c_{i, p(i)})


87 


88  int[] rowAssign = new int[length], colAssign = new int[length];


89  double[] dualCol = new double[length], dualRow = new double[length];


90  for (int i = 0; i < length; i++) { // mark all positions as untouched


91  rowAssign[i] = UNASSIGNED;


92  colAssign[i] = UNASSIGNED;


93  }


94 


95  for (int i = 0; i < length; i++) { // find the minimum (base) level for each row


96  double min = costs[i, 0];


97  int minCol = 0;


98  dualCol[0] = min;


99  for (int j = 1; j < length; j++) {


100  if (costs[i, j] <= min) {


101  min = costs[i, j];


102  minCol = j;


103  }


104  if (costs[i, j] > dualCol[j])


105  dualCol[j] = costs[i, j];


106  }


107  dualRow[i] = min; // this will be the value of our dual variable


108  if (colAssign[minCol] == UNASSIGNED) {


109  colAssign[minCol] = i;


110  rowAssign[i] = minCol;


111  }


112  }


113 


114  for (int j = 0; j < length; j++) { // calculate the second dual variable


115  if (colAssign[j] != UNASSIGNED) dualCol[j] = 0;


116  else {


117  int minRow = 0;


118  for (int i = 0; i < length; i++) {


119  if (dualCol[j] > 0 && costs[i, j]  dualRow[i] < dualCol[j]) {


120  dualCol[j] = costs[i, j]  dualRow[i]; // the value is the original costs minus the first dual value


121  minRow = i;


122  }


123  }


124  if (rowAssign[minRow] == UNASSIGNED) {


125  colAssign[j] = minRow;


126  rowAssign[minRow] = j;


127  }


128  }


129  }


130 


131  // at this point costs_ij  dualRow_i  dualColumn_j results in a matrix that has at least one zero in every row and every column


132 


133  for (int i = 0; i < length; i++) { // try to make the remaining assignments


134  if (rowAssign[i] == UNASSIGNED) {


135  double min = dualRow[i];


136  for (int j = 0; j < length; j++) {


137  if (colAssign[j] == UNASSIGNED && (costs[i, j]  min  dualCol[j]).IsAlmost(0.0)) {


138  rowAssign[i] = j;


139  colAssign[j] = i;


140  break;


141  }


142  }


143  }


144  }


145 


146  bool[] marker = new bool[length];


147  double[] dplus = new double[length], dminus = new double[length];


148  int[] rowMarks = new int[length];


149 


150  for (int u = 0; u < length; u++) {


151  if (rowAssign[u] == UNASSIGNED) {


152  for (int i = 0; i < length; i++) {


153  rowMarks[i] = u;


154  marker[i] = false;


155  dplus[i] = double.MaxValue;


156  dminus[i] = costs[u, i]  dualRow[u]  dualCol[i];


157  }


158 


159  dplus[u] = 0;


160  int index = 1;


161  double minD = double.MaxValue;


162  while (true) {


163  minD = double.MaxValue;


164  for (int i = 0; i < length; i++) {


165  if (!marker[i] && dminus[i] < minD) {


166  minD = dminus[i];


167  index = i;


168  }


169  }


170 


171  if (colAssign[index] == UNASSIGNED) break;


172  marker[index] = true;


173  dplus[colAssign[index]] = minD;


174  for (int i = 0; i < length; i++) {


175  if (marker[i]) continue;


176  double compare = minD + costs[colAssign[index], i]  dualCol[i]  dualRow[colAssign[index]];


177  if (dminus[i] > compare) {


178  dminus[i] = compare;


179  rowMarks[i] = colAssign[index];


180  }


181  }


182 


183  } // while(true)


184 


185  while (true) {


186  colAssign[index] = rowMarks[index];


187  var ind = rowAssign[rowMarks[index]];


188  rowAssign[rowMarks[index]] = index;


189  if (rowMarks[index] == u) break;


190 


191  index = ind;


192  }


193 


194  for (int i = 0; i < length; i++) {


195  if (dplus[i] < double.MaxValue)


196  dualRow[i] += minD  dplus[i];


197  if (dminus[i] < minD)


198  dualCol[i] += dminus[i]  minD;


199  }


200  }


201  }


202 


203  quality = 0;


204  for (int i = 0; i < length; i++) {


205  quality += costs[i, rowAssign[i]];


206  }


207  return rowAssign;


208  }


209  }


210  }

