From e57907c22e8a16d2a9b62b70dd04ffcfd0d96b6a Mon Sep 17 00:00:00 2001 From: Piotr Mardziel Date: Sun, 22 Feb 2015 22:58:30 -0500 Subject: added handling of general sparse constraints --- packages/glpk/src/Numeric/LinearProgramming.hs | 40 ++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 5 deletions(-) (limited to 'packages/glpk/src/Numeric') diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs index b0537cc..a54eb59 100644 --- a/packages/glpk/src/Numeric/LinearProgramming.hs +++ b/packages/glpk/src/Numeric/LinearProgramming.hs @@ -49,6 +49,14 @@ constr2 = Dense [ [2,1,0] :<=: 10 ] @ +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,7 @@ Multiple bounds for a variable are not allowed, instead of module Numeric.LinearProgramming( simplex, + sparseOfGeneral, Optimization(..), Constraints(..), Bounds, @@ -82,13 +91,14 @@ import System.IO.Unsafe(unsafePerformIO) 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 +118,29 @@ data Solution = Undefined | 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 m = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in + withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] m)) 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 +154,8 @@ simplex opt (Sparse constr) bnds = extract sg sol where m = length constr (sz, sg, objfun) = adapt opt +simplex opt constr@(General _) bnds = simplex opt (sparseOfGeneral constr) bnds + adapt :: Optimization -> (Int, Double, [Double]) adapt opt = case opt of Maximize x -> (size x, 1 ,x) @@ -162,6 +185,13 @@ obj (x :&: _) = x 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 -- cgit v1.2.3 From 456aa8ebb8f8ab67f526b33930a54769b15c138a Mon Sep 17 00:00:00 2001 From: Piotr Mardziel Date: Mon, 23 Feb 2015 17:43:16 -0500 Subject: the rest of the files for glp_exact --- packages/glpk/hmatrix-glpk.cabal | 1 + packages/glpk/src/C/glpk.c | 130 +++++++++++++------------ packages/glpk/src/Numeric/LinearProgramming.hs | 39 +++++++- 3 files changed, 105 insertions(+), 65 deletions(-) (limited to 'packages/glpk/src/Numeric') diff --git a/packages/glpk/hmatrix-glpk.cabal b/packages/glpk/hmatrix-glpk.cabal index 229197f..a9859f9 100644 --- a/packages/glpk/hmatrix-glpk.cabal +++ b/packages/glpk/hmatrix-glpk.cabal @@ -20,6 +20,7 @@ extra-source-files: examples/simplex1.hs examples/simplex2.hs examples/simplex3.hs examples/simplex4.hs + examples/simplex5.hs library Build-Depends: base <5, hmatrix >= 0.16, containers >= 0.5.4.0 diff --git a/packages/glpk/src/C/glpk.c b/packages/glpk/src/C/glpk.c index bfbb435..86b1277 100644 --- a/packages/glpk/src/C/glpk.c +++ b/packages/glpk/src/C/glpk.c @@ -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); +#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; \ + } \ - //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; -} +C_X_SPARSE(simplex); +C_X_SPARSE(exact); diff --git a/packages/glpk/src/Numeric/LinearProgramming.hs b/packages/glpk/src/Numeric/LinearProgramming.hs index a54eb59..6a0c47d 100644 --- a/packages/glpk/src/Numeric/LinearProgramming.hs +++ b/packages/glpk/src/Numeric/LinearProgramming.hs @@ -75,6 +75,7 @@ Multiple bounds for a variable are not allowed, instead of module Numeric.LinearProgramming( simplex, + exact, sparseOfGeneral, Optimization(..), Constraints(..), @@ -132,8 +133,8 @@ sparseOfGeneral :: Constraints -> Constraints sparseOfGeneral (General cs) = Sparse $ map (\bl -> let cl = obj bl in - let m = foldr (\(c,t) m -> Map.insertWith (+) t c m) Map.empty cl in - withObj bl (Map.foldrWithKey' (\t c l -> (c#t) : l) [] m)) cs + 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 @@ -156,6 +157,27 @@ simplex opt (Sparse constr) bnds = extract sg sol where 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) @@ -265,6 +287,19 @@ simplexSparse m n c b = unsafePerformIO $ do 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 glpFR = 0 glpLO = 1 -- cgit v1.2.3