hmatrix-glpk (empty) → 0.1.0
raw patch · 9 files changed
+469/−0 lines, 9 filesdep +basedep +hmatrixsetup-changed
Dependencies added: base, hmatrix
Files
- LICENSE +2/−0
- Setup.lhs +4/−0
- examples/simplex1.hs +20/−0
- examples/simplex2.hs +18/−0
- examples/simplex3.hs +25/−0
- examples/simplex4.hs +25/−0
- hmatrix-glpk.cabal +35/−0
- lib/Numeric/LinearProgramming.hs +264/−0
- lib/Numeric/LinearProgramming/glpk.c +76/−0
+ LICENSE view
@@ -0,0 +1,2 @@+Copyright Alberto Ruiz 2010+GPL license
+ Setup.lhs view
@@ -0,0 +1,4 @@+#! /usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain
+ examples/simplex1.hs view
@@ -0,0 +1,20 @@+-- first example in glpk manual++import Numeric.LinearProgramming++objFun = Maximize [10, 6, 4]++constr = Dense [ [1,1,1] :<: 100+ , [10,4,5] :<: 600+ , [2,2,6] :<: 300 ]++-- default bounds+bnds = [ 1 :>: 0+ , 2 :>: 0+ , 3 :>: 0 ]++main = do+ print $ simplex objFun constr []+ print $ simplex objFun constr bnds+ print $ simplex objFun constr [Free 3]+ print $ simplex objFun constr [ 2 :<: 50 ]
+ examples/simplex2.hs view
@@ -0,0 +1,18 @@+import Numeric.LinearProgramming++prob = Maximize [4, -3, 2]++constr1 = Sparse [ [2#1, 1#2] :<: 10+ , [1#2, 5#3] :<: 20+ ]++constr2 = Dense [ [2,1,0] :<: 10+ , [0,1,5] :<: 20+ ]++main = do+ print $ simplex prob constr1 []+ print $ simplex prob constr2 []+ print $ simplex prob constr2 [ 2 :>: 1, 3 :&: (2,7)]+ print $ simplex prob constr2 [ Free 2 ]+
+ examples/simplex3.hs view
@@ -0,0 +1,25 @@+-- compare with+-- $ glpsol --cpxlp /usr/share/doc/glpk-utils/examples/plan.lp -o result.txt++import Numeric.LinearProgramming++prob = Minimize [0.03, 0.08, 0.17, 0.12, 0.15, 0.21, 0.38]++constr = Dense+ [ [1,1,1,1,1,1,1] :==: 2000+ , [0.15, 0.04, 0.02, 0.04, 0.2,0.01, 0.03] :<: 60+ , [0.03, 0.05, 0.08, 0.02, 0.06, 0.01, 0] :<: 100+ , [0.02, 0.04, 0.01, 0.02, 0.02, 0, 0] :<: 40+ , [0.02, 0.03, 0, 0, 0.01, 0, 0] :<: 30+ , [0.7, 0.75, 0.8, 0.75, 0.8, 0.97, 0] :>: 1500+ , [0.02, 0.06, 0.08, 0.12, 0.02, 0.01, 0.97] :&: (250,300)+ ]++bounds = [ 1 :&: (0,200)+ , 2 :&: (0,2500)+ , 3 :&: (400,800)+ , 4 :&: (100,700)+ , 5 :&: (0,1500) ]++main = print $ simplex prob constr bounds+
+ examples/simplex4.hs view
@@ -0,0 +1,25 @@+-- compare with+-- $ glpsol --cpxlp /usr/share/doc/glpk-utils/examples/plan.lp -o result.txt++import Numeric.LinearProgramming++prob = Minimize [0.03, 0.08, 0.17, 0.12, 0.15, 0.21, 0.38]++constr = Sparse+ [ [1#1,1#2,1#3,1#4,1#5,1#6,1#7] :==: 2000+ , [0.15#1, 0.04#2, 0.02#3, 0.04#4, 0.2#5,0.01#6, 0.03#7] :<: 60+ , [0.03#1, 0.05#2, 0.08#3, 0.02#4, 0.06#5, 0.01#6] :<: 100+ , [0.02#1, 0.04#2, 0.01#3, 0.02#4, 0.02#5] :<: 40+ , [0.02#1, 0.03#2, 0.01#5] :<: 30+ , [0.7#1, 0.75#2, 0.8#3, 0.75#4, 0.8#5, 0.97#6] :>: 1500+ , [0.02#1, 0.06#2, 0.08#3, 0.12#4, 0.02#5, 0.01#6, 0.97#7] :&: (250,300)+ ]++bounds = [ 1 :&: (0,200)+ , 2 :&: (0,2500)+ , 3 :&: (400,800)+ , 4 :&: (100,700)+ , 5 :&: (0,1500) ]++main = print $ simplex prob constr bounds+
+ hmatrix-glpk.cabal view
@@ -0,0 +1,35 @@+Name: hmatrix-glpk+Version: 0.1.0+License: GPL+License-file: LICENSE+Author: Alberto Ruiz+Maintainer: Alberto Ruiz <aruiz@um.es>+Stability: experimental+Homepage: http://code.haskell.org/hmatrix+Synopsis: Linear Programming based on GLPK+Description:+ Simple interface to linear programming functions provided by GLPK.++Category: Math+tested-with: GHC ==6.10.4++cabal-version: >=1.2+build-type: Simple++extra-source-files: examples/simplex1.hs+ examples/simplex2.hs+ examples/simplex3.hs+ examples/simplex4.hs++library+ Build-Depends: base >= 3 && < 5, hmatrix >= 0.8.3 && < 0.9++ hs-source-dirs: lib++ Exposed-modules: Numeric.LinearProgramming++ c-sources: lib/Numeric/LinearProgramming/glpk.c++ ghc-options: -Wall++ extra-libraries: glpk
+ lib/Numeric/LinearProgramming.hs view
@@ -0,0 +1,264 @@+{-# LANGUAGE ForeignFunctionInterface #-}++{- |+Module : Numeric.LinearProgramming+Copyright : (c) Alberto Ruiz 2010+License : GPL++Maintainer : Alberto Ruiz (aruiz at um dot es)+Stability : provisional++This module provides an interface to the standard simplex algorithm.++For example, the following LP problem++@maximize 4 x_1 - 3 x_2 + 2 x_3+subject to++2 x_1 + x_2 <= 10+ x_3 + 5 x_4 <= 20++and++x_i >= 0@++can be solved as follows:++@import Numeric.LinearProgramming++prob = Maximize [4, -3, 2]++constr1 = Sparse [ [2\#1, 1\#2] :<: 10+ , [1\#2, 5\#3] :<: 20+ ]++\> simplex prob constr1 []+Optimal (28.0,[5.0,0.0,4.0])@++The coefficients of the constraint matrix can also be given in dense format:++@constr2 = Dense [ [2,1,0] :<: 10+ , [0,1,5] :<: 20+ ]@++By default all variables are bounded as @x_i >= 0@, but this can be+changed:++@\> simplex prob constr2 [ 2 :>: 1, 3 :&: (2,7)]+Optimal (22.6,[4.5,1.0,3.8])++\> simplex prob constr2 [Free 2]+Unbounded@++The given bound for a variable completely replaces the default,+so @0 <= x_i <= b@ must be explicitly given as @i :&: (0,b)@.+Multiple bounds for a variable are not allowed, instead of+@[i :>: a, i:<: b]@ use @i :&: (a,b)@.++-}++module Numeric.LinearProgramming(+ simplex,+ Optimization(..),+ Constraints(..),+ Bounds,+ Bound(..),+ (#),+ Solution(..)+) where++import Numeric.LinearAlgebra hiding (i)+import Data.Packed.Development+import Foreign(Ptr,unsafePerformIO)+import Foreign.C.Types(CInt)+import Data.List((\\),sortBy,nub)+import Data.Function(on)++--import Debug.Trace+--debug x = trace (show x) x++-----------------------------------------------------++-- | Coefficient of a variable for a sparse representation of constraints.+(#) :: Double -> Int -> (Double,Int)+infixl 5 #+(#) = (,)++data Bound x = x :<: Double+ | x :>: Double+ | x :&: (Double,Double)+ | x :==: Double+ | Free x+ deriving Show++data Solution = Undefined+ | Feasible (Double, [Double])+ | Infeasible (Double, [Double])+ | NoFeasible+ | Optimal (Double, [Double])+ | Unbounded+ deriving Show++data Constraints = Dense [ Bound [Double] ]+ | Sparse [ Bound [(Double,Int)] ]++data Optimization = Maximize [Double]+ | Minimize [Double]++type Bounds = [Bound Int]++simplex :: Optimization -> Constraints -> Bounds -> Solution++simplex opt (Dense []) bnds = simplex opt (Sparse []) bnds+simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds++simplex opt (Dense constr) bnds = extract sg sol where+ sol = simplexSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)+ n = length objfun+ m = length constr+ (sz, sg, objfun) = adapt opt++simplex opt (Sparse constr) bnds = extract sg sol where+ sol = simplexSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)+ n = length objfun+ m = length constr+ (sz, sg, objfun) = adapt opt++adapt :: Optimization -> (Int, Double, [Double])+adapt opt = case opt of+ Maximize x -> (size x, 1 ,x)+ Minimize x -> (size x, -1, (map negate x))+ where size x | null x = error "simplex: objective function with zero variables"+ | otherwise = length x++extract :: Double -> Vector Double -> Solution+extract sg sol = r where+ z = sg * (sol@>1)+ v = toList $ subVector 2 (dim sol -2) sol+ r = case round(sol@>0)::Int of+ 1 -> Undefined+ 2 -> Feasible (z,v)+ 3 -> Infeasible (z,v)+ 4 -> NoFeasible+ 5 -> Optimal (z,v)+ 6 -> Unbounded+ _ -> error "simplex: solution type unknown"++-----------------------------------------------------++obj :: Bound t -> t+obj (x :<: _) = x+obj (x :>: _) = x+obj (x :&: _) = x+obj (x :==: _) = x+obj (Free x) = x++tb :: Bound t -> Double+tb (_ :<: _) = glpUP+tb (_ :>: _) = glpLO+tb (_ :&: _) = glpDB+tb (_ :==: _) = glpFX+tb (Free _) = glpFR++lb :: Bound t -> Double+lb (_ :<: _) = 0+lb (_ :>: a) = a+lb (_ :&: (a,_)) = a+lb (_ :==: a) = a+lb (Free _) = 0++ub :: Bound t -> Double+ub (_ :<: a) = a+ub (_ :>: _) = 0+ub (_ :&: (_,a)) = a+ub (_ :==: a) = a+ub (Free _) = 0++mkBound1 :: Bound t -> [Double]+mkBound1 b = [tb b, lb b, ub b]++mkBound2 :: Bound t -> (t, [Double])+mkBound2 b = (obj b, mkBound1 b)++mkBounds :: Int -> [Bound [a]] -> [Bound Int] -> Matrix Double+mkBounds n b1 b2 = fromLists (cb++vb) where+ gv' = map obj b2+ gv | nub gv' == gv' = gv'+ | otherwise = error $ "simplex: duplicate bounds for vars " ++ show (gv'\\nub gv')+ rv | null gv || minimum gv >= 0 && maximum gv <= n = [1..n] \\ gv+ | otherwise = error $ "simplex: bounds: variables "++show gv++" not in 1.."++show n+ vb = map snd $ sortBy (compare `on` fst) $ map (mkBound2 . (:>: 0)) rv ++ map mkBound2 b2+ cb = map mkBound1 b1++mkConstrD :: Int -> [Double] -> [Bound [Double]] -> Matrix Double+mkConstrD n f b1 | ok = fromLists (ob ++ co)+ | otherwise = error $ "simplex: dense constraints require "++show n+ ++" variables, given " ++ show ls+ where+ cs = map obj b1+ ls = map length cs+ ok = all (==n) ls+ den = fromLists cs+ ob = map (([0,0]++).return) f+ co = [[fromIntegral i, fromIntegral j,den@@>(i-1,j-1)]| i<-[1 ..rows den], j<-[1 .. cols den]]++mkConstrS :: Int -> [Double] -> [Bound [(Double, Int)]] -> Matrix Double+mkConstrS n objfun b1 = fromLists (ob ++ co) where+ ob = map (([0,0]++).return) objfun+ co = concat $ zipWith f [1::Int ..] cs+ cs = map obj b1+ f k = map (g k)+ g k (c,v) | v >=1 && v<= n = [fromIntegral k, fromIntegral v,c]+ | otherwise = error $ "simplex: sparse constraints: variable "++show v++" not in 1.."++show n++-----------------------------------------------------++foreign import ccall "c_simplex_sparse" c_simplex_sparse+ :: CInt -> CInt -- rows and cols+ -> CInt -> CInt -> Ptr Double -- coeffs+ -> CInt -> CInt -> Ptr Double -- bounds+ -> CInt -> Ptr Double -- result+ -> IO CInt -- exit code++simplexSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double+simplexSparse m n c b = unsafePerformIO $ do+ s <- createVector (2+n)+ let fi = fromIntegral+ app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"+ return s++glpFR, glpLO, glpUP, glpDB, glpFX :: Double+glpFR = 0+glpLO = 1+glpUP = 2+glpDB = 3+glpFX = 4++{- Raw format of coeffs++simplexSparse++(12><3)+ [ 0.0, 0.0, 10.0+ , 0.0, 0.0, 6.0+ , 0.0, 0.0, 4.0+ , 1.0, 1.0, 1.0+ , 1.0, 2.0, 1.0+ , 1.0, 3.0, 1.0+ , 2.0, 1.0, 10.0+ , 2.0, 2.0, 4.0+ , 2.0, 3.0, 5.0+ , 3.0, 1.0, 2.0+ , 3.0, 2.0, 2.0+ , 3.0, 3.0, 6.0 ]++bounds = (6><3)+ [ glpUP,0,100+ , glpUP,0,600+ , glpUP,0,300+ , glpLO,0,0+ , glpLO,0,0+ , glpLO,0,0 ]++-}+
+ lib/Numeric/LinearProgramming/glpk.c view
@@ -0,0 +1,76 @@+#define DVEC(A) int A##n, double*A##p+#define DMAT(A) int A##r, int A##c, double*A##p++#define AT(M,r,co) (M##p[(r)*M##c+(co)])++#include <stdlib.h>+#include <stdio.h>+#include <glpk.h>+#include <math.h>++/*-----------------------------------------------------*/++int c_simplex_sparse(int m, int n, DMAT(c), DMAT(b), DVEC(s)) {+ glp_prob *lp;+ lp = glp_create_prob();+ glp_set_obj_dir(lp, GLP_MAX);+ int i,j,k;+ int tot = cr - n;+ glp_add_rows(lp, m);+ glp_add_cols(lp, n);++ //printf("%d %d\n",m,n);++ // the first n values+ for (k=1;k<=n;k++) {+ glp_set_obj_coef(lp, k, AT(c, k-1, 2));+ //printf("%d %f\n",k,AT(c, k-1, 2));+ }++ int * ia = malloc((1+tot)*sizeof(int));+ int * ja = malloc((1+tot)*sizeof(int));+ double * ar = malloc((1+tot)*sizeof(double));++ for (k=1; k<= tot; k++) {+ ia[k] = rint(AT(c,k-1+n,0));+ ja[k] = rint(AT(c,k-1+n,1));+ ar[k] = AT(c,k-1+n,2);+ //printf("%d %d %f\n",ia[k],ja[k],ar[k]);+ }+ glp_load_matrix(lp, tot, ia, ja, ar);++ int t;+ for (i=1;i<=m;i++) {+ switch((int)rint(AT(b,i-1,0))) {+ case 0: { t = GLP_FR; break; }+ case 1: { t = GLP_LO; break; }+ case 2: { t = GLP_UP; break; }+ case 3: { t = GLP_DB; break; }+ default: { t = GLP_FX; break; }+ }+ glp_set_row_bnds(lp, i, t , AT(b,i-1,1), AT(b,i-1,2));+ }+ for (j=1;j<=n;j++) {+ switch((int)rint(AT(b,m+j-1,0))) {+ case 0: { t = GLP_FR; break; }+ case 1: { t = GLP_LO; break; }+ case 2: { t = GLP_UP; break; }+ case 3: { t = GLP_DB; break; }+ default: { t = GLP_FX; break; }+ }+ glp_set_col_bnds(lp, j, t , AT(b,m+j-1,1), AT(b,m+j-1,2));+ }+ glp_term_out(0);+ glp_simplex(lp, NULL);+ sp[0] = glp_get_status(lp);+ sp[1] = glp_get_obj_val(lp);+ for (k=1; k<=n; k++) {+ sp[k+1] = glp_get_col_prim(lp, k);+ }+ glp_delete_prob(lp);+ free(ia);+ free(ja);+ free(ar);++ return 0;+}