From f38b4a3076cfae023559ce61cb2a443c809b7a6f Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 21 Feb 2010 18:26:23 +0000 Subject: simple glpk interface --- packages/glpk/LICENSE | 2 + packages/glpk/Setup.lhs | 4 + packages/glpk/examples/simplex1.hs | 21 ++ packages/glpk/examples/simplex2.hs | 16 ++ packages/glpk/examples/simplex3.hs | 24 +++ packages/glpk/examples/simplex4.hs | 24 +++ packages/glpk/glpk.cabal | 35 ++++ packages/glpk/lib/Numeric/LinearProgramming.hs | 211 +++++++++++++++++++++ packages/glpk/lib/Numeric/LinearProgramming/glpk.c | 147 ++++++++++++++ 9 files changed, 484 insertions(+) create mode 100644 packages/glpk/LICENSE create mode 100644 packages/glpk/Setup.lhs create mode 100644 packages/glpk/examples/simplex1.hs create mode 100644 packages/glpk/examples/simplex2.hs create mode 100644 packages/glpk/examples/simplex3.hs create mode 100644 packages/glpk/examples/simplex4.hs create mode 100644 packages/glpk/glpk.cabal create mode 100644 packages/glpk/lib/Numeric/LinearProgramming.hs create mode 100644 packages/glpk/lib/Numeric/LinearProgramming/glpk.c 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 +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 +#include +#include +#include + +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; +} + -- cgit v1.2.3