diff options
Diffstat (limited to 'packages')
-rw-r--r-- | packages/base/THANKS.md | 15 | ||||
-rw-r--r-- | packages/base/hmatrix.cabal | 1 | ||||
-rw-r--r-- | packages/base/src/C/vector-aux.c | 1 | ||||
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Numeric.hs | 2 | ||||
-rw-r--r-- | packages/base/src/Numeric/Chain.hs | 2 | ||||
-rw-r--r-- | packages/base/src/Numeric/LinearAlgebra/Static.hs | 30 | ||||
-rw-r--r-- | packages/glpk/examples/simplex1.hs | 6 | ||||
-rw-r--r-- | packages/glpk/examples/simplex2.hs | 7 | ||||
-rw-r--r-- | packages/glpk/examples/simplex3.hs | 2 | ||||
-rw-r--r-- | packages/glpk/examples/simplex4.hs | 2 | ||||
-rw-r--r-- | packages/glpk/examples/simplex5.hs | 27 | ||||
-rw-r--r-- | packages/glpk/hmatrix-glpk.cabal | 1 | ||||
-rw-r--r-- | packages/glpk/src/C/glpk.c | 130 | ||||
-rw-r--r-- | packages/glpk/src/Numeric/LinearProgramming.hs | 75 |
14 files changed, 200 insertions, 101 deletions
diff --git a/packages/base/THANKS.md b/packages/base/THANKS.md index fdbbe14..a4188eb 100644 --- a/packages/base/THANKS.md +++ b/packages/base/THANKS.md | |||
@@ -92,7 +92,7 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
92 | 92 | ||
93 | - Carter Schonwald helped with the configuration for Homebrew OS X and | 93 | - Carter Schonwald helped with the configuration for Homebrew OS X and |
94 | found a tolerance problem in test "1E5 rots". He also discovered | 94 | found a tolerance problem in test "1E5 rots". He also discovered |
95 | a bug in the signature of cmap. | 95 | a bug in the signature of cmap and fixed the cabal file. |
96 | 96 | ||
97 | - Duncan Coutts reported a problem with configure.hs and contributed | 97 | - Duncan Coutts reported a problem with configure.hs and contributed |
98 | a solution and a simplified Setup.lhs. | 98 | a solution and a simplified Setup.lhs. |
@@ -159,7 +159,8 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
159 | 159 | ||
160 | - Denis Laxalde separated the gsl tests from the base ones. | 160 | - Denis Laxalde separated the gsl tests from the base ones. |
161 | 161 | ||
162 | - Dominic Steinitz (idontgetoutmuch) reported a bug in the static diagonal creation functions. | 162 | - Dominic Steinitz (idontgetoutmuch) reported a bug in the static diagonal creation functions and |
163 | added Cholesky to Static. | ||
163 | 164 | ||
164 | - Dylan Thurston reported an error in the glpk documentation and ambiguity in | 165 | - Dylan Thurston reported an error in the glpk documentation and ambiguity in |
165 | the description of linearSolve. | 166 | the description of linearSolve. |
@@ -179,3 +180,13 @@ module reorganization, monadic mapVectorM, and many other improvements. | |||
179 | 180 | ||
180 | - Kiwamu Ishikura improved randomVector for OSX | 181 | - Kiwamu Ishikura improved randomVector for OSX |
181 | 182 | ||
183 | - C.J. East fixed the examples for simplex. | ||
184 | |||
185 | - Ben Gamari contributed fixes for ghc 7.10 | ||
186 | |||
187 | - Piotr Mardziel added general sparse constraints to simplex and the interface to glp_exact | ||
188 | |||
189 | - Maxim Baz fixed an instance declaration for ghc 7.11 | ||
190 | |||
191 | - Thomas M. DuBuisson fixed a C include file. | ||
192 | |||
diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index f5384f3..350023f 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal | |||
@@ -35,6 +35,7 @@ extra-source-files: src/C/lapack-aux.h | |||
35 | flag openblas | 35 | flag openblas |
36 | description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. | 36 | description: Link with OpenBLAS (https://github.com/xianyi/OpenBLAS) optimized libraries. |
37 | default: False | 37 | default: False |
38 | manual: True | ||
38 | 39 | ||
39 | library | 40 | library |
40 | 41 | ||
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 dc5724d..7a4dd29 100644 --- a/packages/base/src/Data/Packed/Internal/Numeric.hs +++ b/packages/base/src/Data/Packed/Internal/Numeric.hs | |||
@@ -240,7 +240,7 @@ instance Container Vector (Complex Float) | |||
240 | 240 | ||
241 | --------------------------------------------------------------- | 241 | --------------------------------------------------------------- |
242 | 242 | ||
243 | instance (Container Vector a) => Container Matrix a | 243 | instance (Fractional a, Element a, Container Vector a) => Container Matrix a |
244 | where | 244 | where |
245 | size' = size | 245 | size' = size |
246 | scale' x = liftMatrix (scale' x) | 246 | scale' x = liftMatrix (scale' x) |
diff --git a/packages/base/src/Numeric/Chain.hs b/packages/base/src/Numeric/Chain.hs index 55e2df6..64c09c0 100644 --- a/packages/base/src/Numeric/Chain.hs +++ b/packages/base/src/Numeric/Chain.hs | |||
@@ -14,6 +14,8 @@ | |||
14 | -- | 14 | -- |
15 | ----------------------------------------------------------------------------- | 15 | ----------------------------------------------------------------------------- |
16 | 16 | ||
17 | {-# LANGUAGE FlexibleContexts #-} | ||
18 | |||
17 | module Numeric.Chain ( | 19 | module Numeric.Chain ( |
18 | optimiseMult, | 20 | optimiseMult, |
19 | ) where | 21 | ) where |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index a26cc4c..4c3186f 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs | |||
@@ -1,5 +1,3 @@ | |||
1 | #if __GLASGOW_HASKELL__ >= 708 | ||
2 | |||
3 | {-# LANGUAGE DataKinds #-} | 1 | {-# LANGUAGE DataKinds #-} |
4 | {-# LANGUAGE KindSignatures #-} | 2 | {-# LANGUAGE KindSignatures #-} |
5 | {-# LANGUAGE GeneralizedNewtypeDeriving #-} | 3 | {-# LANGUAGE GeneralizedNewtypeDeriving #-} |
@@ -51,7 +49,7 @@ module Numeric.LinearAlgebra.Static( | |||
51 | linSolve, (<\>), | 49 | linSolve, (<\>), |
52 | -- * Factorizations | 50 | -- * Factorizations |
53 | svd, withCompactSVD, svdTall, svdFlat, Eigen(..), | 51 | svd, withCompactSVD, svdTall, svdFlat, Eigen(..), |
54 | withNullspace, qr, | 52 | withNullspace, qr, chol, |
55 | -- * Misc | 53 | -- * Misc |
56 | mean, | 54 | mean, |
57 | Disp(..), Domain(..), | 55 | Disp(..), Domain(..), |
@@ -67,7 +65,7 @@ import Numeric.LinearAlgebra.HMatrix hiding ( | |||
67 | row,col,vector,matrix,linspace,toRows,toColumns, | 65 | row,col,vector,matrix,linspace,toRows,toColumns, |
68 | (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', | 66 | (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', |
69 | eigenvalues,eigenvaluesSH,eigenvaluesSH',build, | 67 | eigenvalues,eigenvaluesSH,eigenvaluesSH',build, |
70 | qr,size,app,mul,dot) | 68 | qr,size,app,mul,dot,chol) |
71 | import qualified Numeric.LinearAlgebra.HMatrix as LA | 69 | import qualified Numeric.LinearAlgebra.HMatrix as LA |
72 | import Data.Proxy(Proxy) | 70 | import Data.Proxy(Proxy) |
73 | import Numeric.LinearAlgebra.Static.Internal | 71 | import Numeric.LinearAlgebra.Static.Internal |
@@ -183,6 +181,7 @@ a ¦ b = tr (tr a —— tr b) | |||
183 | type Sq n = L n n | 181 | type Sq n = L n n |
184 | --type CSq n = CL n n | 182 | --type CSq n = CL n n |
185 | 183 | ||
184 | |||
186 | type GL = forall n m . (KnownNat n, KnownNat m) => L m n | 185 | type GL = forall n m . (KnownNat n, KnownNat m) => L m n |
187 | type GSq = forall n . KnownNat n => Sq n | 186 | type GSq = forall n . KnownNat n => Sq n |
188 | 187 | ||
@@ -305,6 +304,9 @@ instance KnownNat n => Eigen (Sq n) (C n) (M n n) | |||
305 | where | 304 | where |
306 | (l,v) = LA.eig m | 305 | (l,v) = LA.eig m |
307 | 306 | ||
307 | chol :: KnownNat n => Sym n -> Sq n | ||
308 | chol (extract . unSym -> m) = mkL $ LA.cholSH m | ||
309 | |||
308 | -------------------------------------------------------------------------------- | 310 | -------------------------------------------------------------------------------- |
309 | 311 | ||
310 | withNullspace | 312 | withNullspace |
@@ -614,23 +616,3 @@ instance (KnownNat n', KnownNat m') => Testable (L n' m') | |||
614 | where | 616 | where |
615 | checkT _ = test | 617 | checkT _ = test |
616 | 618 | ||
617 | #else | ||
618 | |||
619 | {- | | ||
620 | Module : Numeric.LinearAlgebra.Static | ||
621 | Copyright : (c) Alberto Ruiz 2014 | ||
622 | License : BSD3 | ||
623 | Stability : experimental | ||
624 | |||
625 | Experimental interface with statically checked dimensions. | ||
626 | |||
627 | This module requires GHC >= 7.8 | ||
628 | |||
629 | -} | ||
630 | |||
631 | module Numeric.LinearAlgebra.Static | ||
632 | {-# WARNING "This module requires GHC >= 7.8" #-} | ||
633 | where | ||
634 | |||
635 | #endif | ||
636 | |||
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 |
12 | bnds = [ 1 :=>: 0 | 12 | bnds = [ 1 :>=: 0 |
13 | , 2 :=>: 0 | 13 | , 2 :>=: 0 |
14 | , 3 :=>: 0 ] | 14 | , 3 :>=: 0 ] |
15 | 15 | ||
16 | main = do | 16 | main = 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 | ||
13 | constr3 = General [ [1#1, 1#1, 1#2] :<=: 10 | ||
14 | , [1#2, 5#3] :<=: 20 | ||
15 | ] | ||
16 | |||
13 | main = do | 17 | main = 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 @@ | |||
1 | import Numeric.LinearProgramming | ||
2 | |||
3 | -- This is a linear program from the paper "Picking vs. Guessing Secrets: A Game-theoretic Analysis" | ||
4 | |||
5 | gamma = 100000 :: Double | ||
6 | sigma = 1 :: Double | ||
7 | n = 64 :: Int | ||
8 | cost_fun :: Int -> Double | ||
9 | cost_fun i = (fromIntegral i) / (fromIntegral n) | ||
10 | size_fun :: Int -> Double | ||
11 | size_fun i = 2^(fromIntegral i) | ||
12 | |||
13 | prob = Minimize $ map cost_fun [1..n] | ||
14 | bnds = [i :&: (0,1) | i <- [1..n]] | ||
15 | |||
16 | constr1 = [[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 | |||
25 | main = 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..e64e3e8 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 | ||
24 | library | 25 | library |
25 | Build-Depends: base <5, hmatrix >= 0.17 | 26 | Build-Depends: base <5, hmatrix >= 0.17 |
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 | ||
13 | int 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); | 79 | C_X_SPARSE(simplex); |
23 | 80 | C_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 | ||
52 | Note 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 | @ | ||
55 | constr3 = General [ [1\#1, 1\#1, 1\#2] :<=: 10 | ||
56 | , [1\#2, 5\#3] :<=: 20 | ||
57 | ] | ||
58 | @ | ||
59 | |||
52 | By default all variables are bounded as @x_i >= 0@, but this can be | 60 | By default all variables are bounded as @x_i >= 0@, but this can be |
53 | changed: | 61 | changed: |
54 | 62 | ||
@@ -67,6 +75,8 @@ Multiple bounds for a variable are not allowed, instead of | |||
67 | 75 | ||
68 | module Numeric.LinearProgramming( | 76 | module 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) | |||
82 | import Foreign.C.Types | 92 | import Foreign.C.Types |
83 | import Data.List((\\),sortBy,nub) | 93 | import Data.List((\\),sortBy,nub) |
84 | import Data.Function(on) | 94 | import Data.Function(on) |
95 | import 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) |
93 | infixl 5 # | 104 | infixl 5 # |
94 | (#) = (,) | 105 | (#) = (,) |
@@ -108,18 +119,29 @@ data Solution = Undefined | |||
108 | | Unbounded | 119 | | Unbounded |
109 | deriving Show | 120 | deriving Show |
110 | 121 | ||
111 | data Constraints = Dense [ Bound [Double] ] | 122 | data Constraints = Dense [ Bound [Double] ] |
112 | | Sparse [ Bound [(Double,Int)] ] | 123 | | Sparse [ Bound [(Double,Int)] ] |
124 | | General [ Bound [(Double,Int)] ] | ||
113 | 125 | ||
114 | data Optimization = Maximize [Double] | 126 | data Optimization = Maximize [Double] |
115 | | Minimize [Double] | 127 | | Minimize [Double] |
116 | 128 | ||
117 | type Bounds = [Bound Int] | 129 | type Bounds = [Bound Int] |
118 | 130 | ||
131 | -- | Convert a system of General constraints to one with unique coefficients. | ||
132 | sparseOfGeneral :: Constraints -> Constraints | ||
133 | sparseOfGeneral (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 | ||
138 | sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected" | ||
139 | |||
119 | simplex :: Optimization -> Constraints -> Bounds -> Solution | 140 | simplex :: Optimization -> Constraints -> Bounds -> Solution |
120 | 141 | ||
121 | simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds | 142 | simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds |
122 | simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds | 143 | simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds |
144 | simplex opt (General []) bnds = simplex opt (Sparse [Free [0#1]]) bnds | ||
123 | 145 | ||
124 | simplex opt (Dense constr) bnds = extract sg sol where | 146 | simplex 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 | ||
158 | simplex 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. | ||
161 | exact :: Optimization -> Constraints -> Bounds -> Solution | ||
162 | |||
163 | exact opt (Dense []) bnds = exact opt (Sparse []) bnds | ||
164 | exact opt (Sparse []) bnds = exact opt (Sparse [Free [0#1]]) bnds | ||
165 | exact opt (General []) bnds = exact opt (Sparse [Free [0#1]]) bnds | ||
166 | |||
167 | exact 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 | |||
173 | exact 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 | |||
179 | exact opt constr@(General _) bnds = exact opt (sparseOfGeneral constr) bnds | ||
180 | |||
136 | adapt :: Optimization -> (Int, Double, [Double]) | 181 | adapt :: Optimization -> (Int, Double, [Double]) |
137 | adapt opt = case opt of | 182 | adapt 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 | |||
163 | obj (x :==: _) = x | 208 | obj (x :==: _) = x |
164 | obj (Free x) = x | 209 | obj (Free x) = x |
165 | 210 | ||
211 | withObj :: Bound t -> t -> Bound t | ||
212 | withObj (_ :<=: b) x = (x :<=: b) | ||
213 | withObj (_ :>=: b) x = (x :>=: b) | ||
214 | withObj (_ :&: b) x = (x :&: b) | ||
215 | withObj (_ :==: b) x = (x :==: b) | ||
216 | withObj (Free _) x = Free x | ||
217 | |||
166 | tb :: Bound t -> Double | 218 | tb :: Bound t -> Double |
167 | tb (_ :<=: _) = glpUP | 219 | tb (_ :<=: _) = glpUP |
168 | tb (_ :>=: _) = glpLO | 220 | tb (_ :>=: _) = 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 | ||
291 | foreign 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 | |||
298 | exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double | ||
299 | exactSparse 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 | |||
239 | glpFR, glpLO, glpUP, glpDB, glpFX :: Double | 304 | glpFR, glpLO, glpUP, glpDB, glpFX :: Double |
240 | glpFR = 0 | 305 | glpFR = 0 |
241 | glpLO = 1 | 306 | glpLO = 1 |