summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--packages/glpk/examples/simplex5.hs27
-rw-r--r--packages/glpk/hmatrix-glpk.cabal1
-rw-r--r--packages/glpk/src/C/glpk.c130
-rw-r--r--packages/glpk/src/Numeric/LinearProgramming.hs39
4 files changed, 132 insertions, 65 deletions
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 229197f..a9859f9 100644
--- a/packages/glpk/hmatrix-glpk.cabal
+++ b/packages/glpk/hmatrix-glpk.cabal
@@ -20,6 +20,7 @@ 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.16, containers >= 0.5.4.0 26 Build-Depends: base <5, hmatrix >= 0.16, containers >= 0.5.4.0
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 a54eb59..6a0c47d 100644
--- a/packages/glpk/src/Numeric/LinearProgramming.hs
+++ b/packages/glpk/src/Numeric/LinearProgramming.hs
@@ -75,6 +75,7 @@ Multiple bounds for a variable are not allowed, instead of
75 75
76module Numeric.LinearProgramming( 76module Numeric.LinearProgramming(
77 simplex, 77 simplex,
78 exact,
78 sparseOfGeneral, 79 sparseOfGeneral,
79 Optimization(..), 80 Optimization(..),
80 Constraints(..), 81 Constraints(..),
@@ -132,8 +133,8 @@ sparseOfGeneral :: Constraints -> Constraints
132sparseOfGeneral (General cs) = 133sparseOfGeneral (General cs) =
133 Sparse $ map (\bl -> 134 Sparse $ map (\bl ->
134 let cl = obj bl in 135 let cl = obj bl in
135 let m = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in 136 let cl_unique = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in
136 withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] m)) cs 137 withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] cl_unique)) cs
137sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected" 138sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected"
138 139
139simplex :: Optimization -> Constraints -> Bounds -> Solution 140simplex :: Optimization -> Constraints -> Bounds -> Solution
@@ -156,6 +157,27 @@ simplex opt (Sparse constr) bnds = extract sg sol where
156 157
157simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds 158simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds
158 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
159adapt :: Optimization -> (Int, Double, [Double]) 181adapt :: Optimization -> (Int, Double, [Double])
160adapt opt = case opt of 182adapt opt = case opt of
161 Maximize x -> (size x, 1 ,x) 183 Maximize x -> (size x, 1 ,x)
@@ -265,6 +287,19 @@ simplexSparse m n c b = unsafePerformIO $ do
265 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse" 287 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
266 return s 288 return s
267 289
290foreign import ccall unsafe "c_exact_sparse" c_exact_sparse
291 :: CInt -> CInt -- rows and cols
292 -> CInt -> CInt -> Ptr Double -- coeffs
293 -> CInt -> CInt -> Ptr Double -- bounds
294 -> CInt -> Ptr Double -- result
295 -> IO CInt -- exit code
296
297exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
298exactSparse m n c b = unsafePerformIO $ do
299 s <- createVector (2+n)
300 app3 (c_exact_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_exact_sparse"
301 return s
302
268glpFR, glpLO, glpUP, glpDB, glpFX :: Double 303glpFR, glpLO, glpUP, glpDB, glpFX :: Double
269glpFR = 0 304glpFR = 0
270glpLO = 1 305glpLO = 1
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 13:22:30 +0000