summaryrefslogtreecommitdiff
path: root/packages/glpk/lib/Numeric/LinearProgramming.hs
diff options
context:
space:
mode:
Diffstat (limited to 'packages/glpk/lib/Numeric/LinearProgramming.hs')
-rw-r--r--packages/glpk/lib/Numeric/LinearProgramming.hs264
1 files changed, 0 insertions, 264 deletions
diff --git a/packages/glpk/lib/Numeric/LinearProgramming.hs b/packages/glpk/lib/Numeric/LinearProgramming.hs
deleted file mode 100644
index 25cdab2..0000000
--- a/packages/glpk/lib/Numeric/LinearProgramming.hs
+++ /dev/null
@@ -1,264 +0,0 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2
3{- |
4Module : Numeric.LinearProgramming
5Copyright : (c) Alberto Ruiz 2010
6License : GPL
7
8Maintainer : Alberto Ruiz
9Stability : provisional
10
11This module provides an interface to the standard simplex algorithm.
12
13For example, the following LP problem
14
15@maximize 4 x_1 - 3 x_2 + 2 x_3
16subject to
17
182 x_1 + x_2 <= 10
19 x_3 + 5 x_4 <= 20
20
21and
22
23x_i >= 0@
24
25can be solved as follows:
26
27@import Numeric.LinearProgramming
28
29prob = Maximize [4, -3, 2]
30
31constr1 = Sparse [ [2\#1, 1\#2] :<=: 10
32 , [1\#2, 5\#3] :<=: 20
33 ]
34
35\> simplex prob constr1 []
36Optimal (28.0,[5.0,0.0,4.0])@
37
38The coefficients of the constraint matrix can also be given in dense format:
39
40@constr2 = Dense [ [2,1,0] :<=: 10
41 , [0,1,5] :<=: 20
42 ]@
43
44By default all variables are bounded as @x_i >= 0@, but this can be
45changed:
46
47@\> simplex prob constr2 [ 2 :=>: 1, 3 :&: (2,7)]
48Optimal (22.6,[4.5,1.0,3.8])
49
50\> simplex prob constr2 [Free 2]
51Unbounded@
52
53The given bound for a variable completely replaces the default,
54so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@.
55Multiple bounds for a variable are not allowed, instead of
56@[i :=>: a, i:<=: b]@ use @i :&: (a,b)@.
57
58-}
59
60module Numeric.LinearProgramming(
61 simplex,
62 Optimization(..),
63 Constraints(..),
64 Bounds,
65 Bound(..),
66 (#),
67 Solution(..)
68) where
69
70import Data.Packed
71import Data.Packed.Development
72import Foreign(Ptr)
73import System.IO.Unsafe(unsafePerformIO)
74import Foreign.C.Types
75import Data.List((\\),sortBy,nub)
76import Data.Function(on)
77
78--import Debug.Trace
79--debug x = trace (show x) x
80
81-----------------------------------------------------
82
83-- | Coefficient of a variable for a sparse representation of constraints.
84(#) :: Double -> Int -> (Double,Int)
85infixl 5 #
86(#) = (,)
87
88data Bound x = x :<=: Double
89 | x :=>: Double
90 | x :&: (Double,Double)
91 | x :==: Double
92 | Free x
93 deriving Show
94
95data Solution = Undefined
96 | Feasible (Double, [Double])
97 | Infeasible (Double, [Double])
98 | NoFeasible
99 | Optimal (Double, [Double])
100 | Unbounded
101 deriving Show
102
103data Constraints = Dense [ Bound [Double] ]
104 | Sparse [ Bound [(Double,Int)] ]
105
106data Optimization = Maximize [Double]
107 | Minimize [Double]
108
109type Bounds = [Bound Int]
110
111simplex :: Optimization -> Constraints -> Bounds -> Solution
112
113simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds
114simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds
115
116simplex opt (Dense constr) bnds = extract sg sol where
117 sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)
118 n = length objfun
119 m = length constr
120 (sz, sg, objfun) = adapt opt
121
122simplex opt (Sparse constr) bnds = extract sg sol where
123 sol = simplexSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)
124 n = length objfun
125 m = length constr
126 (sz, sg, objfun) = adapt opt
127
128adapt :: Optimization -> (Int, Double, [Double])
129adapt opt = case opt of
130 Maximize x -> (size x, 1 ,x)
131 Minimize x -> (size x, -1, (map negate x))
132 where size x | null x = error "simplex: objective function with zero variables"
133 | otherwise = length x
134
135extract :: Double -> Vector Double -> Solution
136extract sg sol = r where
137 z = sg * (sol@>1)
138 v = toList $ subVector 2 (dim sol -2) sol
139 r = case round(sol@>0)::Int of
140 1 -> Undefined
141 2 -> Feasible (z,v)
142 3 -> Infeasible (z,v)
143 4 -> NoFeasible
144 5 -> Optimal (z,v)
145 6 -> Unbounded
146 _ -> error "simplex: solution type unknown"
147
148-----------------------------------------------------
149
150obj :: Bound t -> t
151obj (x :<=: _) = x
152obj (x :=>: _) = x
153obj (x :&: _) = x
154obj (x :==: _) = x
155obj (Free x) = x
156
157tb :: Bound t -> Double
158tb (_ :<=: _) = glpUP
159tb (_ :=>: _) = glpLO
160tb (_ :&: _) = glpDB
161tb (_ :==: _) = glpFX
162tb (Free _) = glpFR
163
164lb :: Bound t -> Double
165lb (_ :<=: _) = 0
166lb (_ :=>: a) = a
167lb (_ :&: (a,_)) = a
168lb (_ :==: a) = a
169lb (Free _) = 0
170
171ub :: Bound t -> Double
172ub (_ :<=: a) = a
173ub (_ :=>: _) = 0
174ub (_ :&: (_,a)) = a
175ub (_ :==: a) = a
176ub (Free _) = 0
177
178mkBound1 :: Bound t -> [Double]
179mkBound1 b = [tb b, lb b, ub b]
180
181mkBound2 :: Bound t -> (t, [Double])
182mkBound2 b = (obj b, mkBound1 b)
183
184mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double
185mkBounds n b1 b2 = fromLists (cb++vb) where
186 gv' = map obj b2
187 gv | nub gv' == gv' = gv'
188 | otherwise = error $ "simplex: duplicate bounds for vars " ++ show (gv'\\nub gv')
189 rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv
190 | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n
191 vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:=>: 0)) rv ++ map mkBound2 b2
192 cb = map mkBound1 b1
193
194mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double
195mkConstrD n f b1 | ok = fromLists (ob ++ co)
196 | otherwise = error $ "simplex: dense constraints require "++show n
197 ++" variables, given " ++ show ls
198 where
199 cs = map obj b1
200 ls = map length cs
201 ok = all (==n) ls
202 den = fromLists cs
203 ob = map (([0,0]++).return) f
204 co = [[fromIntegral i, fromIntegral j,den@@>(i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]
205
206mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double
207mkConstrS n objfun b1 = fromLists (ob ++ co) where
208 ob = map (([0,0]++).return) objfun
209 co = concat $ zipWith f [1::Int ..] cs
210 cs = map obj b1
211 f k = map (g k)
212 g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c]
213 | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n
214
215-----------------------------------------------------
216
217foreign import ccall unsafe "c_simplex_sparse" c_simplex_sparse
218 :: CInt -> CInt -- rows and cols
219 -> CInt -> CInt -> Ptr Double -- coeffs
220 -> CInt -> CInt -> Ptr Double -- bounds
221 -> CInt -> Ptr Double -- result
222 -> IO CInt -- exit code
223
224simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double
225simplexSparse m n c b = unsafePerformIO $ do
226 s <- createVector (2+n)
227 app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"
228 return s
229
230glpFR, glpLO, glpUP, glpDB, glpFX :: Double
231glpFR = 0
232glpLO = 1
233glpUP = 2
234glpDB = 3
235glpFX = 4
236
237{- Raw format of coeffs
238
239simplexSparse
240
241(12><3)
242 [ 0.0, 0.0, 10.0
243 , 0.0, 0.0, 6.0
244 , 0.0, 0.0, 4.0
245 , 1.0, 1.0, 1.0
246 , 1.0, 2.0, 1.0
247 , 1.0, 3.0, 1.0
248 , 2.0, 1.0, 10.0
249 , 2.0, 2.0, 4.0
250 , 2.0, 3.0, 5.0
251 , 3.0, 1.0, 2.0
252 , 3.0, 2.0, 2.0
253 , 3.0, 3.0, 6.0 ]
254
255bounds = (6><3)
256 [ glpUP,0,100
257 , glpUP,0,600
258 , glpUP,0,300
259 , glpLO,0,0
260 , glpLO,0,0
261 , glpLO,0,0 ]
262
263-}
264