simple glpk interface

2 files changed, 358 insertions, 0 deletions

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;
+}
+