summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2015-03-16 08:47:46 +0100
committerAlberto Ruiz <aruiz@um.es>2015-03-16 08:47:46 +0100
commitc8bf3420beba4d3ec6c024a7a068166f8b5a064f (patch)
treebacf418fa4b25306c53c5b3beca99f3fffc2e013
parent5a6bf1299d3bdad5289f559192d60fdaf33893e7 (diff)
parentc1081dc0641de9ef9fe37a45ae14452dcb8c0722 (diff)
Merge branch 'master' of github.com:albertoruiz/hmatrix
-rw-r--r--packages/base/src/C/vector-aux.c1
-rw-r--r--packages/base/src/Data/Packed/Internal/Numeric.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Static.hs7
-rw-r--r--packages/glpk/examples/simplex2.hs5
-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, 178 insertions, 72 deletions
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index dda47cb..599f69e 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -13,6 +13,7 @@ typedef float complex TCF;
13#include <math.h> 13#include <math.h>
14#include <stdio.h> 14#include <stdio.h>
15#include <stdlib.h> 15#include <stdlib.h>
16#include <stdint.h>
16 17
17#define MACRO(B) do {B} while (0) 18#define MACRO(B) do {B} while (0)
18#define ERROR(CODE) MACRO(return CODE;) 19#define ERROR(CODE) MACRO(return CODE;)
diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs
index 9adc023..257ad73 100644
--- a/packages/base/src/Data/Packed/Internal/Numeric.hs
+++ b/packages/base/src/Data/Packed/Internal/Numeric.hs
@@ -241,7 +241,7 @@ instance Container Vector (Complex Float)
241 241
242--------------------------------------------------------------- 242---------------------------------------------------------------
243 243
244instance (Container Vector a) => Container Matrix a 244instance (Fractional a, Element a, Container Vector a) => Container Matrix a
245 where 245 where
246 size' = size 246 size' = size
247 scale' x = liftMatrix (scale' x) 247 scale' x = liftMatrix (scale' x)
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs
index 037396d..3398e6a 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Static.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs
@@ -52,7 +52,7 @@ module Numeric.LinearAlgebra.Static(
52 linSolve, (<\>), 52 linSolve, (<\>),
53 -- * Factorizations 53 -- * Factorizations
54 svd, withCompactSVD, svdTall, svdFlat, Eigen(..), 54 svd, withCompactSVD, svdTall, svdFlat, Eigen(..),
55 withNullspace, qr, 55 withNullspace, qr, chol,
56 -- * Misc 56 -- * Misc
57 mean, 57 mean,
58 Disp(..), Domain(..), 58 Disp(..), Domain(..),
@@ -68,7 +68,7 @@ import Numeric.LinearAlgebra.HMatrix hiding (
68 row,col,vector,matrix,linspace,toRows,toColumns, 68 row,col,vector,matrix,linspace,toRows,toColumns,
69 (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', 69 (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',
70 eigenvalues,eigenvaluesSH,eigenvaluesSH',build, 70 eigenvalues,eigenvaluesSH,eigenvaluesSH',build,
71 qr,size,app,mul,dot) 71 qr,size,app,mul,dot,chol)
72import qualified Numeric.LinearAlgebra.HMatrix as LA 72import qualified Numeric.LinearAlgebra.HMatrix as LA
73import Data.Proxy(Proxy) 73import Data.Proxy(Proxy)
74import Numeric.LinearAlgebra.Static.Internal 74import Numeric.LinearAlgebra.Static.Internal
@@ -306,6 +306,9 @@ instance KnownNat n => Eigen (Sq n) (C n) (M n n)
306 where 306 where
307 (l,v) = LA.eig m 307 (l,v) = LA.eig m
308 308
309chol :: KnownNat n => Sym n -> Sq n
310chol (extract . unSym -> m) = mkL $ LA.cholSH m
311
309-------------------------------------------------------------------------------- 312--------------------------------------------------------------------------------
310 313
311withNullspace 314withNullspace
diff --git a/packages/glpk/examples/simplex2.hs b/packages/glpk/examples/simplex2.hs
index e9e8859..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 []
20 print $ simplex prob constr3 []
16 print $ simplex prob constr2 [ 2 :>=: 1, 3 :&: (2,7)] 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/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 cd761e0..a9859f9 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.16 26 Build-Depends: base <5, hmatrix >= 0.16, containers >= 0.5.4.0
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 b0537cc..6a0c47d 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 -> (size x, 1 ,x) 183 Maximize x -> (size x, 1 ,x)
@@ -162,6 +207,13 @@ obj (x :&: _) = x
162obj (x :==: _) = x 207obj (x :==: _) = x
163obj (Free x) = x 208obj (Free x) = x
164 209
210withObj :: Bound t -> t -> Bound t
211withObj (_ :<=: b) x = (x :<=: b)
212withObj (_ :>=: b) x = (x :>=: b)
213withObj (_ :&: b) x = (x :&: b)
214withObj (_ :==: b) x = (x :==: b)
215withObj (Free _) x = Free x
216
165tb :: Bound t -> Double 217tb :: Bound t -> Double
166tb (_ :<=: _) = glpUP 218tb (_ :<=: _) = glpUP
167tb (_ :>=: _) = glpLO 219tb (_ :>=: _) = glpLO
@@ -235,6 +287,19 @@ simplexSparse m n c b = unsafePerformIO $ do
235 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"
236 return s 288 return s
237 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
238glpFR, glpLO, glpUP, glpDB, glpFX :: Double 303glpFR, glpLO, glpUP, glpDB, glpFX :: Double
239glpFR = 0 304glpFR = 0
240glpLO = 1 305glpLO = 1