Linear assignment problems, part III
Here's the code. See the previous posting for background.
using System;
namespace NathanBrixius.AssignmentLib
{
/// <summary>Linear assignment problem solver.<summary>
/// <remarks>
/// This code is a port of Fortran code by R.E. Burkard and U. Derigs.
/// Assignment and matching problems: Solution methods with fortran programs,
/// volume 184 of Lecture Notes in Economics and Mathematical Systems.
/// Springer, Berlin, 1980. The original code is available on the QAPLIB site.
/// </remarks>
public class LinearAssignment
{
public static void Main(params string[] args)
{
double[][] L = new double[];
for (int i = 0; i < L.Length; i++)
{
L = new double;
for (int j = 0; j < L.Length; j++)
{
L = i + j;
}
}
int[] p = new int[L.Length];
int[] ip = new int[L.Length];
double[] ys = new double[L.Length];
double[] yt = new double[L.Length];
Solve(L, p, ip, ys, yt);
}
/// <summary>Evaluate a permutation against the given LAP.<summary>
/// <param name="L">LAP cost matrix.</param>
/// <param name="p">Permutation.</param>
/// <returns>The cost of the permutation.</returns>
public static double Evaluate(double[][] L, int[] p)
{
double cost = 0.0;
for (int i = 0; i < p.Length; i++)
{
cost += L[p];
}
return cost;
}
/// <summary>Calculate the reduced costs matrix for the given LAP and solution.<summary>
/// <param name="L">LAP cost matrix.</param>
/// <param name="U">Reduced cost matrix (filled in by the method)</param>
/// <param name="ys">Dual variables.</param>
/// <param name="yt">Dual variables.</param>
public static void ReducedCosts(
double[][] L,
double[][] U,
double[] ys,
double[] yt)
{
for (int i = 0; i < L.Length; i++)
{
for (int j = 0; j < L.Length; j++)
{
U = L - (yt + ys);
}
}
}
/// <summary>Solve the specified linear assignment problem.<summary>
/// <param name="L">LAP cost matrix.</param>
/// <param name="p">The solution permutation.</param>
/// <param name="ip">The inverse of the solution permutation.</param>
/// <param name="ys">Dual variables.</param>
/// <param name="yt">Dual variables.</param>
/// <returns></returns>
public static double Solve(
double[][] L,
int[] p,
int[] ip,
double[] ys,
double[] yt
)
{
if (L.Length == 0)
{
return 0.0;
}
return SolveCore(L, ip, p, ys, yt, 1.0e-7);
}
/// <summary>Linear sum assignment problem with double precision costs.<summary>
/// <param name="L">Cost matrix.</param>
/// <param name="ip">Inverse of optimal assignment.</param>
/// <param name="p">Optimal assignment.</param>
/// <param name="ys">Optimal dual (row) variables.</param>
/// <param name="yt">Optimal dual (column) variables.</param>
/// <param name="tolerance">Machine accuracy.</param>
/// <returns>The optimal cost.</returns>
private static double SolveCore(
double[][] L,
int[] ip,
int[] p,
double[] ys,
double[] yt,
double tolerance
)
{
const int BadIndex = -1;
int[] vor = new int[L.Length];
bool[] label = new bool[L.Length];
double[] dMinus = new double[L.Length];
double[] dPlus = new double[L.Length];
// Start: Construction of an initial (partial) assignment.
// (partial --> p/ip may be BadIndex for some i,j)
double maxValue = 1.0e15;
for (int i = 0; i < p.Length; i++)
{
ip = BadIndex;
p = BadIndex;
vor = BadIndex;
ys = 0;
yt = 0;
}
double ui = 0.0;
int jOpt = -1; // optimal j index
for (int i = 0; i < L.Length; i++)
{
for (int j = 0; j < L.Length; j++)
{
if ((j == 0) || (L - ui < tolerance))
{
ui = L;
jOpt = j;
}
}
// guarantee: jOpt >= 0
ys = ui;
if (ip == BadIndex)
{
ip = i;
p = jOpt;
}
}
for (int j = 0; j < ip.Length; j++)
{
if (ip == BadIndex)
{
yt = maxValue;
}
}
for (int i = 0; i < ys.Length; i++)
{
ui = ys;
for (int j = 0; j < yt.Length; j++)
{
double vj = yt;
if (vj > tolerance)
{
double cc = L - ui;
if (cc + tolerance < vj)
{
yt = cc;
vor = i;
}
}
}
}
for (int j = 0; j < vor.Length; j++)
{
int i = vor;
if ((i != BadIndex) && (p == BadIndex))
{
p = j;
ip = i;
}
}
for (int i = 0; i < p.Length; i++)
{
if (p == BadIndex)
{
ui = ys;
for (int j = 0; j < ip.Length; j++)
{
if (ip == BadIndex)
{
double lij = L;
if ((lij - ui - yt + tolerance) <= 0)
{
p = j;
ip = i;
}
}
}
}
}
// Construction of the optimal assignment
double d = 0.0;
int w = BadIndex;
int index = BadIndex;
for (int u = 0; u < p.Length; u++)
{
if (p == BadIndex)
{
// Shortest path computation
for (int i = 0; i < vor.Length; i++)
{
vor = u;
label = false;
dPlus = maxValue;
dMinus = L - ys - yt;
}
dPlus = 0;
bool done = false;
while (!done)
{
d = maxValue;
for (int i = 0; i < dMinus.Length; i++)
{
if (!label && (dMinus + tolerance < d))
{
d = dMinus;
index = i;
}
}
// note: index >= 0
if (ip != BadIndex)
{
label = true;
w = ip;
dPlus = d;
for (int i = 0; i < label.Length; i++)
{
if (!label)
{
double vgl = d + L - ys - yt;
if (dMinus > vgl + tolerance)
{
dMinus = vgl;
vor = w;
}
}
}
}
else
{
done = true;
}
}
// Augmentation
done = false;
while (!done)
{
w = vor;
ip = w;
int tempIndex = p;
p = index;
if (w == u)
{
done = true;
}
else
{
index = tempIndex;
}
}
// Transformation
for (int i = 0; i < dPlus.Length; i++)
{
if (dPlus != maxValue)
{
ys += d - dPlus;
}
if (dMinus + tolerance < d)
{
yt += dMinus - d;
}
}
}
}
// Computation of the optimal value
return Evaluate(L, p);
}
}
}