[7873] | 1 | #region License Information
|
---|
| 2 | /* HeuristicLab
|
---|
[16565] | 3 | * Copyright (C) 2002-2019 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
|
---|
[7873] | 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;
|
---|
[8022] | 23 | using HeuristicLab.Core;
|
---|
[7873] | 24 | using HeuristicLab.Data;
|
---|
[8022] | 25 | using HeuristicLab.Encodings.PermutationEncoding;
|
---|
| 26 | using HeuristicLab.Operators;
|
---|
[11970] | 27 | using HeuristicLab.Optimization;
|
---|
[8022] | 28 | using HeuristicLab.Parameters;
|
---|
[16565] | 29 | using HEAL.Attic;
|
---|
[7873] | 30 |
|
---|
| 31 | namespace HeuristicLab.Problems.LinearAssignment {
|
---|
[8022] | 32 | [Item("LinearAssignmentProblemSolver", "Uses the hungarian algorithm to solve linear assignment problems.")]
|
---|
[16565] | 33 | [StorableType("ABF202CC-44E4-4208-9EBF-1A104806358F")]
|
---|
[11970] | 34 | public sealed class LinearAssignmentProblemSolver : SingleSuccessorOperator, ISingleObjectiveOperator {
|
---|
[7873] | 35 | private const int UNASSIGNED = -1;
|
---|
| 36 |
|
---|
[11087] | 37 | public IValueLookupParameter<BoolValue> MaximizationParameter {
|
---|
| 38 | get { return (IValueLookupParameter<BoolValue>)Parameters["Maximization"]; }
|
---|
| 39 | }
|
---|
[8022] | 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]
|
---|
[16565] | 51 | private LinearAssignmentProblemSolver(StorableConstructorFlag _) : base(_) { }
|
---|
[8022] | 52 | private LinearAssignmentProblemSolver(LinearAssignmentProblemSolver original, Cloner cloner) : base(original, cloner) { }
|
---|
| 53 | public LinearAssignmentProblemSolver()
|
---|
| 54 | : base() {
|
---|
[11087] | 55 | Parameters.Add(new ValueLookupParameter<BoolValue>("Maximization", "Whether the costs should be maximized or minimized."));
|
---|
[8022] | 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 |
|
---|
[11087] | 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 |
|
---|
[8022] | 74 | public override IOperation Apply() {
|
---|
| 75 | var costs = CostsParameter.ActualValue;
|
---|
[11087] | 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 | }
|
---|
[8022] | 83 | double quality;
|
---|
| 84 | var solution = Solve(costs, out quality);
|
---|
| 85 |
|
---|
| 86 | AssignmentParameter.ActualValue = new Permutation(PermutationTypes.Absolute, solution);
|
---|
[11087] | 87 | if (maximization) quality = -quality;
|
---|
[8022] | 88 | QualityParameter.ActualValue = new DoubleValue(quality);
|
---|
| 89 |
|
---|
| 90 | return base.Apply();
|
---|
| 91 | }
|
---|
| 92 |
|
---|
[7873] | 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 | }
|
---|