diff options
Diffstat (limited to 'packages/glpk/src/Numeric')
-rw-r--r-- | packages/glpk/src/Numeric/LinearProgramming.hs | 264 | ||||
-rw-r--r-- | packages/glpk/src/Numeric/LinearProgramming/L1.hs | 67 |
2 files changed, 331 insertions, 0 deletions
diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs new file mode 100644 index 0000000..25cdab2 --- /dev/null +++ b/packages/glpk/src/Numeric/LinearProgramming.hs | |||
@@ -0,0 +1,264 @@ | |||
1 | {-# LANGUAGE ForeignFunctionInterface #-} | ||
2 | |||
3 | {- | | ||
4 | Module : Numeric.LinearProgramming | ||
5 | Copyright : (c) Alberto Ruiz 2010 | ||
6 | License : GPL | ||
7 | |||
8 | Maintainer : Alberto Ruiz | ||
9 | Stability : provisional | ||
10 | |||
11 | This module provides an interface to the standard simplex algorithm. | ||
12 | |||
13 | For example, the following LP problem | ||
14 | |||
15 | @maximize 4 x_1 - 3 x_2 + 2 x_3 | ||
16 | subject to | ||
17 | |||
18 | 2 x_1 + x_2 <= 10 | ||
19 | x_3 + 5 x_4 <= 20 | ||
20 | |||
21 | and | ||
22 | |||
23 | x_i >= 0@ | ||
24 | |||
25 | can be solved as follows: | ||
26 | |||
27 | @import Numeric.LinearProgramming | ||
28 | |||
29 | prob = Maximize [4, -3, 2] | ||
30 | |||
31 | constr1 = Sparse [ [2\#1, 1\#2] :<=: 10 | ||
32 | , [1\#2, 5\#3] :<=: 20 | ||
33 | ] | ||
34 | |||
35 | \> simplex prob constr1 [] | ||
36 | Optimal (28.0,[5.0,0.0,4.0])@ | ||
37 | |||
38 | The 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 | |||
44 | By default all variables are bounded as @x_i >= 0@, but this can be | ||
45 | changed: | ||
46 | |||
47 | @\> simplex prob constr2 [ 2 :=>: 1, 3 :&: (2,7)] | ||
48 | Optimal (22.6,[4.5,1.0,3.8]) | ||
49 | |||
50 | \> simplex prob constr2 [Free 2] | ||
51 | Unbounded@ | ||
52 | |||
53 | The given bound for a variable completely replaces the default, | ||
54 | so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@. | ||
55 | Multiple bounds for a variable are not allowed, instead of | ||
56 | @[i :=>: a, i:<=: b]@ use @i :&: (a,b)@. | ||
57 | |||
58 | -} | ||
59 | |||
60 | module Numeric.LinearProgramming( | ||
61 | simplex, | ||
62 | Optimization(..), | ||
63 | Constraints(..), | ||
64 | Bounds, | ||
65 | Bound(..), | ||
66 | (#), | ||
67 | Solution(..) | ||
68 | ) where | ||
69 | |||
70 | import Data.Packed | ||
71 | import Data.Packed.Development | ||
72 | import Foreign(Ptr) | ||
73 | import System.IO.Unsafe(unsafePerformIO) | ||
74 | import Foreign.C.Types | ||
75 | import Data.List((\\),sortBy,nub) | ||
76 | import 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) | ||
85 | infixl 5 # | ||
86 | (#) = (,) | ||
87 | |||
88 | data Bound x = x :<=: Double | ||
89 | | x :=>: Double | ||
90 | | x :&: (Double,Double) | ||
91 | | x :==: Double | ||
92 | | Free x | ||
93 | deriving Show | ||
94 | |||
95 | data Solution = Undefined | ||
96 | | Feasible (Double, [Double]) | ||
97 | | Infeasible (Double, [Double]) | ||
98 | | NoFeasible | ||
99 | | Optimal (Double, [Double]) | ||
100 | | Unbounded | ||
101 | deriving Show | ||
102 | |||
103 | data Constraints = Dense [ Bound [Double] ] | ||
104 | | Sparse [ Bound [(Double,Int)] ] | ||
105 | |||
106 | data Optimization = Maximize [Double] | ||
107 | | Minimize [Double] | ||
108 | |||
109 | type Bounds = [Bound Int] | ||
110 | |||
111 | simplex :: Optimization -> Constraints -> Bounds -> Solution | ||
112 | |||
113 | simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds | ||
114 | simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds | ||
115 | |||
116 | simplex 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 | |||
122 | simplex 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 | |||
128 | adapt :: Optimization -> (Int, Double, [Double]) | ||
129 | adapt 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 | |||
135 | extract :: Double -> Vector Double -> Solution | ||
136 | extract 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 | |||
150 | obj :: Bound t -> t | ||
151 | obj (x :<=: _) = x | ||
152 | obj (x :=>: _) = x | ||
153 | obj (x :&: _) = x | ||
154 | obj (x :==: _) = x | ||
155 | obj (Free x) = x | ||
156 | |||
157 | tb :: Bound t -> Double | ||
158 | tb (_ :<=: _) = glpUP | ||
159 | tb (_ :=>: _) = glpLO | ||
160 | tb (_ :&: _) = glpDB | ||
161 | tb (_ :==: _) = glpFX | ||
162 | tb (Free _) = glpFR | ||
163 | |||
164 | lb :: Bound t -> Double | ||
165 | lb (_ :<=: _) = 0 | ||
166 | lb (_ :=>: a) = a | ||
167 | lb (_ :&: (a,_)) = a | ||
168 | lb (_ :==: a) = a | ||
169 | lb (Free _) = 0 | ||
170 | |||
171 | ub :: Bound t -> Double | ||
172 | ub (_ :<=: a) = a | ||
173 | ub (_ :=>: _) = 0 | ||
174 | ub (_ :&: (_,a)) = a | ||
175 | ub (_ :==: a) = a | ||
176 | ub (Free _) = 0 | ||
177 | |||
178 | mkBound1 :: Bound t -> [Double] | ||
179 | mkBound1 b = [tb b, lb b, ub b] | ||
180 | |||
181 | mkBound2 :: Bound t -> (t, [Double]) | ||
182 | mkBound2 b = (obj b, mkBound1 b) | ||
183 | |||
184 | mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double | ||
185 | mkBounds 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 | |||
194 | mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double | ||
195 | mkConstrD 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 | |||
206 | mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double | ||
207 | mkConstrS 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 | |||
217 | foreign 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 | |||
224 | simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double | ||
225 | simplexSparse 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 | |||
230 | glpFR, glpLO, glpUP, glpDB, glpFX :: Double | ||
231 | glpFR = 0 | ||
232 | glpLO = 1 | ||
233 | glpUP = 2 | ||
234 | glpDB = 3 | ||
235 | glpFX = 4 | ||
236 | |||
237 | {- Raw format of coeffs | ||
238 | |||
239 | simplexSparse | ||
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 | |||
255 | bounds = (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 | |||
diff --git a/packages/glpk/src/Numeric/LinearProgramming/L1.hs b/packages/glpk/src/Numeric/LinearProgramming/L1.hs new file mode 100644 index 0000000..2e8f94a --- /dev/null +++ b/packages/glpk/src/Numeric/LinearProgramming/L1.hs | |||
@@ -0,0 +1,67 @@ | |||
1 | {- | | ||
2 | Module : Numeric.LinearProgramming.L1 | ||
3 | Copyright : (c) Alberto Ruiz 2011-14 | ||
4 | Stability : provisional | ||
5 | |||
6 | Linear system solvers in the L_1 norm using linear programming. | ||
7 | |||
8 | -} | ||
9 | ----------------------------------------------------------------------------- | ||
10 | |||
11 | module Numeric.LinearProgramming.L1 ( | ||
12 | l1SolveO, lInfSolveO, | ||
13 | l1SolveU, | ||
14 | ) where | ||
15 | |||
16 | import Numeric.LinearAlgebra | ||
17 | import Numeric.LinearProgramming | ||
18 | |||
19 | -- | L_Inf solution of overconstrained system Ax=b. | ||
20 | -- | ||
21 | -- Find argmin x ||Ax-b||_inf | ||
22 | lInfSolveO :: Matrix Double -> Vector Double -> Vector Double | ||
23 | lInfSolveO a b = fromList (take n x) | ||
24 | where | ||
25 | n = cols a | ||
26 | as = toRows a | ||
27 | bs = toList b | ||
28 | c1 = zipWith (mk (1)) as bs | ||
29 | c2 = zipWith (mk (-1)) as bs | ||
30 | mk sign a_i b_i = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(n+1)]) :<=: (sign * b_i) | ||
31 | p = Sparse (c1++c2) | ||
32 | Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ [1])) p (map Free [1..(n+1)]) | ||
33 | |||
34 | |||
35 | -- | L_1 solution of overconstrained system Ax=b. | ||
36 | -- | ||
37 | -- Find argmin x ||Ax-b||_1. | ||
38 | l1SolveO :: Matrix Double -> Vector Double -> Vector Double | ||
39 | l1SolveO a b = fromList (take n x) | ||
40 | where | ||
41 | n = cols a | ||
42 | m = rows a | ||
43 | as = toRows a | ||
44 | bs = toList b | ||
45 | ks = [1..] | ||
46 | c1 = zipWith3 (mk (1)) as bs ks | ||
47 | c2 = zipWith3 (mk (-1)) as bs ks | ||
48 | mk sign a_i b_i k = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(k+n)]) :<=: (sign * b_i) | ||
49 | p = Sparse (c1++c2) | ||
50 | Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ replicate m 1)) p (map Free [1..(n+m)]) | ||
51 | |||
52 | |||
53 | |||
54 | -- | L1 solution of underconstrained linear system Ax=b. | ||
55 | -- | ||
56 | -- Find argmin x ||x||_1 such that Ax=b. | ||
57 | l1SolveU :: Matrix Double -> Vector Double -> Vector Double | ||
58 | l1SolveU a y = fromList (take n x) | ||
59 | where | ||
60 | n = cols a | ||
61 | c1 = map (\k -> [ 1#k, -1#k+n] :<=: 0) [1..n] | ||
62 | c2 = map (\k -> [-1#k, -1#k+n] :<=: 0) [1..n] | ||
63 | c3 = zipWith (:==:) (map sp $ toRows a) (toList y) | ||
64 | sp v = zipWith (#) (toList v) [1..] | ||
65 | p = Sparse (c1 ++ c2 ++ c3) | ||
66 | Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ replicate n 1)) p (map Free [1..(2*n)]) | ||
67 | |||