1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022013 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 IValueLookupParameter<BoolValue> MaximizationParameter {


37  get { return (IValueLookupParameter<BoolValue>)Parameters["Maximization"]; }


38  }


39  public ILookupParameter<DoubleMatrix> CostsParameter {


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


41  }


42  public ILookupParameter<Permutation> AssignmentParameter {


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


44  }


45  public ILookupParameter<DoubleValue> QualityParameter {


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


47  }


48 


49  [StorableConstructor]


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


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


52  public LinearAssignmentProblemSolver()


53  : base() {


54  Parameters.Add(new ValueLookupParameter<BoolValue>("Maximization", "Whether the costs should be maximized or minimized."));


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


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


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


58  }


59 


60  public override IDeepCloneable Clone(Cloner cloner) {


61  return new LinearAssignmentProblemSolver(this, cloner);


62  }


63 


64  [StorableHook(HookType.AfterDeserialization)]


65  private void AfterDeserialization() {


66  // BackwardsCompatibility3.3


67  #region Backwards compatible code, remove with 3.4


68  if (!Parameters.ContainsKey("Maximization"))


69  Parameters.Add(new ValueLookupParameter<BoolValue>("Maximization", "Whether the costs should be maximized or minimized."));


70  #endregion


71  }


72 


73  public override IOperation Apply() {


74  var costs = CostsParameter.ActualValue;


75  var maximization = MaximizationParameter.ActualValue.Value;


76  if (maximization) {


77  costs = (DoubleMatrix)costs.Clone();


78  for (int i = 0; i < costs.Rows; i++)


79  for (int j = 0; j < costs.Rows; j++)


80  costs[i, j] = costs[i, j];


81  }


82  double quality;


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


84 


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


86  if (maximization) quality = quality;


87  QualityParameter.ActualValue = new DoubleValue(quality);


88 


89  return base.Apply();


90  }


91 


92  /// <summary>


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


94  /// 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.


95  ///


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


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


98  /// </summary>


99  /// <remarks>


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


101  /// </remarks>


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


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


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


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


106  int length = costs.Rows;


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


108 


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


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


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


112  rowAssign[i] = UNASSIGNED;


113  colAssign[i] = UNASSIGNED;


114  }


115 


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


117  double min = costs[i, 0];


118  int minCol = 0;


119  dualCol[0] = min;


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


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


122  min = costs[i, j];


123  minCol = j;


124  }


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


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


127  }


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


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


130  colAssign[minCol] = i;


131  rowAssign[i] = minCol;


132  }


133  }


134 


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


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


137  else {


138  int minRow = 0;


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


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


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


142  minRow = i;


143  }


144  }


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


146  colAssign[j] = minRow;


147  rowAssign[minRow] = j;


148  }


149  }


150  }


151 


152  // 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


153 


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


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


156  double min = dualRow[i];


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


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


159  rowAssign[i] = j;


160  colAssign[j] = i;


161  break;


162  }


163  }


164  }


165  }


166 


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


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


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


170 


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


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


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


174  rowMarks[i] = u;


175  marker[i] = false;


176  dplus[i] = double.MaxValue;


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


178  }


179 


180  dplus[u] = 0;


181  int index = 1;


182  double minD = double.MaxValue;


183  while (true) {


184  minD = double.MaxValue;


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


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


187  minD = dminus[i];


188  index = i;


189  }


190  }


191 


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


193  marker[index] = true;


194  dplus[colAssign[index]] = minD;


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


196  if (marker[i]) continue;


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


198  if (dminus[i] > compare) {


199  dminus[i] = compare;


200  rowMarks[i] = colAssign[index];


201  }


202  }


203 


204  } // while(true)


205 


206  while (true) {


207  colAssign[index] = rowMarks[index];


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


209  rowAssign[rowMarks[index]] = index;


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


211 


212  index = ind;


213  }


214 


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


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


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


218  if (dminus[i] < minD)


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


220  }


221  }


222  }


223 


224  quality = 0;


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


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


227  }


228  return rowAssign;


229  }


230  }


231  }

