error checking for simplex

5 files changed, 74 insertions, 64 deletions

diff --git a/packages/glpk/examples/simplex1.hs b/packages/glpk/examples/simplex1.hs
index 9639f37..ec72bdc 100644
--- a/packages/glpk/examples/simplex1.hs
+++ b/packages/glpk/examples/simplex1.hs
@@ -5,17 +5,16 @@ import Numeric.LinearProgramming
 objFun = Maximize [10, 6, 4]
 
 constr = Dense [ [1,1,1]  :<: 100
-               , [10,4,5] :<: 600 
+               , [10,4,5] :<: 600
                , [2,2,6]  :<: 300 ]
 
 -- default bounds
 bnds = [ 1 :>: 0
        , 2 :>: 0
        , 3 :>: 0 ]
-         
+
 main = do
     print $ simplex objFun constr []
     print $ simplex objFun constr bnds
     print $ simplex objFun constr [Free 3]
     print $ simplex objFun constr [ 2 :<: 50 ]
-
diff --git a/packages/glpk/examples/simplex2.hs b/packages/glpk/examples/simplex2.hs
index 40e4bbd..0b865a6 100644
--- a/packages/glpk/examples/simplex2.hs
+++ b/packages/glpk/examples/simplex2.hs
@@ -1,17 +1,18 @@
 import Numeric.LinearProgramming
 
-prob = Maximize [4, 3, -2, 7]
+prob = Maximize [4, -3, 2]
 
-constr1 = Sparse [ [1#1, 1#2] :<: 10
-                 , [1#3, 1#4] :<: 10 
+constr1 = Sparse [ [2#1, 1#2] :<: 10
+                 , [1#2, 5#3] :<: 20
                  ]
 
-constr2 = Dense [ [1,1,0,0] :<: 10
-                , [0,0,1,1] :<: 10 
+constr2 = Dense [ [2,1,0] :<: 10
+                , [0,1,5] :<: 20
                 ]
 
 main = do
     print $ simplex prob constr1 []
-    print $ simplex prob constr2 [ 2 :>: 1, 4 :&: (2,7)]
-    print $ simplex prob constr2 [ Free 3 ]
+    print $ simplex prob constr2 []
+    print $ simplex prob constr2 [ 2 :>: 1, 3 :&: (2,7)]
+    print $ simplex prob constr2 [ Free 2 ]
 
diff --git a/packages/glpk/examples/simplex3.hs b/packages/glpk/examples/simplex3.hs
index 3a7e8e8..1df6b28 100644
--- a/packages/glpk/examples/simplex3.hs
+++ b/packages/glpk/examples/simplex3.hs
@@ -7,12 +7,12 @@ prob = Minimize [0.03, 0.08, 0.17, 0.12, 0.15, 0.21, 0.38]
 
 constr = Dense
     [ [1,1,1,1,1,1,1]                           :==: 2000
-    , [0.15, 0.04, 0.02, 0.04, 0.2,0.01, 0.03]   :<:  60 
-    , [0.03, 0.05, 0.08, 0.02, 0.06, 0.01, 0]    :<:  100 
-    , [0.02, 0.04, 0.01, 0.02, 0.02, 0,    0]    :<:  40 
+    , [0.15, 0.04, 0.02, 0.04, 0.2,0.01, 0.03]   :<:  60
+    , [0.03, 0.05, 0.08, 0.02, 0.06, 0.01, 0]    :<:  100
+    , [0.02, 0.04, 0.01, 0.02, 0.02, 0,    0]    :<:  40
     , [0.02, 0.03, 0,    0,    0.01, 0,    0]    :<:  30
-    , [0.7,  0.75, 0.8,  0.75, 0.8,  0.97,  0]   :>:  1500 
-    , [0.02, 0.06, 0.08, 0.12, 0.02, 0.01, 0.97] :&: (250,300) 
+    , [0.7,  0.75, 0.8,  0.75, 0.8,  0.97,  0]   :>:  1500
+    , [0.02, 0.06, 0.08, 0.12, 0.02, 0.01, 0.97] :&: (250,300)
     ]
 
 bounds = [ 1 :&: (0,200)
diff --git a/packages/glpk/examples/simplex4.hs b/packages/glpk/examples/simplex4.hs
index 8496f6d..dd65002 100644
--- a/packages/glpk/examples/simplex4.hs
+++ b/packages/glpk/examples/simplex4.hs
@@ -7,12 +7,12 @@ prob = Minimize [0.03, 0.08, 0.17, 0.12, 0.15, 0.21, 0.38]
 
 constr = Sparse
     [ [1#1,1#2,1#3,1#4,1#5,1#6,1#7]                          :==: 2000
-    , [0.15#1, 0.04#2, 0.02#3, 0.04#4, 0.2#5,0.01#6, 0.03#7] :<: 60 
-    , [0.03#1, 0.05#2, 0.08#3, 0.02#4, 0.06#5, 0.01#6]       :<: 100 
-    , [0.02#1, 0.04#2, 0.01#3, 0.02#4, 0.02#5]               :<:  40 
+    , [0.15#1, 0.04#2, 0.02#3, 0.04#4, 0.2#5,0.01#6, 0.03#7] :<: 60
+    , [0.03#1, 0.05#2, 0.08#3, 0.02#4, 0.06#5, 0.01#6]       :<: 100
+    , [0.02#1, 0.04#2, 0.01#3, 0.02#4, 0.02#5]               :<:  40
     , [0.02#1, 0.03#2, 0.01#5]                               :<:  30
-    , [0.7#1,  0.75#2, 0.8#3,  0.75#4, 0.8#5,  0.97#6]       :>:  1500 
-    , [0.02#1, 0.06#2, 0.08#3, 0.12#4, 0.02#5, 0.01#6, 0.97#7] :&: (250,300) 
+    , [0.7#1,  0.75#2, 0.8#3,  0.75#4, 0.8#5,  0.97#6]       :>:  1500
+    , [0.02#1, 0.06#2, 0.08#3, 0.12#4, 0.02#5, 0.01#6, 0.97#7] :&: (250,300)
     ]
 
 bounds = [ 1 :&: (0,200)
diff --git a/packages/glpk/lib/Numeric/LinearProgramming.hs b/packages/glpk/lib/Numeric/LinearProgramming.hs
index fff304f..a888214 100644
--- a/packages/glpk/lib/Numeric/LinearProgramming.hs
+++ b/packages/glpk/lib/Numeric/LinearProgramming.hs
@@ -10,13 +10,13 @@ Stability   :  provisional
 
 This module provides an interface to the standard simplex algorithm.
 
-For example, the following linear programming problem
+For example, the following LP problem
 
-@maximize 4 x_1 + 3 x_2 - 2 x_3 + 7 x_4
+@maximize 4 x_1 - 3 x_2 + 2 x_3
 subject to
 
-x_1 + x_2 <= 10
-x_3 + x_4 <= 10
+2 x_1 +   x_2 <= 10
+  x_3 + 5 x_4 <= 20
 
 and
 
@@ -26,31 +26,33 @@ can be solved as follows:
 
 @import Numeric.LinearProgramming
 
-prob = Maximize [4, 3, -2, 7]
+prob = Maximize [4, -3, 2]
 
-constr1 = Sparse [ [1\#1, 1\#2] :<: 10
-                 , [1\#3, 1\#4] :<: 10 
+constr1 = Sparse [ [2#1, 1#2] :<: 10
+                 , [1#2, 5#3] :<: 20
                  ]
-                 
 
 \> simplex prob constr1 []
-Optimal (110.0,[10.0,0.0,0.0,10.0])@
+Optimal (28.0,[5.0,0.0,4.0])@
 
 The coefficients of the constraint matrix can also be given in dense format:
 
-@constr2 = Dense [ [1,1,0,0] :<: 10
-                , [0,0,1,1] :<: 10 
+@constr2 = Dense [ [2,1,0] :<: 10
+                , [0,1,5] :<: 20
                 ]@
 
-By default all variables are bounded as @x_i <= 0@, but this can be
+By default all variables are bounded as @x_i >= 0@, but this can be
 changed:
 
-@\> simplex prob constr2 [2 :>: 1, 4 :&: (2,7)]
-Optimal (88.0,[9.0,1.0,0.0,7.0])
+@\> simplex prob constr2 [ 2 :>: 1, 3 :&: (2,7)]
+Optimal (22.6,[4.5,1.0,3.8])
 
-\> simplex prob constr2 [Free 3]
+\> simplex prob constr2 [Free 2]
 Unbounded@
 
+The given bound for a variable completely replaces the default,
+so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@.
+
 -}
 
 module Numeric.LinearProgramming(
@@ -94,7 +96,7 @@ data Solution = Undefined
               | Optimal (Double, [Double])
               | Unbounded
               deriving Show
-            
+
 data Constraints = Dense  [ Bound [Double] ]
                  | Sparse [ Bound [(Double,Int)] ]
 
@@ -105,19 +107,25 @@ type Bounds = [Bound Int]
 
 simplex :: Optimization -> Constraints -> Bounds -> Solution
 
+simplex opt (Dense  []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+
 simplex opt (Dense constr) bnds = extract sg sol where
-    sol = simplexDense (mkConstrD objfun constr) (mkBoundsD constr bnds)
-    (sg, objfun) = case opt of
-        Maximize x -> (1 ,x)
-        Minimize x -> (-1, (map negate x))
+    sol = simplexDense (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
+    (sz, sg, objfun) = adapt opt
 
 simplex opt (Sparse constr) bnds = extract sg sol where
-    sol = simplexSparse m n (mkConstrS objfun constr) (mkBoundsS constr bnds)
+    sol = simplexSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)
     n = length objfun
     m = length constr
-    (sg, objfun) = case opt of
-        Maximize x -> (1 ,x)
-        Minimize x -> (-1, (map negate x))
+    (sz, sg, objfun) = adapt opt
+
+adapt :: Optimization -> (Int, Double, [Double])
+adapt opt = case opt of
+    Maximize x -> (size x, 1 ,x)
+    Minimize x -> (size x, -1, (map negate x))
+ where size x | null x = error "simplex: objective function with zero variables"
+              | otherwise = length x
 
 extract :: Double -> Vector Double -> Solution
 extract sg sol = r where
@@ -168,31 +176,33 @@ mkBound1 b = [tb b, lb b, ub b]
 mkBound2 :: Bound t -> (t, [Double])
 mkBound2 b = (obj b, mkBound1 b)
 
-mkBoundsD :: [Bound [a]] -> [Bound Int] -> Matrix Double
-mkBoundsD b1 b2 = fromLists (cb++vb) where
-    c = length (obj (head b1))
+mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double
+mkBounds n b1 b2 = fromLists (cb++vb) where
     gv = map obj b2
-    rv = [1..c] \\ gv
+    rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv
+       | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n
     vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2
     cb = map mkBound1 b1
 
-mkConstrD :: [Double] -> [Bound [Double]] -> Matrix Double
-mkConstrD f b1 = fromLists (f : map obj b1)
-
-mkBoundsS :: [Bound [(Double, Int)]] -> [Bound Int] -> Matrix Double
-mkBoundsS b1 b2 = fromLists (cb++vb) where
-    c = maximum $ map snd $ concatMap obj b1
-    gv = map obj b2
-    rv = [1..c] \\ gv
-    vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2
-    cb = map mkBound1 b1
-
-mkConstrS :: [Double] -> [Bound [(Double, Int)]] -> Matrix Double
-mkConstrS objfun constr = fromLists (ob ++ co) where
+mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double
+mkConstrD n f b1 | ok = fromLists (f : cs)
+                 | otherwise = error $ "simplex: dense constraints require "++show n
+                                     ++" variables, given " ++ show ls
+    where
+       cs | null b1   = error $ "simplex: dense: empty constraints"
+          | otherwise = map obj b1
+       ls = map length cs
+       ok = all (==n) ls
+
+mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double
+mkConstrS n objfun b1 = fromLists (ob ++ co) where
     ob = map (([0,0]++).return) objfun
-    co = concat $ zipWith f [1::Int ..] (map obj constr)
+    co = concat $ zipWith f [1::Int ..] cs
+    cs    | null b1   = error $ "simplex: sparse: empty constraints"
+          | otherwise = map obj b1
     f k = map (g k)
-    g k (c,v) = [fromIntegral k, fromIntegral v,c]
+    g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c]
+              | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n
 
 -----------------------------------------------------
 
@@ -262,6 +272,6 @@ bounds = (6><3)
   , glpLO,0,0
   , glpLO,0,0
   , glpLO,0,0 ]
-  
+
 -}