Home > Optimization, QAP > Simple matrix, vector, permutation utilities

## 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 } `