summaryrefslogtreecommitdiff
path: root/packages/glpk/lib/Numeric/LinearProgramming.hs
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2010-02-21 18:26:23 +0000
committerAlberto Ruiz <aruiz@um.es>2010-02-21 18:26:23 +0000
commitf38b4a3076cfae023559ce61cb2a443c809b7a6f (patch)
tree022c127181fb65c34705cdcf44221b4ac89ba50b /packages/glpk/lib/Numeric/LinearProgramming.hs
parenta3d1bb34ae7b1f97b7e9900fc38f145094fe4777 (diff)
simple glpk interface
Diffstat (limited to 'packages/glpk/lib/Numeric/LinearProgramming.hs')
-rw-r--r--packages/glpk/lib/Numeric/LinearProgramming.hs211
1 files changed, 211 insertions, 0 deletions
diff --git a/packages/glpk/lib/Numeric/LinearProgramming.hs b/packages/glpk/lib/Numeric/LinearProgramming.hs
new file mode 100644
index 0000000..1b14080
--- /dev/null
+++ b/packages/glpk/lib/Numeric/LinearProgramming.hs
@@ -0,0 +1,211 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2
3module Numeric.LinearProgramming(
4 Optimization(..),
5 Bound(..),
6 (#),
7 Coeffs(..),
8 Solution(..),
9 simplex,
10) where
11
12import Numeric.LinearAlgebra
13import Data.Packed.Development
14import Foreign(Ptr,unsafePerformIO)
15import Foreign.C.Types(CInt)
16import Data.List((\\),sortBy)
17import Data.Function(on)
18
19--import Debug.Trace
20--debug x = trace (show x) x
21
22-----------------------------------------------------
23
24-- | Coefficient of a variable
25(#) :: Double -> Int -> (Double,Int)
26infixl 5 #
27(#) = (,)
28
29data Bound x = x :<: Double
30 | x :>: Double
31 | x :&: (Double,Double)
32 | x :==: Double
33 | Free x
34 deriving Show
35
36data Solution = Undefined
37 | Feasible (Double, [Double])
38 | Infeasible (Double, [Double])
39 | NoFeasible
40 | Optimal (Double, [Double])
41 | Unbounded
42 deriving Show
43
44data Coeffs = Dense [ Bound [Double] ]
45 | Sparse [ Bound [(Double,Int)] ]
46
47data Optimization = Maximize [Double]
48 | Minimize [Double]
49
50simplex :: Optimization -> Coeffs -> [Bound Int] -> Solution
51
52simplex opt (Dense constr) bnds = extract sg sol where
53 sol = simplexDense (mkConstrD objfun constr) (mkBoundsD constr bnds)
54 (sg, objfun) = case opt of
55 Maximize x -> (1 ,x)
56 Minimize x -> (-1, (map negate x))
57
58simplex opt (Sparse constr) bnds = extract sg sol where
59 sol = simplexSparse m n (mkConstrS objfun constr) (mkBoundsS constr bnds)
60 n = length objfun
61 m = length constr
62 (sg, objfun) = case opt of
63 Maximize x -> (1 ,x)
64 Minimize x -> (-1, (map negate x))
65
66extract :: Double -> Vector Double -> Solution
67extract sg sol = r where
68 z = sg * (sol@>1)
69 v = toList $ subVector 2 (dim sol -2) sol
70 r = case round(sol@>0)::Int of
71 1 -> Undefined
72 2 -> Feasible (z,v)
73 3 -> Infeasible (z,v)
74 4 -> NoFeasible
75 5 -> Optimal (z,v)
76 6 -> Unbounded
77 _ -> error "simplex: solution type unknown"
78
79-----------------------------------------------------
80
81obj :: Bound t -> t
82obj (x :<: _) = x
83obj (x :>: _) = x
84obj (x :&: _) = x
85obj (x :==: _) = x
86obj (Free x) = x
87
88tb :: Bound t -> Double
89tb (_ :<: _) = glpUP
90tb (_ :>: _) = glpLO
91tb (_ :&: _) = glpDB
92tb (_ :==: _) = glpFX
93tb (Free _) = glpFR
94
95lb :: Bound t -> Double
96lb (_ :<: _) = 0
97lb (_ :>: a) = a
98lb (_ :&: (a,_)) = a
99lb (_ :==: a) = a
100lb (Free _) = 0
101
102ub :: Bound t -> Double
103ub (_ :<: a) = a
104ub (_ :>: _) = 0
105ub (_ :&: (_,a)) = a
106ub (_ :==: a) = a
107ub (Free _) = 0
108
109mkBound1 :: Bound t -> [Double]
110mkBound1 b = [tb b, lb b, ub b]
111
112mkBound2 :: Bound t -> (t, [Double])
113mkBound2 b = (obj b, mkBound1 b)
114
115mkBoundsD :: [Bound [a]] -> [Bound Int] -> Matrix Double
116mkBoundsD b1 b2 = fromLists (cb++vb) where
117 c = length (obj (head b1))
118 gv = map obj b2
119 rv = [1..c] \\ gv
120 vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2
121 cb = map mkBound1 b1
122
123mkConstrD :: [Double] -> [Bound [Double]] -> Matrix Double
124mkConstrD f b1 = fromLists (f : map obj b1)
125
126mkBoundsS :: [Bound [(Double, Int)]] -> [Bound Int] -> Matrix Double
127mkBoundsS b1 b2 = fromLists (cb++vb) where
128 c = maximum $ map snd $ concatMap obj b1
129 gv = map obj b2
130 rv = [1..c] \\ gv
131 vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2
132 cb = map mkBound1 b1
133
134mkConstrS :: [Double] -> [Bound [(Double, Int)]] -> Matrix Double
135mkConstrS objfun constr = fromLists (ob ++ co) where
136 ob = map (([0,0]++).return) objfun
137 co = concat $ zipWith f [1::Int ..] (map obj constr)
138 f k = map (g k)
139 g k (c,v) = [fromIntegral k, fromIntegral v,c]
140
141-----------------------------------------------------
142
143foreign import ccall "c_simplex_dense" c_simplex_dense
144 :: CInt -> CInt -> Ptr Double -- coeffs
145 -> CInt -> CInt -> Ptr Double -- bounds
146 -> CInt -> Ptr Double -- result
147 -> IO CInt -- exit code
148
149simplexDense :: Matrix Double -> Matrix Double -> Vector Double
150simplexDense c b = unsafePerformIO $ do
151 s <- createVector (2+cols c)
152 app3 c_simplex_dense mat (cmat c) mat (cmat b) vec s "c_simplex_dense"
153 return s
154
155foreign import ccall "c_simplex_sparse" c_simplex_sparse
156 :: CInt -> CInt -- rows and cols
157 -> CInt -> CInt -> Ptr Double -- coeffs
158 -> CInt -> CInt -> Ptr Double -- bounds
159 -> CInt -> Ptr Double -- result
160 -> IO CInt -- exit code
161
162simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
163simplexSparse m n c b = unsafePerformIO $ do
164 s <- createVector (2+n)
165 let fi = fromIntegral
166 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
167 return s
168
169glpFR, glpLO, glpUP, glpDB, glpFX :: Double
170glpFR = 0
171glpLO = 1
172glpUP = 2
173glpDB = 3
174glpFX = 4
175
176{- Raw format of coeffs
177
178simplexDense
179
180((3+1) >< 3)
181 [ 10, 6, 4
182 , 1, 1, 1
183 , 10, 4, 5
184 , 2, 2, 6 :: Double]
185
186simplexSparse
187
188(12><3)
189 [ 0.0, 0.0, 10.0
190 , 0.0, 0.0, 6.0
191 , 0.0, 0.0, 4.0
192 , 1.0, 1.0, 1.0
193 , 1.0, 2.0, 1.0
194 , 1.0, 3.0, 1.0
195 , 2.0, 1.0, 10.0
196 , 2.0, 2.0, 4.0
197 , 2.0, 3.0, 5.0
198 , 3.0, 1.0, 2.0
199 , 3.0, 2.0, 2.0
200 , 3.0, 3.0, 6.0 ]
201
202bounds = (6><3)
203 [ glpUP,0,100
204 , glpUP,0,600
205 , glpUP,0,300
206 , glpLO,0,0
207 , glpLO,0,0
208 , glpLO,0,0 ]
209
210-}
211