summaryrefslogtreecommitdiff
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
parentdcc03a4a764cb8683b80758af97fcbcc9aadba73 (diff)
merge recent
Diffstat
-rw-r--r--examples/inplace.hs39
-rw-r--r--packages/base/THANKS.md15
-rw-r--r--packages/base/hmatrix.cabal1
-rw-r--r--packages/base/src/C/vector-aux.c1
-rw-r--r--packages/base/src/Data/Packed/IO.hs1
-rw-r--r--packages/base/src/Data/Packed/Internal/Numeric.hs2
-rw-r--r--packages/base/src/Numeric/Chain.hs2
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Static.hs30
-rw-r--r--packages/base/src/Numeric/Vectorized.hs1
-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
17 files changed, 221 insertions, 123 deletions
diff --git a/examples/inplace.hs b/examples/inplace.hs
index dcfff56..574aa44 100644
--- a/examples/inplace.hs
+++ b/examples/inplace.hs
@@ -1,9 +1,8 @@
1-- some tests of the interface for pure 1-- some tests of the interface for pure
2-- computations with inplace updates 2-- computations with inplace updates
3 3
4import Numeric.LinearAlgebra 4import Numeric.LinearAlgebra.HMatrix
5import Data.Packed.ST 5import Numeric.LinearAlgebra.Devel
6import Data.Packed.Convert
7 6
8import Data.Array.Unboxed 7import Data.Array.Unboxed
9import Data.Array.ST 8import Data.Array.ST
@@ -15,15 +14,13 @@ main = sequence_[
15 print test2, 14 print test2,
16 print test3, 15 print test3,
17 print test4, 16 print test4,
18 test5, 17-- test5,
19 test6, 18-- test6,
20 print test7, 19-- print test7,
21 test8, 20 test8,
22 test0] 21 test0]
23 22
24-- helper functions 23
25vector l = fromList l :: Vector Double
26norm v = pnorm PNorm2 v
27 24
28-- hmatrix vector and matrix 25-- hmatrix vector and matrix
29v = vector [1..10] 26v = vector [1..10]
@@ -34,16 +31,16 @@ m = (5><10) [1..50::Double]
34-- vector creation by in-place updates on a copy of the argument 31-- vector creation by in-place updates on a copy of the argument
35test1 = fun v 32test1 = fun v
36 33
37fun :: Element t => Vector t -> Vector t 34-- fun :: (Num t, Element t, Container) => Vector t -> Vector t
38fun x = runSTVector $ do 35fun x = runSTVector $ do
39 a <- thawVector x 36 a <- thawVector x
40 mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1] 37 mapM_ (flip (modifyVector a) (+57)) [0 .. size x `div` 2 - 1]
41 return a 38 return a
42 39
43-- another example: creation of an antidiagonal matrix from a list 40-- another example: creation of an antidiagonal matrix from a list
44test2 = antiDiag 5 8 [1..] :: Matrix Double 41test2 = antiDiag 5 8 [1..] :: Matrix Double
45 42
46antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b 43-- antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
47antiDiag r c l = runSTMatrix $ do 44antiDiag r c l = runSTMatrix $ do
48 m <- newMatrix 0 r c 45 m <- newMatrix 0 r c
49 let d = min r c - 1 46 let d = min r c - 1
@@ -55,21 +52,23 @@ test3 = g1 v
55 52
56g1 x = runST $ do 53g1 x = runST $ do
57 a <- thawVector x 54 a <- thawVector x
58 writeVector a (dim x -1) 0 55 writeVector a (size x -1) 0
59 b <- freezeVector a 56 b <- freezeVector a
60 return (norm b) 57 return (norm_2 b)
61 58
62-- another possibility: 59-- another possibility:
63test4 = g2 v 60test4 = g2 v
64 61
65g2 x = runST $ do 62g2 x = runST $ do
66 a <- thawVector x 63 a <- thawVector x
67 writeVector a (dim x -1) 0 64 writeVector a (size x -1) 0
68 t <- liftSTVector norm a 65 t <- liftSTVector norm_2 a
69 return t 66 return t
70 67
71-------------------------------------------------------------- 68--------------------------------------------------------------
72 69
70{-
71
73-- haskell arrays 72-- haskell arrays
74hv = listArray (0,9) [1..10::Double] 73hv = listArray (0,9) [1..10::Double]
75hm = listArray ((0,0),(4,9)) [1..50::Double] 74hm = listArray ((0,0),(4,9)) [1..50::Double]
@@ -78,8 +77,8 @@ hm = listArray ((0,0),(4,9)) [1..50::Double]
78 77
79-- conversion from standard Haskell arrays 78-- conversion from standard Haskell arrays
80test5 = do 79test5 = do
81 print $ norm (vectorFromArray hv) 80 print $ norm_2 (vectorFromArray hv)
82 print $ norm v 81 print $ norm_2 v
83 print $ rcond (matrixFromArray hm) 82 print $ rcond (matrixFromArray hm)
84 print $ rcond m 83 print $ rcond m
85 84
@@ -101,10 +100,11 @@ test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)
101 100
102unitary v = runSTUArray $ do 101unitary v = runSTUArray $ do
103 a <- thaw v 102 a <- thaw v
104 n <- norm `fmap` vectorFromMArray a 103 n <- norm_2 `fmap` vectorFromMArray a
105 b <- mapArray (/n) a 104 b <- mapArray (/n) a
106 return b 105 return b
107 106
107-}
108------------------------------------------------- 108-------------------------------------------------
109 109
110-- (just to check that they are not affected) 110-- (just to check that they are not affected)
@@ -150,3 +150,4 @@ test8 = do
150 print $ histoCheck2 ds 150 print $ histoCheck2 ds
151 print $ histoCheck2 ds 151 print $ histoCheck2 ds
152 putStrLn "----------------------" 152 putStrLn "----------------------"
153
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
35flag openblas 35flag 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
39library 40library
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/IO.hs b/packages/base/src/Data/Packed/IO.hs
index 85f1b37..b0999a8 100644
--- a/packages/base/src/Data/Packed/IO.hs
+++ b/packages/base/src/Data/Packed/IO.hs
@@ -22,7 +22,6 @@ import Text.Printf(printf)
22import Data.List(intersperse) 22import Data.List(intersperse)
23import Data.Complex 23import Data.Complex
24import Numeric.Vectorized(vectorScan,saveMatrix) 24import Numeric.Vectorized(vectorScan,saveMatrix)
25import Control.Applicative((<$>))
26import Data.Packed.Internal 25import Data.Packed.Internal
27 26
28{- | Creates a string from a matrix given a separator and a function to show each entry. Using 27{- | Creates a string from a matrix given a separator and a function to show each entry. Using
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
243instance (Container Vector a) => Container Matrix a 243instance (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
17module Numeric.Chain ( 19module 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)
71import qualified Numeric.LinearAlgebra.HMatrix as LA 69import qualified Numeric.LinearAlgebra.HMatrix as LA
72import Data.Proxy(Proxy) 70import Data.Proxy(Proxy)
73import Numeric.LinearAlgebra.Static.Internal 71import Numeric.LinearAlgebra.Static.Internal
@@ -183,6 +181,7 @@ a ¦ b = tr (tr a —— tr b)
183type Sq n = L n n 181type Sq n = L n n
184--type CSq n = CL n n 182--type CSq n = CL n n
185 183
184
186type GL = forall n m . (KnownNat n, KnownNat m) => L m n 185type GL = forall n m . (KnownNat n, KnownNat m) => L m n
187type GSq = forall n . KnownNat n => Sq n 186type 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
307chol :: KnownNat n => Sym n -> Sq n
308chol (extract . unSym -> m) = mkL $ LA.cholSH m
309
308-------------------------------------------------------------------------------- 310--------------------------------------------------------------------------------
309 311
310withNullspace 312withNullspace
@@ -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{- |
620Module : Numeric.LinearAlgebra.Static
621Copyright : (c) Alberto Ruiz 2014
622License : BSD3
623Stability : experimental
624
625Experimental interface with statically checked dimensions.
626
627This module requires GHC >= 7.8
628
629-}
630
631module Numeric.LinearAlgebra.Static
632{-# WARNING "This module requires GHC >= 7.8" #-}
633where
634
635#endif
636
diff --git a/packages/base/src/Numeric/Vectorized.hs b/packages/base/src/Numeric/Vectorized.hs
index 6f0d240..405ae01 100644
--- a/packages/base/src/Numeric/Vectorized.hs
+++ b/packages/base/src/Numeric/Vectorized.hs
@@ -37,7 +37,6 @@ import Foreign.C.String
37import System.IO.Unsafe(unsafePerformIO) 37import System.IO.Unsafe(unsafePerformIO)
38 38
39import Control.Monad(when) 39import Control.Monad(when)
40import Control.Applicative((<$>))
41 40
42 41
43 42
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-11-10 21:48:20 +0000