Simplex maximization algorithm in C#

simplex.cs

using System;
using System.Diagnostics;
using System.Collections.Generic;

namespace Simplex {
  class Simplex {
    private double[] c;
    private double[,] A;
    private double[] b;
    private HashSet<int> N = new HashSet<int>();
    private HashSet<int> B = new HashSet<int>();
    private double v = 0;

    public Simplex(double[] c, double[,] A, double[] b) {
      int vars = c.Length, constraints = b.Length;

      if (vars != A.GetLength(1)) {
        throw new Exception("Number of variables in c doesn't match number in A.");
      }

      if (constraints != A.GetLength(0)) {
        throw new Exception("Number of constraints in A doesn't match number in b.");
      }

      // Extend max fn coefficients vector with 0 padding
      this.c = new double[vars+constraints];
      Array.Copy(c, this.c, vars);

      // Extend coefficient matrix with 0 padding
      this.A = new double[vars+constraints,vars+constraints];
      for (int i = 0; i < constraints; i++) {
        for (int j = 0; j < vars; j++) {
          this.A[i+vars, j] = A[i,j];
        }
      }

      // Extend constraint right-hand side vector with 0 padding
      this.b = new double[vars+constraints];
      Array.Copy(b, 0, this.b, vars, constraints);

      // Populate non-basic and basic sets
      for (int i = 0; i < vars; i++) {
        N.Add(i);
      }

      for (int i = 0; i < constraints; i++) {
        B.Add(vars + i);
      }
    }

    public Tuple<double, double[]> maximize() {
      while (true) {
        // Find highest coefficient for entering var
        int e = -1;
        double ce = 0;
        foreach (var _e in N) {
          if (c[_e] > ce) {
            ce = c[_e];
            e = _e;
          }
        }

        // If no coefficient > 0, there's no more maximizing to do, and we're almost done
        if (e == -1) break;

        // Find lowest check ratio
        double minRatio = double.PositiveInfinity;
        int l = -1;
        foreach (var i in B) {
          if (A[i, e] > 0) {
            double r = b[i] / A[i, e];
            if (r < minRatio) {
              minRatio = r;
              l = i;
            }
          }
        }

        // Unbounded
        if (double.IsInfinity(minRatio)) {
          return Tuple.Create<double, double[]>(double.PositiveInfinity, null);
        }

        pivot(e, l);
      }

      // Extract amounts and slack for optimal solution
      double[] x = new double[b.Length];
      int n = b.Length;
      for (var i = 0; i < n; i++) {
        x[i] = B.Contains(i) ? b[i] : 0;
      }

      // Return max and variables
      return Tuple.Create<double, double[]>(v, x);
    }

    private void pivot(int e, int l) {
      N.Remove(e);
      B.Remove(l);

      b[e] = b[l] / A[l, e];

      foreach (var j in N) {
        A[e, j] = A[l, j] / A[l, e];
      }

      A[e, l] = 1 / A[l, e];

      foreach (var i in B) {
        b[i] = b[i] - A[i, e] * b[e];

        foreach (var j in N) {
          A[i, j] = A[i, j] - A[i, e] * A[e, j];
        }

        A[i, l] = -1 * A[i, e] * A[e, l];
      }

      v = v + c[e] * b[e];

      foreach (var j in N) {
        c[j] = c[j] - c[e] * A[e, j];
      }

      c[l] = -1 * c[e] * A[e, l];

      N.Add(l);
      B.Add(e);
    }
  }

  class MainClass {
    static void Main(string[] args) {
      var s = new Simplex(
        new [] {10.2, 422.3, 6.91, 853},
        new [,] {
          {0.1, 0.5, 0.333333, 1},
          {30, 15, 19, 12},
          {1000, 6000, 4100, 9100},
          {50, 20, 21, 10},
          {4, 6, 19, 30}
        },
        new double[] {2000, 1000, 1000000, 640, 432}
      );

      var answer = s.maximize();
      Console.WriteLine(answer.Item1);
      Console.WriteLine(string.Join(", ", answer.Item2));
    }
  }
}

It doesn't seem to support negative constraints (b). Is there a workaround?

For example this works:

// We are asked to maximize revenue where
//   Revenue = 45*numChairs + 80*numTables
//   Subject to constraints:
//     5*numChairs + 20*numTables <= 400 wood units
//     10*numChairs + 15*numTables <= 450 labor hours
//     numChairs + numTables <= 35 as per Soviet rules
//   Implicit: numChairs >= 0 && numTables >= 0
var s = new Simplex(
    new double[] { 45, 80 }, // revenue for each product
    new double[,] { // constraint coefficients
        { 5, 20 },
        { 10, 15 },
        { 1, 1 },
    },
    new double[] { 400, 450, 35 } // limits
);
// Solution: (maximum: 2100, production: [20, 15], remaining numbers: [0, 25, 0])

but this doesn't work:

// We are asked to maximize revenue where
//   Revenue = 45*numChairs + 80*numTables
//   Subject to constraints:
//     5*numChairs + 20*numTables <= 400 wood units
//     10*numChairs + 15*numTables <= 450 labor hours
//     numChairs + numTables >= 40 as per Soviet rules
//   Implicit: numChairs >= 0 && numTables >= 0
var s = new Simplex(
    new double[] { 45, 80 }, // revenue for each product
    new double[,] { // constraint coefficients
        { 5, 20 },
        { 10, 15 },
        { -1, -1 },
    },
    new double[] { 400, 450, -40 } // limits
);
// Incorrect solution: (maximum: 2200, production: [24, 14], remaining numbers: [0, 0, -2])

Edit: it looks like the author has only implemented "one-phase" simplex, but support for negative limits (or >= inequalities, which are equivalent) requires "two-phase simplex".