source: stable/HeuristicLab.Problems.LinearAssignment/3.3/LinearAssignmentProblemSolver.cs @ 10032

Last change on this file since 10032 was 9456, checked in by swagner, 9 years ago

Updated copyright year and added some missing license headers (#1889)

File size: 8.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 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}
Note: See TracBrowser for help on using the repository browser.