added handling of general sparse constraints

3 files changed, 41 insertions, 6 deletions

diff --git a/packages/glpk/examples/simplex2.hs b/packages/glpk/examples/simplex2.hs
index e9e8859..0d83ca9 100644
--- a/packages/glpk/examples/simplex2.hs
+++ b/packages/glpk/examples/simplex2.hs
@@ -10,9 +10,14 @@ constr2 = Dense [ [2,1,0] :<=: 10
                 , [0,1,5] :<=: 20
                 ]
 
+constr3 = General [ [1#1, 1#1, 1#2] :<=: 10
+                  , [1#2, 5#3] :<=: 20
+                  ]
+
 main = do
     print $ simplex prob constr1 []
     print $ simplex prob constr2 []
+    print $ simplex prob constr3 []
     print $ simplex prob constr2 [ 2 :>=: 1, 3 :&: (2,7)]
     print $ simplex prob constr2 [ Free 2 ]
 
diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal
index cd761e0..229197f 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -22,7 +22,7 @@ extra-source-files:     examples/simplex1.hs
                         examples/simplex4.hs
 
 library
-    Build-Depends:      base <5, hmatrix >= 0.16
+    Build-Depends:      base <5, hmatrix >= 0.16, containers >= 0.5.4.0
 
     hs-source-dirs:     src
 
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
index b0537cc..a54eb59 100644
--- a/packages/glpk/src/Numeric/LinearProgramming.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
@@ -49,6 +49,14 @@ constr2 = Dense [ [2,1,0] :<=: 10
                 ]
 @
 
+Note that when using sparse constraints, coefficients cannot appear more than once in each constraint. You can alternatively use General which will automatically sum any duplicate coefficients when necessary.
+
+@
+constr3 = General [ [1\#1, 1\#1, 1\#2] :<=: 10
+                  , [1\#2, 5\#3] :<=: 20
+                  ]
+@
+
 By default all variables are bounded as @x_i >= 0@, but this can be
 changed:
 
@@ -67,6 +75,7 @@ Multiple bounds for a variable are not allowed, instead of
 
 module Numeric.LinearProgramming(
     simplex,
+    sparseOfGeneral,
     Optimization(..),
     Constraints(..),
     Bounds,
@@ -82,13 +91,14 @@ import System.IO.Unsafe(unsafePerformIO)
 import Foreign.C.Types
 import Data.List((\\),sortBy,nub)
 import Data.Function(on)
+import qualified Data.Map.Strict as Map
 
 --import Debug.Trace
 --debug x = trace (show x) x
 
 -----------------------------------------------------
 
--- | Coefficient of a variable for a sparse representation of constraints.
+-- | Coefficient of a variable for a sparse and general representations of constraints.
 (#) :: Double -> Int -> (Double,Int)
 infixl 5 #
 (#) = (,)
@@ -108,18 +118,29 @@ data Solution = Undefined
               | Unbounded
               deriving Show
 
-data Constraints = Dense  [ Bound [Double] ]
-                 | Sparse [ Bound [(Double,Int)] ]
+data Constraints = Dense   [ Bound [Double] ]
+                 | Sparse  [ Bound [(Double,Int)] ]
+                 | General [ Bound [(Double,Int)] ]
 
 data Optimization = Maximize [Double]
                   | Minimize [Double]
 
 type Bounds = [Bound Int]
 
+-- | Convert a system of General constraints to one with unique coefficients. 
+sparseOfGeneral :: Constraints -> Constraints
+sparseOfGeneral (General cs) =
+    Sparse $ map (\bl -> 
+                      let cl = obj bl in
+                      let m = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in
+                      withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] m)) cs
+sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected"
+
 simplex :: Optimization -> Constraints -> Bounds -> Solution
 
-simplex opt (Dense  []) bnds = simplex opt (Sparse []) bnds
-simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+simplex opt (Dense   []) bnds = simplex opt (Sparse []) bnds
+simplex opt (Sparse  []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
+simplex opt (General []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
 
 simplex opt (Dense constr) bnds = extract sg sol where
     sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
@@ -133,6 +154,8 @@ simplex opt (Sparse constr) bnds = extract sg sol where
     m = length constr
     (sz, sg, objfun) = adapt opt
 
+simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds
+
 adapt :: Optimization -> (Int, Double, [Double])
 adapt opt = case opt of
     Maximize x -> (size x, 1 ,x)
@@ -162,6 +185,13 @@ obj (x :&: _)  = x
 obj (x :==: _) = x
 obj (Free x)   = x
 
+withObj :: Bound t -> t -> Bound t
+withObj (_ :<=: b) x = (x :<=: b)
+withObj (_ :>=: b) x = (x :>=: b)
+withObj (_ :&: b) x = (x :&: b)
+withObj (_ :==: b) x = (x :==: b)
+withObj (Free _) x = Free x
+
 tb :: Bound t -> Double
 tb (_ :<=: _)  = glpUP
 tb (_ :>=: _)  = glpLO