Pivot correctly when infeasible

README.md

@@ -9,7 +9,6 @@ a tutorial, by Simon Peyton Jones and David R Lester](http://research.microsoft.
 Limitations:
 
 - Only <= and >= constraints are supported
-- Origin must be feasible
 - Continuous variables only
 - No bounds expressions; 0 <= all vars <= infinity
 - No constants except on RHS of constraints

js/simplex.js

@@ -5,18 +5,29 @@ function getByVar(mat, colIdx, sym) {
     return mat.rows[colIdx][mat.varIndices[sym]];
 }
 
+function padValue(v) {
+    var s = '';
+    if (v < 10 && v > -10) {
+        s += ' ';
+    }
+    if (v >= 0) {
+        s += ' ';
+    }
+    return s + v;
+}
+
 function printMat(mat) {
-    var s = "";
+    var s = "\n";
     for (var c = 0; c < mat.vars.length; c++) {
-        s += mat.vars[c] + " ";
+        s += "  " + mat.vars[c] + " ";
     }
     s += "\n";
     for (var r = 0; r < mat.rows.length; r++) {
         var row = mat.rows[r];
         for (c = 0; c < mat.vars.length; c++) {
-            s += row[c] + " ";
+            s += padValue(row[c]) + " ";
         }
-        s += " : " + mat.rhs[r] + "\n";
+        s += " : " + padValue(mat.rhs[r]) + "\n";
     }
     return s;
 }
@@ -66,12 +77,33 @@ function firstInfeasibility(mat) {
             }
         }
         if (val < 0) {
-            return {
-                row: row,
-                col: c,
-                sym: mat.vars[c]
-            };
+            break;
+        } else {
+            row = null;
+        }
+    }
+
+    if (row != null) {
+        var col = null;
+        var max = null;
+        for (c = 0; c < mat.vars.length; c++) {
+            coef = mat.rows[row][c];
+            if (coef > 0) {
+                if (max === null || coef > max) {
+                    max = coef;
+                    col = c;
+                }
+            }
+        }
+
+        if (col === null) {
+            throw new Error('Failed to find a suitable pivot var for infeasible row ' + row);
         }
+        return {
+            row: row,
+            col: col,
+            sym: mat.vars[col]
+        };
     }
 
     return null;
@@ -80,15 +112,15 @@ function firstInfeasibility(mat) {
 function pivotVar(mat) {
     var p = {
         sym: mat.vars[0],
-        index: 0
+        col: 0
     };
     var obj = mat.rows[mat.rows.length - 1] ;
     var min = obj[0];
     for (var i = 1; i < mat.vars.length; i++) {
         if (obj[i] < min) {
             p = {
                 sym: mat.vars[i],
-                index: i
+                col: i
             }
             min = obj[i];
         }
@@ -105,7 +137,7 @@ function pivotRow(mat, pivotVar) {
     var min = null;
 
     for (var i = 0; i < mat.rows.length - 1; i++) {
-        var v = mat.rows[i][pivotVar.index];
+        var v = mat.rows[i][pivotVar.col];
         if (v <= 0) {
             continue;
         }

test/simplexSpec.js

@@ -106,7 +106,7 @@ describe('simplex', function () {
             var pv = pivotVar(mat);
             expect(pv).toEqual({
                 sym: 'b',
-                index: 1
+                col: 1
             });
         });
 
@@ -143,7 +143,7 @@ describe('simplex', function () {
 
             var pv = {
                 sym: 'a',
-                index: 0
+                col: 0
             };
             var pr = pivotRow(mat, pv);
             expect(pr).toBe(2);
@@ -162,7 +162,7 @@ describe('simplex', function () {
 
             var pv = {
                 sym: 'a',
-                index: 0
+                col: 0
             };
             var pr = pivotRow(mat, pv);
             expect(pr).toBe(null);
@@ -275,8 +275,8 @@ describe('simplex', function () {
                 rhs: [1, 1, 1, 1]
             };
 
-            var row = firstInfeasibility(mat);
-            expect(row).toBe(null);
+            var inf = firstInfeasibility(mat);
+            expect(inf).toBe(null);
         });
 
         it('should return the first infeasible row', function () {
@@ -285,18 +285,18 @@ describe('simplex', function () {
                 varIndices: {},
                 rows: [
                     [1,  0,  4,  0],
-                    [1, -3,  0,  0],
+                    [1, -3,  2,  0],
                     [1,  0,  0, -1],
                     [1,  0, -8,  0]
                 ],
                 rhs: [1, 1, 1, 1]
             };
 
-            var row = firstInfeasibility(mat);
-            expect(row).toEqual({
+            var inf = firstInfeasibility(mat);
+            expect(inf).toEqual({
                 row: 1,
-                col: 1,
-                sym: 'b'
+                col: 2,
+                sym: 'c'
             });
         });
     });
@@ -322,13 +322,47 @@ describe('simplex', function () {
                     break;
                 }
                 var pr = pivotRow(mat, pv);
-                mat = pivot(mat, pr, pv.index);
+                mat = pivot(mat, pr, pv.col);
             }
 
             var sol = solution(mat);
             expect(sol.x).toBe(30);
             expect(sol.y).toBe(40);
             expect(sol._obj).toBeCloseTo(410);
         });
+
+        it('should solve a problem which is initially infeasible', function () {
+            var toks = clex(fromArray(
+                "maximise\n\
+                2x + 3y + z\n\
+                subject to\n\
+                x + y + z <= 40\n\
+                2x + y - z >= 10\n\
+                -y + z >= 10\n\
+                end\n"
+            ));
+            var result = pLp(toks);
+            var prob = takeFirstParse(result);
+            var canon = canonical(prob);
+            var mat = toMatrix(canon);
+
+            while (true) {
+                var pv = firstInfeasibility(mat);
+                if (pv === null) {
+                    pv = pivotVar(mat);
+                }
+                if (pv === null) {
+                    break;
+                }
+                var pr = pivotRow(mat, pv);
+                mat = pivot(mat, pr, pv.col);
+            }
+
+            var sol = solution(mat);
+            expect(sol.x).toBe(10);
+            expect(sol.y).toBe(10);
+            expect(sol.z).toBe(20);
+            expect(sol._obj).toBe(70);
+        });
     });
 });