Simple matrix, vector, permutation utilities

I am going to be rolling out the rest of my branch-and-bound algorithm in the next few posts.  To make that easier, in this post I introduce some common matrix, vector, and permutation methods.  It turns out that for technical computing applications, C#’s extension methods (introduced in 3.0) are awesome.  With vectors it’s great because you retain the control of having direct array access, but you get the nice object-oriented notation.

I’ve left out comments and error-checking, and it’s not super-optimized, but you get the idea.  None of these methods end up being the bottleneck.  Most of the methods are obvious, but I do want to comment on the permutation methods.  The goal of QAP is to find an optimal assignment of facilities to locations, represented as a permutation.  In the course of our branch-and-bound algorithm, we’ll be making partial assignments – that is, only some of the entries in the permutation will be filled in.  By convention, p[i] == -1 will mean that facility i is unassigned.  An index where p[i] < 0 is called “unused”.  When we branch, we’ll want to pick an unused facility (or location), and try assigning all unused locations (or facilities) to it.   The use of C#’s iterator and yield concepts really help here.

That said, here’s the code.

public static class MatrixUtilities {
  #region Vector
  public static void ConstantFill(this T[] data, T val) 
{
for (int i = 0; i < data.Length; i++) {
data[i] = val;
}
}

public static int ArgMin(this double[] data, out double best) {
if (data.Length == 0) {
best = Double.MinValue;
return -1;
}
int iBest = 0;
best = data[0];
for (int i = 1; i < data.Length; i++) {
if (best > data[i]) {
best = data[i];
iBest = i;
}
}
return iBest;
}
  public static int ArgMax(this double[] data, out double best) { 
if (data.Length == 0) { best = Double.MinValue; return -1; }
int iBest = 0;
best = data[0];
for (int i = 1; i < data.Length; i++) {
if (best < data[i]) {
best = data[i];
iBest = i;
}
}
return iBest;
}
#endregion

#region Permutation
public static void Swap(this int[] p, int i, int j) {
int temp = p[i];
p[i] = p[j];
p[j] = temp;
}

public static int FindUnused(this int[] data, int index) {
foreach (int unused in data.UnusedIndices()) {
if (index-- <= 0) {
return unused;
}
}
return -1;
}

public static IEnumerable UnusedIndices(this int[] data) {
for (int i = 0; i < data.Length; i++) {
if (data[i] < 0) {
yield return i;
}
}
}

public static string PermutationToString(this int[] p, bool oneBased) {
StringBuilder build = new StringBuilder(p.Length * 4);
int width = ((int)Math.Log10(p.Length)) + 2;
build.Append("[");
for (int i = 0; i < p.Length; i++) {
int index = oneBased ? p[i] + 1 : p[i];
build.Append(index.ToString().PadLeft(width));
}
build.Append("]");
return build.ToString();
}

#endregion

#region Matrix
public static string MatrixToString(this double[][] A) {
if (A != null) {
StringBuilder build = new StringBuilder(A.Length * A.Length * 3);
for (int i = 0; i < A.Length; i++) {
if (A[i] != null) {
for (int j = 0; j < A[i].Length; j++) {
build.AppendFormat("{0,4}", A[i][j]);
build.Append(" ");
}
build.AppendLine();
}
}
return build.ToString();
}
return null;
}

public static double[][] NewMatrix(int m, int n) {
double[][] M = new double[m][];
for (int i = 0; i < M.Length; i++) {
M[i] = new double[n];
}
return M;
}

public static double[][] NewMatrix(int n) {
return NewMatrix(n, n);
}
#endregion
}

Author: natebrix

Follow me on twitter at @natebrix.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s