hmatrix-glpk 0.4.0.2 → 0.4.1.0
raw patch · 8 files changed
+178/−76 lines, 8 filesdep +containersPVP ok
version bump matches the API change (PVP)
Dependencies added: containers
API changes (from Hackage documentation)
+ Numeric.LinearProgramming: General :: [Bound [(Double, Int)]] -> Constraints
+ Numeric.LinearProgramming: exact :: Optimization -> Constraints -> Bounds -> Solution
+ Numeric.LinearProgramming: sparseOfGeneral :: Constraints -> Constraints
Files
- examples/simplex1.hs +3/−3
- examples/simplex2.hs +6/−1
- examples/simplex3.hs +1/−1
- examples/simplex4.hs +1/−1
- examples/simplex5.hs +27/−0
- hmatrix-glpk.cabal +3/−2
- src/C/glpk.c +67/−63
- src/Numeric/LinearProgramming.hs +70/−5
examples/simplex1.hs view
@@ -9,9 +9,9 @@ , [2,2,6] :<=: 300 ] -- default bounds-bnds = [ 1 :=>: 0- , 2 :=>: 0- , 3 :=>: 0 ]+bnds = [ 1 :>=: 0+ , 2 :>=: 0+ , 3 :>=: 0 ] main = do print $ simplex objFun constr []
examples/simplex2.hs view
@@ -10,9 +10,14 @@ , [0,1,5] :<=: 20 ] +constr3 = General [ [1#1, 1#1, 1#2] :<=: 10+ , [1#2, 5#3] :<=: 20+ ]+ main = do print $ simplex prob constr1 [] print $ simplex prob constr2 []- print $ simplex prob constr2 [ 2 :=>: 1, 3 :&: (2,7)]+ print $ simplex prob constr3 []+ print $ simplex prob constr2 [ 2 :>=: 1, 3 :&: (2,7)] print $ simplex prob constr2 [ Free 2 ]
examples/simplex3.hs view
@@ -11,7 +11,7 @@ , [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.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) ]
examples/simplex4.hs view
@@ -11,7 +11,7 @@ , [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.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) ]
+ examples/simplex5.hs view
@@ -0,0 +1,27 @@+import Numeric.LinearProgramming++-- This is a linear program from the paper "Picking vs. Guessing Secrets: A Game-theoretic Analysis"++gamma = 100000 :: Double+sigma = 1 :: Double+n = 64 :: Int+cost_fun :: Int -> Double+cost_fun i = (fromIntegral i) / (fromIntegral n)+size_fun :: Int -> Double+size_fun i = 2^(fromIntegral i)++prob = Minimize $ map cost_fun [1..n]+bnds = [i :&: (0,1) | i <- [1..n]]++constr1 = [[1 # i | i <- [1..n]] :==: 1] +++ [[1/(size_fun i) # i,+ -1/(size_fun (i+1)) # i+1] :>=: 0 | i <- [1..n-1]] +++ [(+ [gamma#i | i <- [1..k]] +++ (concat [[sigma*(size_fun i) # j | j <- [1..i-1]] | i <- [1..k]]) +++ [((size_fun i) - 1)/2 # i | i <- [1..k]])+ :<=: (sigma * (foldr (+) 0 (map size_fun [1..k]))) | k <- [1..n]]++main = do+ print $ simplex prob (General constr1) bnds -- NoFeasible+ print $ exact prob (General constr1) bnds -- solution found
hmatrix-glpk.cabal view
@@ -1,5 +1,5 @@ Name: hmatrix-glpk-Version: 0.4.0.2+Version: 0.4.1.0 License: GPL License-file: LICENSE Author: Alberto Ruiz@@ -20,9 +20,10 @@ examples/simplex2.hs examples/simplex3.hs examples/simplex4.hs+ examples/simplex5.hs library- Build-Depends: base <5, hmatrix >= 0.16+ Build-Depends: base <5, hmatrix >= 0.16, containers >= 0.5.4.0 hs-source-dirs: src
src/C/glpk.c view
@@ -10,67 +10,71 @@ /*-----------------------------------------------------*/ -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);+#define C_X_SPARSE(X) \+ int c_##X##_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_##X(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; \+ } \ - return 0;-}+C_X_SPARSE(simplex);+C_X_SPARSE(exact);
src/Numeric/LinearProgramming.hs view
@@ -49,6 +49,14 @@ ] @ +Note 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.++@+constr3 = General [ [1\#1, 1\#1, 1\#2] :<=: 10+ , [1\#2, 5\#3] :<=: 20+ ]+@+ By default all variables are bounded as @x_i >= 0@, but this can be changed: @@ -67,6 +75,8 @@ module Numeric.LinearProgramming( simplex,+ exact,+ sparseOfGeneral, Optimization(..), Constraints(..), Bounds,@@ -82,13 +92,14 @@ import Foreign.C.Types import Data.List((\\),sortBy,nub) import Data.Function(on)+import qualified Data.Map.Strict as Map --import Debug.Trace --debug x = trace (show x) x ----------------------------------------------------- --- | Coefficient of a variable for a sparse representation of constraints.+-- | Coefficient of a variable for a sparse and general representations of constraints. (#) :: Double -> Int -> (Double,Int) infixl 5 # (#) = (,)@@ -108,18 +119,29 @@ | Unbounded deriving Show -data Constraints = Dense [ Bound [Double] ]- | Sparse [ Bound [(Double,Int)] ]+data Constraints = Dense [ Bound [Double] ]+ | Sparse [ Bound [(Double,Int)] ]+ | General [ Bound [(Double,Int)] ] data Optimization = Maximize [Double] | Minimize [Double] type Bounds = [Bound Int] +-- | Convert a system of General constraints to one with unique coefficients. +sparseOfGeneral :: Constraints -> Constraints+sparseOfGeneral (General cs) =+ Sparse $ map (\bl -> + let cl = obj bl in+ let cl_unique = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in+ withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] cl_unique)) cs+sparseOfGeneral _ = error "sparseOfGeneral: a general system of constraints was expected"+ 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 []) bnds = simplex opt (Sparse []) bnds+simplex opt (Sparse []) bnds = simplex opt (Sparse [Free [0#1]]) bnds+simplex opt (General []) 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)@@ -133,6 +155,29 @@ m = length constr (sz, sg, objfun) = adapt opt +simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds++-- | Simplex method with exact internal arithmetic. See glp_exact in glpk documentation for more information.+exact :: Optimization -> Constraints -> Bounds -> Solution++exact opt (Dense []) bnds = exact opt (Sparse []) bnds+exact opt (Sparse []) bnds = exact opt (Sparse [Free [0#1]]) bnds+exact opt (General []) bnds = exact opt (Sparse [Free [0#1]]) bnds++exact opt (Dense constr) bnds = extract sg sol where+ sol = exactSparse m n (mkConstrD sz objfun constr) (mkBounds sz constr bnds)+ n = length objfun+ m = length constr+ (sz, sg, objfun) = adapt opt++exact opt (Sparse constr) bnds = extract sg sol where+ sol = exactSparse m n (mkConstrS sz objfun constr) (mkBounds sz constr bnds)+ n = length objfun+ m = length constr+ (sz, sg, objfun) = adapt opt++exact opt constr@(General _) bnds = exact opt (sparseOfGeneral constr) bnds+ adapt :: Optimization -> (Int, Double, [Double]) adapt opt = case opt of Maximize x -> (size x, 1 ,x)@@ -162,6 +207,13 @@ obj (x :==: _) = x obj (Free x) = x +withObj :: Bound t -> t -> Bound t+withObj (_ :<=: b) x = (x :<=: b)+withObj (_ :>=: b) x = (x :>=: b)+withObj (_ :&: b) x = (x :&: b)+withObj (_ :==: b) x = (x :==: b)+withObj (Free _) x = Free x+ tb :: Bound t -> Double tb (_ :<=: _) = glpUP tb (_ :>=: _) = glpLO@@ -233,6 +285,19 @@ simplexSparse m n c b = unsafePerformIO $ do s <- createVector (2+n) app3 (c_simplex_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_simplex_sparse"+ return s++foreign import ccall unsafe "c_exact_sparse" c_exact_sparse+ :: CInt -> CInt -- rows and cols+ -> CInt -> CInt -> Ptr Double -- coeffs+ -> CInt -> CInt -> Ptr Double -- bounds+ -> CInt -> Ptr Double -- result+ -> IO CInt -- exit code++exactSparse :: Int -> Int -> Matrix Double -> Matrix Double -> Vector Double+exactSparse m n c b = unsafePerformIO $ do+ s <- createVector (2+n)+ app3 (c_exact_sparse (fi m) (fi n)) mat (cmat c) mat (cmat b) vec s "c_exact_sparse" return s glpFR, glpLO, glpUP, glpDB, glpFX :: Double