summaryrefslogtreecommitdiff
path: root/packages/glpk
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2015-04-10 09:53:23 +0200
committerAlberto Ruiz <aruiz@um.es>2015-04-10 09:53:23 +0200
commite499467d2c67f214b0871376b068c5b6fc7896a1 (patch)
tree166021773cbb2c066811ee27830dd5618114ab10 /packages/glpk
parentdcc03a4a764cb8683b80758af97fcbcc9aadba73 (diff)
merge recent
Diffstat (limited to 'packages/glpk')
-rw-r--r--packages/glpk/examples/simplex1.hs6
-rw-r--r--packages/glpk/examples/simplex2.hs7
-rw-r--r--packages/glpk/examples/simplex3.hs2
-rw-r--r--packages/glpk/examples/simplex4.hs2
-rw-r--r--packages/glpk/examples/simplex5.hs27
-rw-r--r--packages/glpk/hmatrix-glpk.cabal3
-rw-r--r--packages/glpk/src/C/glpk.c130
-rw-r--r--packages/glpk/src/Numeric/LinearProgramming.hs75
8 files changed, 177 insertions, 75 deletions
diff --git a/packages/glpk/examples/simplex1.hs b/packages/glpk/examples/simplex1.hs
index e7aeaa9..a326555 100644
--- a/packages/glpk/examples/simplex1.hs
+++ b/packages/glpk/examples/simplex1.hs
@@ -9,9 +9,9 @@ constr = Dense [ [1,1,1] :<=: 100
9 , [2,2,6] :<=: 300 ] 9 , [2,2,6] :<=: 300 ]
10 10
11-- default bounds 11-- default bounds
12bnds = [ 1 :=>: 0 12bnds = [ 1 :>=: 0
13 , 2 :=>: 0 13 , 2 :>=: 0
14 , 3 :=>: 0 ] 14 , 3 :>=: 0 ]
15 15
16main = do 16main = do
17 print $ simplex objFun constr [] 17 print $ simplex objFun constr []
diff --git a/packages/glpk/examples/simplex2.hs b/packages/glpk/examples/simplex2.hs
index f4e27fd..0d83ca9 100644
--- a/packages/glpk/examples/simplex2.hs
+++ b/packages/glpk/examples/simplex2.hs
@@ -10,9 +10,14 @@ constr2 = Dense [ [2,1,0] :<=: 10
10 , [0,1,5] :<=: 20 10 , [0,1,5] :<=: 20
11 ] 11 ]
12 12
13constr3 = General [ [1#1, 1#1, 1#2] :<=: 10
14 , [1#2, 5#3] :<=: 20
15 ]
16
13main = do 17main = do
14 print $ simplex prob constr1 [] 18 print $ simplex prob constr1 []
15 print $ simplex prob constr2 [] 19 print $ simplex prob constr2 []
16 print $ simplex prob constr2 [ 2 :=>: 1, 3 :&: (2,7)] 20 print $ simplex prob constr3 []
21 print $ simplex prob constr2 [ 2 :>=: 1, 3 :&: (2,7)]
17 print $ simplex prob constr2 [ Free 2 ] 22 print $ simplex prob constr2 [ Free 2 ]
18 23
diff --git a/packages/glpk/examples/simplex3.hs b/packages/glpk/examples/simplex3.hs
index e093124..0997320 100644
--- a/packages/glpk/examples/simplex3.hs
+++ b/packages/glpk/examples/simplex3.hs
@@ -11,7 +11,7 @@ constr = Dense
11 , [0.03, 0.05, 0.08, 0.02, 0.06, 0.01, 0] :<=: 100 11 , [0.03, 0.05, 0.08, 0.02, 0.06, 0.01, 0] :<=: 100
12 , [0.02, 0.04, 0.01, 0.02, 0.02, 0, 0] :<=: 40 12 , [0.02, 0.04, 0.01, 0.02, 0.02, 0, 0] :<=: 40
13 , [0.02, 0.03, 0, 0, 0.01, 0, 0] :<=: 30 13 , [0.02, 0.03, 0, 0, 0.01, 0, 0] :<=: 30
14 , [0.7, 0.75, 0.8, 0.75, 0.8, 0.97, 0] :=>: 1500 14 , [0.7, 0.75, 0.8, 0.75, 0.8, 0.97, 0] :>=: 1500
15 , [0.02, 0.06, 0.08, 0.12, 0.02, 0.01, 0.97] :&: (250,300) 15 , [0.02, 0.06, 0.08, 0.12, 0.02, 0.01, 0.97] :&: (250,300)
16 ] 16 ]
17 17
diff --git a/packages/glpk/examples/simplex4.hs b/packages/glpk/examples/simplex4.hs
index 9a205ad..22b131c 100644
--- a/packages/glpk/examples/simplex4.hs
+++ b/packages/glpk/examples/simplex4.hs
@@ -11,7 +11,7 @@ constr = Sparse
11 , [0.03#1, 0.05#2, 0.08#3, 0.02#4, 0.06#5, 0.01#6] :<=: 100 11 , [0.03#1, 0.05#2, 0.08#3, 0.02#4, 0.06#5, 0.01#6] :<=: 100
12 , [0.02#1, 0.04#2, 0.01#3, 0.02#4, 0.02#5] :<=: 40 12 , [0.02#1, 0.04#2, 0.01#3, 0.02#4, 0.02#5] :<=: 40
13 , [0.02#1, 0.03#2, 0.01#5] :<=: 30 13 , [0.02#1, 0.03#2, 0.01#5] :<=: 30
14 , [0.7#1, 0.75#2, 0.8#3, 0.75#4, 0.8#5, 0.97#6] :=>: 1500 14 , [0.7#1, 0.75#2, 0.8#3, 0.75#4, 0.8#5, 0.97#6] :>=: 1500
15 , [0.02#1, 0.06#2, 0.08#3, 0.12#4, 0.02#5, 0.01#6, 0.97#7] :&: (250,300) 15 , [0.02#1, 0.06#2, 0.08#3, 0.12#4, 0.02#5, 0.01#6, 0.97#7] :&: (250,300)
16 ] 16 ]
17 17
diff --git a/packages/glpk/examples/simplex5.hs b/packages/glpk/examples/simplex5.hs
new file mode 100644
index 0000000..ecbcdaa
--- /dev/null
+++ b/packages/glpk/examples/simplex5.hs
@@ -0,0 +1,27 @@
1import Numeric.LinearProgramming
2
3-- This is a linear program from the paper "Picking vs. Guessing Secrets: A Game-theoretic Analysis"
4
5gamma = 100000 :: Double
6sigma = 1 :: Double
7n = 64 :: Int
8cost_fun :: Int -> Double
9cost_fun i = (fromIntegral i) / (fromIntegral n)
10size_fun :: Int -> Double
11size_fun i = 2^(fromIntegral i)
12
13prob = Minimize $ map cost_fun [1..n]
14bnds = [i :&: (0,1) | i <- [1..n]]
15
16constr1 = [[1 # i | i <- [1..n]] :==: 1] ++
17 [[1/(size_fun i) # i,
18 -1/(size_fun (i+1)) # i+1] :>=: 0 | i <- [1..n-1]] ++
19 [(
20 [gamma#i | i <- [1..k]] ++
21 (concat [[sigma*(size_fun i) # j | j <- [1..i-1]] | i <- [1..k]]) ++
22 [((size_fun i) - 1)/2 # i | i <- [1..k]])
23 :<=: (sigma * (foldr (+) 0 (map size_fun [1..k]))) | k <- [1..n]]
24
25main = do
26 print $ simplex prob (General constr1) bnds -- NoFeasible
27 print $ exact prob (General constr1) bnds -- solution found
diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal
index 4764bfe..8593e0a 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -20,9 +20,10 @@ extra-source-files: examples/simplex1.hs
20 examples/simplex2.hs 20 examples/simplex2.hs
21 examples/simplex3.hs 21 examples/simplex3.hs
22 examples/simplex4.hs 22 examples/simplex4.hs
23 examples/simplex5.hs
23 24
24library 25library
25 Build-Depends: base <5, hmatrix >= 0.17 26 Build-Depends: base <5, hmatrix >= 0.17, containers
26 27
27 hs-source-dirs: src 28 hs-source-dirs: src
28 29
diff --git a/packages/glpk/src/C/glpk.c b/packages/glpk/src/C/glpk.c
index bfbb435..86b1277 100644
--- a/packages/glpk/src/C/glpk.c
+++ b/packages/glpk/src/C/glpk.c
@@ -10,67 +10,71 @@
10 10
11/*-----------------------------------------------------*/ 11/*-----------------------------------------------------*/
12 12
13int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) { 13#define C_X_SPARSE(X) \
14 glp_prob *lp; 14 int c_##X##_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) { \
15 lp = glp_create_prob(); 15 glp_prob *lp; \
16 glp_set_obj_dir(lp, GLP_MAX); 16 lp = glp_create_prob(); \
17 int i,j,k; 17 glp_set_obj_dir(lp, GLP_MAX); \
18 int tot = cr - n; 18 int i,j,k; \
19 glp_add_rows(lp, m); 19 int tot = cr - n; \
20 glp_add_cols(lp, n); 20 glp_add_rows(lp, m); \
21 glp_add_cols(lp, n); \
22 \
23 /*printf("%d %d\n",m,n);*/ \
24 \
25 /* the first n values */ \
26 for (k=1;k<=n;k++) { \
27 glp_set_obj_coef(lp, k, AT(c, k-1, 2)); \
28 /*printf("%d %f\n",k,AT(c, k-1, 2)); */ \
29 } \
30 \
31 int * ia = malloc((1+tot)*sizeof(int)); \
32 int * ja = malloc((1+tot)*sizeof(int)); \
33 double * ar = malloc((1+tot)*sizeof(double)); \
34 \
35 for (k=1; k<= tot; k++) { \
36 ia[k] = rint(AT(c,k-1+n,0)); \
37 ja[k] = rint(AT(c,k-1+n,1)); \
38 ar[k] = AT(c,k-1+n,2); \
39 /*printf("%d %d %f\n",ia[k],ja[k],ar[k]);*/ \
40 } \
41 glp_load_matrix(lp, tot, ia, ja, ar); \
42 \
43 int t; \
44 for (i=1;i<=m;i++) { \
45 switch((int)rint(AT(b,i-1,0))) { \
46 case 0: { t = GLP_FR; break; } \
47 case 1: { t = GLP_LO; break; } \
48 case 2: { t = GLP_UP; break; } \
49 case 3: { t = GLP_DB; break; } \
50 default: { t = GLP_FX; break; } \
51 } \
52 glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2)); \
53 } \
54 for (j=1;j<=n;j++) { \
55 switch((int)rint(AT(b,m+j-1,0))) { \
56 case 0: { t = GLP_FR; break; } \
57 case 1: { t = GLP_LO; break; } \
58 case 2: { t = GLP_UP; break; } \
59 case 3: { t = GLP_DB; break; } \
60 default: { t = GLP_FX; break; } \
61 } \
62 glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2)); \
63 } \
64 glp_term_out(0); \
65 glp_##X(lp, NULL); \
66 sp[0] = glp_get_status(lp); \
67 sp[1] = glp_get_obj_val(lp); \
68 for (k=1; k<=n; k++) { \
69 sp[k+1] = glp_get_col_prim(lp, k); \
70 } \
71 glp_delete_prob(lp); \
72 free(ia); \
73 free(ja); \
74 free(ar); \
75 \
76 return 0; \
77 } \
21 78
22 //printf("%d %d\n",m,n); 79C_X_SPARSE(simplex);
23 80C_X_SPARSE(exact);
24 // the first n values
25 for (k=1;k<=n;k++) {
26 glp_set_obj_coef(lp, k, AT(c, k-1, 2));
27 //printf("%d %f\n",k,AT(c, k-1, 2));
28 }
29
30 int * ia = malloc((1+tot)*sizeof(int));
31 int * ja = malloc((1+tot)*sizeof(int));
32 double * ar = malloc((1+tot)*sizeof(double));
33
34 for (k=1; k<= tot; k++) {
35 ia[k] = rint(AT(c,k-1+n,0));
36 ja[k] = rint(AT(c,k-1+n,1));
37 ar[k] = AT(c,k-1+n,2);
38 //printf("%d %d %f\n",ia[k],ja[k],ar[k]);
39 }
40 glp_load_matrix(lp, tot, ia, ja, ar);
41
42 int t;
43 for (i=1;i<=m;i++) {
44 switch((int)rint(AT(b,i-1,0))) {
45 case 0: { t = GLP_FR; break; }
46 case 1: { t = GLP_LO; break; }
47 case 2: { t = GLP_UP; break; }
48 case 3: { t = GLP_DB; break; }
49 default: { t = GLP_FX; break; }
50 }
51 glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));
52 }
53 for (j=1;j<=n;j++) {
54 switch((int)rint(AT(b,m+j-1,0))) {
55 case 0: { t = GLP_FR; break; }
56 case 1: { t = GLP_LO; break; }
57 case 2: { t = GLP_UP; break; }
58 case 3: { t = GLP_DB; break; }
59 default: { t = GLP_FX; break; }
60 }
61 glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));
62 }
63 glp_term_out(0);
64 glp_simplex(lp, NULL);
65 sp[0] = glp_get_status(lp);
66 sp[1] = glp_get_obj_val(lp);
67 for (k=1; k<=n; k++) {
68 sp[k+1] = glp_get_col_prim(lp, k);
69 }
70 glp_delete_prob(lp);
71 free(ia);
72 free(ja);
73 free(ar);
74
75 return 0;
76}
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs
index d2e9f3c..7bf4279 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
49 ] 49 ]
50@ 50@
51 51
52Note 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.
53
54@
55constr3 = General [ [1\#1, 1\#1, 1\#2] :<=: 10
56 , [1\#2, 5\#3] :<=: 20
57 ]
58@
59
52By default all variables are bounded as @x_i >= 0@, but this can be 60By default all variables are bounded as @x_i >= 0@, but this can be
53changed: 61changed:
54 62
@@ -67,6 +75,8 @@ Multiple bounds for a variable are not allowed, instead of
67 75
68module Numeric.LinearProgramming( 76module Numeric.LinearProgramming(
69 simplex, 77 simplex,
78 exact,
79 sparseOfGeneral,
70 Optimization(..), 80 Optimization(..),
71 Constraints(..), 81 Constraints(..),
72 Bounds, 82 Bounds,
@@ -82,13 +92,14 @@ import System.IO.Unsafe(unsafePerformIO)
82import Foreign.C.Types 92import Foreign.C.Types
83import Data.List((\\),sortBy,nub) 93import Data.List((\\),sortBy,nub)
84import Data.Function(on) 94import Data.Function(on)
95import qualified Data.Map.Strict as Map
85 96
86--import Debug.Trace 97--import Debug.Trace
87--debug x = trace (show x) x 98--debug x = trace (show x) x
88 99
89----------------------------------------------------- 100-----------------------------------------------------
90 101
91-- | Coefficient of a variable for a sparse representation of constraints. 102-- | Coefficient of a variable for a sparse and general representations of constraints.
92(#) :: Double -> Int -> (Double,Int) 103(#) :: Double -> Int -> (Double,Int)
93infixl 5 # 104infixl 5 #
94(#) = (,) 105(#) = (,)
@@ -108,18 +119,29 @@ data Solution = Undefined
108 | Unbounded 119 | Unbounded
109 deriving Show 120 deriving Show
110 121
111data Constraints = Dense [ Bound [Double] ] 122data Constraints = Dense [ Bound [Double] ]
112 | Sparse [ Bound [(Double,Int)] ] 123 | Sparse [ Bound [(Double,Int)] ]
124 | General [ Bound [(Double,Int)] ]
113 125
114data Optimization = Maximize [Double] 126data Optimization = Maximize [Double]
115 | Minimize [Double] 127 | Minimize [Double]
116 128
117type Bounds = [Bound Int] 129type Bounds = [Bound Int]
118 130
131-- | Convert a system of General constraints to one with unique coefficients.
132sparseOfGeneral :: Constraints -> Constraints
133sparseOfGeneral (General cs) =
134 Sparse $ map (\bl ->
135 let cl = obj bl in
136 let cl_unique = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in
137 withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] cl_unique)) cs
138sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected"
139
119simplex :: Optimization -> Constraints -> Bounds -> Solution 140simplex :: Optimization -> Constraints -> Bounds -> Solution
120 141
121simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds 142simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds
122simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds 143simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
144simplex opt (General []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
123 145
124simplex opt (Dense constr) bnds = extract sg sol where 146simplex opt (Dense constr) bnds = extract sg sol where
125 sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds) 147 sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
@@ -133,6 +155,29 @@ simplex opt (Sparse constr) bnds = extract sg sol where
133 m = length constr 155 m = length constr
134 (sz, sg, objfun) = adapt opt 156 (sz, sg, objfun) = adapt opt
135 157
158simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds
159
160-- | Simplex method with exact internal arithmetic. See glp_exact in glpk documentation for more information.
161exact :: Optimization -> Constraints -> Bounds -> Solution
162
163exact opt (Dense []) bnds = exact opt (Sparse []) bnds
164exact opt (Sparse []) bnds = exact opt (Sparse [Free [0#1]]) bnds
165exact opt (General []) bnds = exact opt (Sparse [Free [0#1]]) bnds
166
167exact opt (Dense constr) bnds = extract sg sol where
168 sol = exactSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
169 n = length objfun
170 m = length constr
171 (sz, sg, objfun) = adapt opt
172
173exact opt (Sparse constr) bnds = extract sg sol where
174 sol = exactSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)
175 n = length objfun
176 m = length constr
177 (sz, sg, objfun) = adapt opt
178
179exact opt constr@(General _) bnds = exact opt (sparseOfGeneral constr) bnds
180
136adapt :: Optimization -> (Int, Double, [Double]) 181adapt :: Optimization -> (Int, Double, [Double])
137adapt opt = case opt of 182adapt opt = case opt of
138 Maximize x -> (sz x, 1 ,x) 183 Maximize x -> (sz x, 1 ,x)
@@ -163,6 +208,13 @@ obj (x :&: _) = x
163obj (x :==: _) = x 208obj (x :==: _) = x
164obj (Free x) = x 209obj (Free x) = x
165 210
211withObj :: Bound t -> t -> Bound t
212withObj (_ :<=: b) x = (x :<=: b)
213withObj (_ :>=: b) x = (x :>=: b)
214withObj (_ :&: b) x = (x :&: b)
215withObj (_ :==: b) x = (x :==: b)
216withObj (Free _) x = Free x
217
166tb :: Bound t -> Double 218tb :: Bound t -> Double
167tb (_ :<=: _) = glpUP 219tb (_ :<=: _) = glpUP
168tb (_ :>=: _) = glpLO 220tb (_ :>=: _) = glpLO
@@ -236,6 +288,19 @@ simplexSparse m n c b = unsafePerformIO $ do
236 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse" 288 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
237 return s 289 return s
238 290
291foreign import ccall unsafe "c_exact_sparse" c_exact_sparse
292 :: CInt -> CInt -- rows and cols
293 -> CInt -> CInt -> Ptr Double -- coeffs
294 -> CInt -> CInt -> Ptr Double -- bounds
295 -> CInt -> Ptr Double -- result
296 -> IO CInt -- exit code
297
298exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
299exactSparse m n c b = unsafePerformIO $ do
300 s <- createVector (2+n)
301 app3 (c_exact_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_exact_sparse"
302 return s
303
239glpFR, glpLO, glpUP, glpDB, glpFX :: Double 304glpFR, glpLO, glpUP, glpDB, glpFX :: Double
240glpFR = 0 305glpFR = 0
241glpLO = 1 306glpLO = 1
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-18 12:27:54 +0000