Free cookie consent management tool by TermsFeed Policy Generator

source: branches/ALPS/HeuristicLab.Problems.LinearAssignment/3.3/LinearAssignmentProblemSolver.cs @ 12035

Last change on this file since 12035 was 12018, checked in by pfleck, 10 years ago

#2269

  • merged trunk after 3.3.11 release
  • updated copyright and plugin version in ALPS plugin
  • removed old ALPS samples based on an userdefined alg
File size: 9.0 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 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.Optimization;
28using HeuristicLab.Parameters;
29using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
30
31namespace HeuristicLab.Problems.LinearAssignment {
32  [Item("LinearAssignmentProblemSolver", "Uses the hungarian algorithm to solve linear assignment problems.")]
33  [StorableClass]
34  public sealed class LinearAssignmentProblemSolver : SingleSuccessorOperator, ISingleObjectiveOperator {
35    private const int UNASSIGNED = -1;
36
37    public IValueLookupParameter<BoolValue> MaximizationParameter {
38      get { return (IValueLookupParameter<BoolValue>)Parameters["Maximization"]; }
39    }
40    public ILookupParameter<DoubleMatrix> CostsParameter {
41      get { return (ILookupParameter<DoubleMatrix>)Parameters["Costs"]; }
42    }
43    public ILookupParameter<Permutation> AssignmentParameter {
44      get { return (ILookupParameter<Permutation>)Parameters["Assignment"]; }
45    }
46    public ILookupParameter<DoubleValue> QualityParameter {
47      get { return (ILookupParameter<DoubleValue>)Parameters["Quality"]; }
48    }
49
50    [StorableConstructor]
51    private LinearAssignmentProblemSolver(bool deserializing) : base(deserializing) { }
52    private LinearAssignmentProblemSolver(LinearAssignmentProblemSolver original, Cloner cloner) : base(original, cloner) { }
53    public LinearAssignmentProblemSolver()
54      : base() {
55      Parameters.Add(new ValueLookupParameter<BoolValue>("Maximization", "Whether the costs should be maximized or minimized."));
56      Parameters.Add(new LookupParameter<DoubleMatrix>("Costs", LinearAssignmentProblem.CostsDescription));
57      Parameters.Add(new LookupParameter<Permutation>("Assignment", "The assignment solution to create."));
58      Parameters.Add(new LookupParameter<DoubleValue>("Quality", "The quality value of the solution."));
59    }
60
61    public override IDeepCloneable Clone(Cloner cloner) {
62      return new LinearAssignmentProblemSolver(this, cloner);
63    }
64
65    [StorableHook(HookType.AfterDeserialization)]
66    private void AfterDeserialization() {
67      // BackwardsCompatibility3.3
68      #region Backwards compatible code, remove with 3.4
69      if (!Parameters.ContainsKey("Maximization"))
70        Parameters.Add(new ValueLookupParameter<BoolValue>("Maximization", "Whether the costs should be maximized or minimized."));
71      #endregion
72    }
73
74    public override IOperation Apply() {
75      var costs = CostsParameter.ActualValue;
76      var maximization = MaximizationParameter.ActualValue.Value;
77      if (maximization) {
78        costs = (DoubleMatrix)costs.Clone();
79        for (int i = 0; i < costs.Rows; i++)
80          for (int j = 0; j < costs.Rows; j++)
81            costs[i, j] = -costs[i, j];
82      }
83      double quality;
84      var solution = Solve(costs, out quality);
85
86      AssignmentParameter.ActualValue = new Permutation(PermutationTypes.Absolute, solution);
87      if (maximization) quality = -quality;
88      QualityParameter.ActualValue = new DoubleValue(quality);
89
90      return base.Apply();
91    }
92
93    /// <summary>
94    /// Uses the Hungarian algorithm to solve the linear assignment problem (LAP).
95    /// 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.
96    ///
97    /// The runtime complexity of the algorithm is O(n^3). The algorithm is deterministic and terminates
98    /// returning one of the optimal solutions and the corresponding quality.
99    /// </summary>
100    /// <remarks>
101    /// The algorithm is written similar to the fortran implementation given in http://www.seas.upenn.edu/qaplib/code.d/qapglb.f
102    /// </remarks>
103    /// <param name="costs">An NxN costs matrix.</param>
104    /// <param name="quality">The quality value of the optimal solution.</param>
105    /// <returns>The optimal solution.</returns>
106    public static int[] Solve(DoubleMatrix costs, out double quality) {
107      int length = costs.Rows;
108      // solve the linear assignment problem f(p) = Sum(i = 1..|p|, c_{i, p(i)})
109
110      int[] rowAssign = new int[length], colAssign = new int[length];
111      double[] dualCol = new double[length], dualRow = new double[length];
112      for (int i = 0; i < length; i++) { // mark all positions as untouched
113        rowAssign[i] = UNASSIGNED;
114        colAssign[i] = UNASSIGNED;
115      }
116
117      for (int i = 0; i < length; i++) { // find the minimum (base) level for each row
118        double min = costs[i, 0];
119        int minCol = 0;
120        dualCol[0] = min;
121        for (int j = 1; j < length; j++) {
122          if (costs[i, j] <= min) {
123            min = costs[i, j];
124            minCol = j;
125          }
126          if (costs[i, j] > dualCol[j])
127            dualCol[j] = costs[i, j];
128        }
129        dualRow[i] = min; // this will be the value of our dual variable
130        if (colAssign[minCol] == UNASSIGNED) {
131          colAssign[minCol] = i;
132          rowAssign[i] = minCol;
133        }
134      }
135
136      for (int j = 0; j < length; j++) { // calculate the second dual variable
137        if (colAssign[j] != UNASSIGNED) dualCol[j] = 0;
138        else {
139          int minRow = 0;
140          for (int i = 0; i < length; i++) {
141            if (dualCol[j] > 0 && costs[i, j] - dualRow[i] < dualCol[j]) {
142              dualCol[j] = costs[i, j] - dualRow[i]; // the value is the original costs minus the first dual value
143              minRow = i;
144            }
145          }
146          if (rowAssign[minRow] == UNASSIGNED) {
147            colAssign[j] = minRow;
148            rowAssign[minRow] = j;
149          }
150        }
151      }
152
153      // 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
154
155      for (int i = 0; i < length; i++) { // try to make the remaining assignments
156        if (rowAssign[i] == UNASSIGNED) {
157          double min = dualRow[i];
158          for (int j = 0; j < length; j++) {
159            if (colAssign[j] == UNASSIGNED && (costs[i, j] - min - dualCol[j]).IsAlmost(0.0)) {
160              rowAssign[i] = j;
161              colAssign[j] = i;
162              break;
163            }
164          }
165        }
166      }
167
168      bool[] marker = new bool[length];
169      double[] dplus = new double[length], dminus = new double[length];
170      int[] rowMarks = new int[length];
171
172      for (int u = 0; u < length; u++) {
173        if (rowAssign[u] == UNASSIGNED) {
174          for (int i = 0; i < length; i++) {
175            rowMarks[i] = u;
176            marker[i] = false;
177            dplus[i] = double.MaxValue;
178            dminus[i] = costs[u, i] - dualRow[u] - dualCol[i];
179          }
180
181          dplus[u] = 0;
182          int index = -1;
183          double minD = double.MaxValue;
184          while (true) {
185            minD = double.MaxValue;
186            for (int i = 0; i < length; i++) {
187              if (!marker[i] && dminus[i] < minD) {
188                minD = dminus[i];
189                index = i;
190              }
191            }
192
193            if (colAssign[index] == UNASSIGNED) break;
194            marker[index] = true;
195            dplus[colAssign[index]] = minD;
196            for (int i = 0; i < length; i++) {
197              if (marker[i]) continue;
198              double compare = minD + costs[colAssign[index], i] - dualCol[i] - dualRow[colAssign[index]];
199              if (dminus[i] > compare) {
200                dminus[i] = compare;
201                rowMarks[i] = colAssign[index];
202              }
203            }
204
205          } // while(true)
206
207          while (true) {
208            colAssign[index] = rowMarks[index];
209            var ind = rowAssign[rowMarks[index]];
210            rowAssign[rowMarks[index]] = index;
211            if (rowMarks[index] == u) break;
212
213            index = ind;
214          }
215
216          for (int i = 0; i < length; i++) {
217            if (dplus[i] < double.MaxValue)
218              dualRow[i] += minD - dplus[i];
219            if (dminus[i] < minD)
220              dualCol[i] += dminus[i] - minD;
221          }
222        }
223      }
224
225      quality = 0;
226      for (int i = 0; i < length; i++) {
227        quality += costs[i, rowAssign[i]];
228      }
229      return rowAssign;
230    }
231  }
232}
Note: See TracBrowser for help on using the repository browser.