simple glpk interface

9 files changed, 484 insertions, 0 deletions

diff --git a/packages/glpk/LICENSE b/packages/glpk/LICENSE
new file mode 100644
index 0000000..f2125ec
--- /dev/null
+++ b/packages/glpk/LICENSE
@@ -0,0 +1,2 @@
+Copyright Alberto Ruiz 2010
+GPL license
diff --git a/packages/glpk/Setup.lhs b/packages/glpk/Setup.lhs
new file mode 100644
index 0000000..6b32049
--- /dev/null
+++ b/packages/glpk/Setup.lhs
@@ -0,0 +1,4 @@
+#! /usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
diff --git a/packages/glpk/examples/simplex1.hs b/packages/glpk/examples/simplex1.hs
new file mode 100644
index 0000000..4609524
--- /dev/null
+++ b/packages/glpk/examples/simplex1.hs
@@ -0,0 +1,21 @@
+-- first example in glpk manual
+
+import Numeric.LinearProgramming
+
+objFun = Maximize [10, 6, 4]
+
+constr = Dense [ [1,1,1]  :<: 100
+              , [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
new file mode 100644
index 0000000..76a53df
--- /dev/null
+++ b/packages/glpk/examples/simplex2.hs
@@ -0,0 +1,16 @@
+import Numeric.LinearProgramming
+
+prob = Maximize [1,2,3,4]
+
+constr1 = Sparse [ [1#1, 1#2] :<: 10
+                 , [1#3, 1#4] :<: 10 
+                 ]
+
+constr2 = Dense [ [1,1,0,0] :<: 10
+                , [0,0,1,1] :<: 10 
+                ]
+
+main = do
+    print $ simplex prob constr1 []
+    print $ simplex prob constr2 []
+
diff --git a/packages/glpk/examples/simplex3.hs b/packages/glpk/examples/simplex3.hs
new file mode 100644
index 0000000..0f787cb
--- /dev/null
+++ b/packages/glpk/examples/simplex3.hs
@@ -0,0 +1,24 @@
+import Numeric.LinearProgramming
+
+-- $ glpsol --cpxlp  /usr/share/doc/glpk-utils/examples/plan.lp  -o result.txt
+
+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.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) 
+    ]
+
+bounds = [ 1 :&: (0,200)
+         , 2 :&: (0,2500)
+         , 3 :&: (400,800)
+         , 4 :&: (100,700)
+         , 5 :&: (0,1500) ]
+
+main = print $ simplex prob constr bounds
+
diff --git a/packages/glpk/examples/simplex4.hs b/packages/glpk/examples/simplex4.hs
new file mode 100644
index 0000000..3b9c060
--- /dev/null
+++ b/packages/glpk/examples/simplex4.hs
@@ -0,0 +1,24 @@
+import Numeric.LinearProgramming
+
+-- $ glpsol --cpxlp  /usr/share/doc/glpk-utils/examples/plan.lp  -o result.txt
+
+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.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) 
+    ]
+
+bounds = [ 1 :&: (0,200)
+         , 2 :&: (0,2500)
+         , 3 :&: (400,800)
+         , 4 :&: (100,700)
+         , 5 :&: (0,1500) ]
+
+main = print $ simplex prob constr bounds
+
diff --git a/packages/glpk/glpk.cabal b/packages/glpk/glpk.cabal
new file mode 100644
index 0000000..5df45ca
--- /dev/null
+++ b/packages/glpk/glpk.cabal
@@ -0,0 +1,35 @@
+Name:               hmatrix-glpk
+Version:            0.1.0
+License:            GPL
+License-file:       LICENSE
+Author:             Alberto Ruiz
+Maintainer:         Alberto Ruiz <aruiz@um.es>
+Stability:          experimental
+Homepage:           http://code.haskell.org/hmatrix
+Synopsis:           Linear Programming based on GLPK
+Description:
+ Simple interface to linear programming functions provided by GLPK.
+
+Category:           Math
+tested-with:        GHC ==6.10.4
+
+cabal-version:      >=1.2
+build-type:         Simple
+
+extra-source-files:     examples/simplex1.hs
+                        examples/simplex2.hs
+                        examples/simplex3.hs
+                        examples/simplex4.hs
+
+library
+    Build-Depends:      base >= 3 && < 5, hmatrix >= 0.8.3
+
+    hs-source-dirs:     lib
+
+    Exposed-modules:    Numeric.LinearProgramming
+
+    c-sources:          lib/Numeric/LinearProgramming/glpk.c
+
+    ghc-options:        -Wall
+
+    extra-libraries:    glpk
diff --git a/packages/glpk/lib/Numeric/LinearProgramming.hs b/packages/glpk/lib/Numeric/LinearProgramming.hs
new file mode 100644
index 0000000..1b14080
--- /dev/null
+++ b/packages/glpk/lib/Numeric/LinearProgramming.hs
@@ -0,0 +1,211 @@
+{-# LANGUAGE ForeignFunctionInterface #-}
+
+module Numeric.LinearProgramming(
+    Optimization(..),
+    Bound(..),
+    (#),
+    Coeffs(..),
+    Solution(..),
+    simplex,
+) where
+
+import Numeric.LinearAlgebra
+import Data.Packed.Development
+import Foreign(Ptr,unsafePerformIO)
+import Foreign.C.Types(CInt)
+import Data.List((\\),sortBy)
+import Data.Function(on)
+
+--import Debug.Trace
+--debug x = trace (show x) x
+
+-----------------------------------------------------
+
+-- | Coefficient of a variable
+(#) :: Double -> Int -> (Double,Int)
+infixl 5 #
+(#) = (,)
+
+data Bound x =  x :<: Double 
+             |  x :>: Double
+             |  x :&: (Double,Double)
+             |  x :==: Double
+             |  Free x
+             deriving Show
+
+data Solution = Undefined
+              | Feasible (Double, [Double])
+              | Infeasible (Double, [Double])
+              | NoFeasible
+              | Optimal (Double, [Double])
+              | Unbounded
+              deriving Show
+            
+data Coeffs = Dense  [ Bound [Double] ]
+            | Sparse [ Bound [(Double,Int)] ]
+
+data Optimization = Maximize [Double]
+                  | Minimize [Double]
+
+simplex :: Optimization -> Coeffs -> [Bound Int] -> Solution
+
+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))
+
+simplex opt (Sparse constr) bnds = extract sg sol where
+    sol = simplexSparse m n (mkConstrS objfun constr) (mkBoundsS constr bnds)
+    n = length objfun
+    m = length constr
+    (sg, objfun) = case opt of
+        Maximize x -> (1 ,x)
+        Minimize x -> (-1, (map negate x))
+
+extract :: Double -> Vector Double -> Solution
+extract sg sol = r where
+    z = sg * (sol@>1)
+    v = toList $ subVector 2 (dim sol -2) sol
+    r = case round(sol@>0)::Int of
+          1 -> Undefined
+          2 -> Feasible (z,v)
+          3 -> Infeasible (z,v)
+          4 -> NoFeasible
+          5 -> Optimal (z,v)
+          6 -> Unbounded
+          _ -> error "simplex: solution type unknown"
+
+-----------------------------------------------------
+
+obj :: Bound t -> t
+obj (x :<: _)  = x
+obj (x :>: _)  = x
+obj (x :&: _)  = x
+obj (x :==: _) = x
+obj (Free x)   = x
+
+tb :: Bound t -> Double
+tb (_ :<: _)  = glpUP
+tb (_ :>: _)  = glpLO
+tb (_ :&: _)  = glpDB
+tb (_ :==: _) = glpFX
+tb (Free _)   = glpFR
+
+lb :: Bound t -> Double
+lb (_ :<: _)     = 0
+lb (_ :>: a)     = a
+lb (_ :&: (a,_)) = a
+lb (_ :==: a)    = a
+lb (Free _)      = 0
+
+ub :: Bound t -> Double
+ub (_ :<: a)     = a
+ub (_ :>: _)     = 0
+ub (_ :&: (_,a)) = a
+ub (_ :==: a)    = a
+ub (Free _)      = 0
+
+mkBound1 :: Bound t -> [Double]
+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))
+    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
+
+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
+    ob = map (([0,0]++).return) objfun
+    co = concat $ zipWith f [1::Int ..] (map obj constr)
+    f k = map (g k)
+    g k (c,v) = [fromIntegral k, fromIntegral v,c]
+
+-----------------------------------------------------
+
+foreign import ccall "c_simplex_dense" c_simplex_dense
+    :: CInt -> CInt -> Ptr Double    -- coeffs
+    -> CInt -> CInt -> Ptr Double    -- bounds
+    -> CInt -> Ptr Double            -- result
+    -> IO CInt                       -- exit code
+
+simplexDense :: Matrix Double -> Matrix Double -> Vector Double
+simplexDense c b = unsafePerformIO $ do
+    s <- createVector (2+cols c)
+    app3 c_simplex_dense mat (cmat c) mat (cmat b) vec s "c_simplex_dense"
+    return s
+
+foreign import ccall "c_simplex_sparse" c_simplex_sparse
+    :: CInt -> CInt                  -- rows and cols
+    -> CInt -> CInt -> Ptr Double    -- coeffs
+    -> CInt -> CInt -> Ptr Double    -- bounds
+    -> CInt -> Ptr Double            -- result
+    -> IO CInt                       -- exit code
+
+simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
+simplexSparse m n c b = unsafePerformIO $ do
+    s <- createVector (2+n)
+    let fi = fromIntegral
+    app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
+    return s
+
+glpFR, glpLO, glpUP, glpDB, glpFX :: Double
+glpFR = 0
+glpLO = 1
+glpUP = 2
+glpDB = 3
+glpFX = 4
+
+{- Raw format of coeffs
+
+simplexDense
+
+((3+1) >< 3)
+  [ 10, 6, 4
+  ,  1, 1, 1
+  , 10, 4, 5
+  ,  2, 2, 6 :: Double]
+
+simplexSparse
+
+(12><3)
+ [ 0.0, 0.0, 10.0
+ , 0.0, 0.0,  6.0
+ , 0.0, 0.0,  4.0
+ , 1.0, 1.0,  1.0
+ , 1.0, 2.0,  1.0
+ , 1.0, 3.0,  1.0
+ , 2.0, 1.0, 10.0
+ , 2.0, 2.0,  4.0
+ , 2.0, 3.0,  5.0
+ , 3.0, 1.0,  2.0
+ , 3.0, 2.0,  2.0
+ , 3.0, 3.0,  6.0 ]
+  
+bounds = (6><3)
+  [ glpUP,0,100
+  , glpUP,0,600
+  , glpUP,0,300
+  , glpLO,0,0
+  , glpLO,0,0
+  , glpLO,0,0 ]
+  
+-}
+
diff --git a/packages/glpk/lib/Numeric/LinearProgramming/glpk.c b/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
new file mode 100644
index 0000000..0f4ad79
--- /dev/null
+++ b/packages/glpk/lib/Numeric/LinearProgramming/glpk.c
@@ -0,0 +1,147 @@
+typedef struct { double r, i; } doublecomplex;
+
+#define DVEC(A) int A##n, double*A##p
+#define CVEC(A) int A##n, doublecomplex*A##p
+#define DMAT(A) int A##r, int A##c, double*A##p
+#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
+
+#define AT(M,r,co) (M##p[(r)*M##c+(co)])
+
+/*-----------------------------------------------------*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <glpk.h>
+#include <math.h>
+
+int c_simplex_dense(DMAT(c), DMAT(b), DVEC(s)) {
+    glp_prob *lp;
+    lp = glp_create_prob();
+    glp_set_obj_dir(lp, GLP_MAX);
+
+    int m = cr-1;
+    int n = cc;
+    glp_add_rows(lp, m);
+    glp_add_cols(lp, n);
+    int * ia = malloc((1+m*n)*sizeof(int));
+    int * ja = malloc((1+m*n)*sizeof(int));
+    double * ar = malloc((1+m*n)*sizeof(double));
+    int i,j,k;
+    k=0;
+    for (i=1;i<=m;i++) {
+        for(j=1;j<=n;j++) {
+            k++;
+            ia[k] = i;
+            ja[k] = j;
+            ar[k] = AT(c,i,j-1);
+            //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
+        }
+    }
+    glp_load_matrix(lp, k, ia, ja, ar);
+    for (j=1;j<=n;j++) {
+        glp_set_obj_coef(lp, j, AT(c,0,j-1));
+    }
+    int t;
+    for (i=1;i<=m;i++) {
+    switch((int)rint(AT(b,i-1,0))) {
+        case 0: { t = GLP_FR; break; } 
+        case 1: { t = GLP_LO; break; }
+        case 2: { t = GLP_UP; break; }
+        case 3: { t = GLP_DB; break; }
+       default: { t = GLP_FX; break; }  
+    }
+    glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
+    }
+    for (j=1;j<=n;j++) {
+    switch((int)rint(AT(b,m+j-1,0))) {
+        case 0: { t = GLP_FR; break; } 
+        case 1: { t = GLP_LO; break; }
+        case 2: { t = GLP_UP; break; }
+        case 3: { t = GLP_DB; break; }
+       default: { t = GLP_FX; break; }  
+    }
+    glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
+    }
+    glp_term_out(0);
+    glp_simplex(lp, NULL);
+    sp[0] = glp_get_status(lp);
+    sp[1] = glp_get_obj_val(lp);
+    for (k=1; k<=n; k++) {
+        sp[k+1] = glp_get_col_prim(lp, k);
+    }
+    
+    glp_delete_prob(lp);
+    free(ia);
+    free(ja);
+    free(ar);
+
+    return 0;
+}
+
+/* ---------------------------------------------------- */
+
+int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {
+    glp_prob *lp;
+    lp = glp_create_prob();
+    glp_set_obj_dir(lp, GLP_MAX);
+    int i,j,k;
+    int tot = cr - n;
+    glp_add_rows(lp, m);
+    glp_add_cols(lp, n);
+
+    //printf("%d %d\n",m,n);
+
+    // the first n values
+    for (k=1;k<=n;k++) {
+        glp_set_obj_coef(lp, k, AT(c, k-1, 2));
+        //printf("%d %f\n",k,AT(c, k-1, 2));
+    }
+
+    int * ia = malloc((1+tot)*sizeof(int));
+    int * ja = malloc((1+tot)*sizeof(int));
+    double * ar = malloc((1+tot)*sizeof(double));
+
+    for (k=1; k<= tot; k++) {
+        ia[k] = rint(AT(c,k-1+n,0));
+        ja[k] = rint(AT(c,k-1+n,1));
+        ar[k] =      AT(c,k-1+n,2);
+        //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
+    }
+    glp_load_matrix(lp, tot, ia, ja, ar);
+    
+    int t;
+    for (i=1;i<=m;i++) {
+    switch((int)rint(AT(b,i-1,0))) {
+        case 0: { t = GLP_FR; break; } 
+        case 1: { t = GLP_LO; break; }
+        case 2: { t = GLP_UP; break; }
+        case 3: { t = GLP_DB; break; }
+       default: { t = GLP_FX; break; }  
+    }
+    glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
+    }
+    for (j=1;j<=n;j++) {
+    switch((int)rint(AT(b,m+j-1,0))) {
+        case 0: { t = GLP_FR; break; } 
+        case 1: { t = GLP_LO; break; }
+        case 2: { t = GLP_UP; break; }
+        case 3: { t = GLP_DB; break; }
+       default: { t = GLP_FX; break; }  
+    }
+    glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
+    }
+    glp_term_out(0);
+    glp_simplex(lp, NULL);
+    sp[0] = glp_get_status(lp);
+    sp[1] = glp_get_obj_val(lp);
+    for (k=1; k<=n; k++) {
+        sp[k+1] = glp_get_col_prim(lp, k);
+    }
+    glp_delete_prob(lp);
+    free(ia);
+    free(ja);
+    free(ar);
+
+    return 0;
+}
+