source: trunk/sources/HeuristicLab.Problems.LinearAssignment/3.3/LinearAssignmentProblemSolver.cs @ 11087

Last change on this file since 11087 was 11087, checked in by abeham, 5 years ago

#2131: Fixed maximization in LAP (by using negative cost matrix)

File size: 9.0 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2013 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
22using HeuristicLab.Common;
23using HeuristicLab.Core;
24using HeuristicLab.Data;
25using HeuristicLab.Encodings.PermutationEncoding;
26using HeuristicLab.Operators;
27using HeuristicLab.Parameters;
28using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
29
30namespace 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}
Note: See TracBrowser for help on using the repository browser.