comfort-glpk 0.0.0.1 → 0.0.1

raw patch · 8 files changed

+822/−146 lines, 8 filesdep +randomdep +transformersdep ~non-emptydep ~utility-htPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: random, transformers

Dependency ranges changed: non-empty, utility-ht

API changes (from Hackage documentation)

Files

• comfort-glpk.cabal +12/−4
• src/Numeric/GLPK.hs +227/−128
• src/Numeric/GLPK/Monad.hs +110/−0
• src/Numeric/GLPK/Private.hs +177/−0
• test/Main.hs +2/−0
• test/Test/Numeric/GLPK.hs +60/−14
• test/Test/Numeric/GLPK/Generator.hs +208/−0
• test/Test/Numeric/GLPK/Monad.hs +26/−0

@@ -1,6 +1,6 @@ Cabal-Version:    2.2 Name:             comfort-glpk-Version:          0.0.0.1+Version:          0.0.1 License:          BSD-3-Clause License-File:     LICENSE Author:           Henning Thielemann <haskell@henning-thielemann.de>@@ -30,10 +30,10 @@   >           printf "absol %f,  pos %f, neg %f\n" absol pos neg   >         _ -> fail "GLPK solver failed"   .-  Alternatives: @hmatrix-glpk@, @glpk-hs@+  Alternatives: @coinor-clp@, @hmatrix-glpk@, @glpk-hs@  Source-Repository this-  Tag:         0.0.0.1+  Tag:         0.0.1   Type:        darcs   Location:    https://hub.darcs.net/thielema/comfort-glpk/ @@ -46,6 +46,7 @@     glpk-headers >=0.4.1 && <0.6,     comfort-array >=0.4 && <0.6,     deepseq >=1.3 && <1.5,+    transformers >=0.3 && <0.7,     non-empty >=0.3.2 && <0.4,     utility-ht >=0.0.16 && <0.1,     base >=4.5 && <5@@ -56,16 +57,21 @@   Extra-Libraries: glpk   Exposed-modules:     Numeric.GLPK+    Numeric.GLPK.Monad+  Other-Modules:+    Numeric.GLPK.Private  Test-Suite comfort-glpk-test   Type:             exitcode-stdio-1.0   Build-Depends:     comfort-glpk,     comfort-array >=0.5.2,-    utility-ht,+    non-empty,+    utility-ht >=0.0.17,     doctest-exitcode-stdio >=0.0 && <0.1,     doctest-lib >=0.1 && <0.2,     QuickCheck >=2.1 && <3,+    random >=1.0 && <1.3,     base >=4.5 && <5    GHC-Options:      -Wall@@ -73,4 +79,6 @@   Default-Language: Haskell98   Main-Is: Main.hs   Other-Modules:+    Test.Numeric.GLPK.Generator+    Test.Numeric.GLPK.Monad     Test.Numeric.GLPK

@@ -38,34 +38,49 @@    objectiveFromTerms,    simplex,    simplexMulti,+   simplexSuccessive,    exact,+   exactMulti,+   exactSuccessive,    interior,    interiorMulti,+   interiorSuccessive,++   solveSuccessive,++   FormatIdentifier,+   formatMathProg,    ) where  import qualified Math.Programming.Glpk.Header as FFI+import Numeric.GLPK.Private  import qualified Data.Array.Comfort.Storable.Mutable as Mutable import qualified Data.Array.Comfort.Storable as Array import qualified Data.Array.Comfort.Shape as Shape import qualified Data.NonEmpty as NonEmpty+import qualified Data.List as List import Data.Array.Comfort.Storable (Array) import Data.Tuple.HT (mapFst, mapSnd)-import Data.Maybe (fromMaybe)-import Data.Foldable (traverse_, for_)+import Data.Traversable (for)+import Data.Foldable (for_) -import Control.Monad (void)+import Text.Printf (printf)++import qualified Control.Monad.Trans.Except as ME+import qualified Control.Monad.Trans.State as MS+import Control.Monad (void, when) import Control.Applicative (liftA2) import Control.Exception (bracket)-import Control.DeepSeq (NFData, rnf)  import System.IO.Unsafe (unsafePerformIO)  import qualified Foreign import Foreign.Ptr (nullPtr)-import Foreign.C.Types (CDouble) + {- $setup+>>> import qualified Test.Numeric.GLPK.Generator as TestLP >>> import qualified Numeric.GLPK as LP >>> import Numeric.GLPK ((.*), (<=.), (==.)) >>>@@ -97,10 +112,7 @@ (.*) :: Double -> ix -> Term ix (.*) = Term -data Term ix = Term Double ix-   deriving (Show) - infix 4 <=., >=., >=<., ==.  (<=.), (>=.), (==.) :: x -> Double -> Inequality x@@ -114,88 +126,21 @@ free :: x -> Inequality x free x = Inequality x Free -data Inequality x = Inequality x Bound-   deriving Show -data Bound =-     LessEqual Double-   | GreaterEqual Double-   | Between Double Double-   | Equal Double-   | Free-   deriving Show -instance Functor Inequality where-   fmap f (Inequality x bnd)  =  Inequality (f x) bnd--data NoSolutionType =-     Undefined-   | NoFeasible-   | Unbounded-   deriving (Eq, Show)--data SolutionType =-     Feasible-   | Infeasible-   | Optimal-   deriving (Eq, Show)--instance NFData NoSolutionType where-    rnf NoFeasible = ()-    rnf _ = ()--instance NFData SolutionType where-    rnf Optimal = ()-    rnf _ = ()--type Solution sh =-      Either NoSolutionType (SolutionType, (Double, Array sh Double))--type Constraints ix = [Inequality [Term ix]]--data Direction = Minimize | Maximize--type Objective sh = Array sh Double--type Bounds ix = [Inequality ix]--- objectiveFromTerms ::    (Shape.Indexed sh, Shape.Index sh ~ ix) => sh -> [Term ix] -> Objective sh objectiveFromTerms sh =    Array.fromAssociations 0 sh . map (\(Term x ix) -> (ix,x))  -prepareBounds :: Inequality a -> (a, (FFI.GlpkConstraintType, CDouble, CDouble))-prepareBounds (Inequality x bnd) =-   (,) x $-   (\(bndType,lo,up) -> (bndType, realToFrac lo, realToFrac up)) $-   case bnd of-      LessEqual up    -> (FFI.glpkLT,      0,  up)-      GreaterEqual lo -> (FFI.glpkGT,      lo, 0)-      Between lo up   -> (FFI.glpkBounded, lo, up)-      Equal y         -> (FFI.glpkFixed,   y,  y)-      Free            -> (FFI.glpkFree,    0,  0)--columnIndex :: (Shape.Indexed sh, Shape.Index sh ~ ix) => sh -> ix -> FFI.Column-columnIndex shape var = 1 + fromIntegral (Shape.offset shape var)--deferredColumnIndex :: Shape.DeferredIndex ix -> FFI.Column-deferredColumnIndex (Shape.DeferredIndex var) = 1 + fromIntegral var--allocaArray :: (Foreign.Storable a) => Int -> (FFI.GlpkArray a -> IO b) -> IO b-allocaArray n f = Foreign.allocaArray (n+1) $ f . FFI.GlpkArray--pokeElem :: (Foreign.Storable a) => FFI.GlpkArray a -> Int -> a -> IO ()-pokeElem (FFI.GlpkArray ptr) k a = Foreign.pokeElemOff ptr k a-- {- | >>> case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double) Right (Optimal,(28.0,(5.0,0.0,4.0))) +prop> \target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (1,1) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (TestLP.approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False prop> \(QC.Positive posWeight) (QC.Positive negWeight) target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (posWeight,negWeight) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False+prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, _) -> True; _ -> False -} simplex ::    (Shape.Indexed sh, Shape.Index sh ~ ix) =>@@ -203,6 +148,9 @@    (Direction, Objective sh) -> Solution sh simplex = solve (flip FFI.glp_simplex nullPtr) +{-# DEPRECATED simplexMulti "use GLPK.Monad instead" #-}+{-# DEPRECATED exactMulti "use GLPK.Monad instead" #-}+{-# DEPRECATED interiorMulti "run 'interior' in Either monad instead" #-} {- | Optimize for one objective after another. That is, if the first optimization succeeds@@ -227,12 +175,19 @@ because an added objective can make the system infeasible due to rounding errors. E.g. a non-negative objective can become very small but negative.+++prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case LP.simplexMulti bounds constrs (Array.shape origin) objs of (_, Right (LP.Optimal, _)) -> QC.property True; result -> QC.counterexample (show result) False++The same property fails for 'exactMulti' and 'interiorMulti'.+I guess, due to rounding errors. -}-simplexMulti, interiorMulti ::+simplexMulti, exactMulti, interiorMulti ::    (Shape.Indexed sh, Shape.Index sh ~ ix) =>    Bounds ix -> Constraints ix ->    sh -> NonEmpty.T [] (Direction, [Term ix]) -> ([Double], Solution sh) simplexMulti = solveMulti . simplex+exactMulti = solveMulti . exact interiorMulti = solveMulti . interior  solveMulti ::@@ -245,16 +200,78 @@          mapFst (opt:) $          let extConstrs = (curObj==.opt) : constrs in          go extConstrs obj objs $-         solver extConstrs-            (dir, objectiveFromTerms sh obj)+         solver extConstrs (dir, objectiveFromTerms sh obj)        go _ _ _ sol = ([], sol)    in go constrs0 (snd obj0) objs0 $       solver constrs0 $ mapSnd (objectiveFromTerms sh) obj0  +{-# DEPRECATED simplexSuccessive "use GLPK.Monad instead" #-}+{-# DEPRECATED exactSuccessive "use GLPK.Monad instead" #-}+{-# DEPRECATED interiorSuccessive "run 'interior' in Either monad instead" #-} {- |+Like the @Multi@ functions,+but allows not only to fix the previously+found optimal solution as constraint,+but allows constraints with a tolerance.+This is necessary in the presence of rounding errors.++prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.simplexSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False+prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.exactSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False+-}+simplexSuccessive, exactSuccessive, interiorSuccessive ::+   (Traversable f, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) ->+   f ((SolutionType, (Double, Array sh Double)) -> Constraints ix,+      (Direction, Objective sh)) ->+   Either NoSolutionType+      (NonEmpty.T f (SolutionType, (Double, Array sh Double)))+simplexSuccessive = solveSuccessiveInPlace (flip FFI.glp_simplex nullPtr)+exactSuccessive = solveSuccessiveInPlace (flip FFI.glp_exact nullPtr)+interiorSuccessive = solveSuccessive . interior++{-# DEPRECATED solveSuccessive+      "run simple solvers in GLPK.Monad or Either monad instead" #-}+{- |+Allows for generic implementation of 'simplexSuccessive' et.al.+without reuse of interim results.++prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.simplexSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.simplex bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB+prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.exactSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.exact bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB+-}+solveSuccessive ::+   (Traversable f, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   (Constraints ix -> (Direction, Objective sh) -> Solution sh) ->+   Constraints ix ->+   (Direction, Objective sh) ->+   f ((SolutionType, (Double, Array sh Double)) -> Constraints ix,+      (Direction, Objective sh)) ->+   Either NoSolutionType+      (NonEmpty.T f (SolutionType, (Double, Array sh Double)))+solveSuccessive solver constrs0 obj0 objs = do+   let checkShape obj =+         if Array.shape (snd obj0) == Array.shape obj+            then obj+            else error "GLPK.solveSuccessive: objective shapes mismatch"+   let solveWithConstraints constrs problem =+         (\sol -> (sol, (constrs,sol))) <$> solver constrs problem+   (sol0,state0) <- solveWithConstraints constrs0 obj0+   NonEmpty.cons sol0 <$>+      MS.evalStateT+         (for objs $+            \(newConstrs,(dir,obj)) -> MS.StateT $ \(constrs,sol) ->+               solveWithConstraints+                  (newConstrs sol ++ constrs)+                  (dir, checkShape obj))+         state0+++{- | >>> case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.exact [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double) Right (Optimal,(28.0,(5.0,0.0,4.0)))++prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.exact bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False -} exact ::    (Shape.Indexed sh, Shape.Index sh ~ ix) =>@@ -273,20 +290,60 @@    bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> do    storeProblem bounds constrs (dir,obj) lp    void $ solver lp-   let examine =-         liftA2 (,)-            (realToFrac <$> FFI.glp_get_obj_val lp)-            (readGLPArray (Array.shape obj) $ \arr ix ->-               Mutable.write arr ix . realToFrac-                  =<< FFI.glp_get_col_prim lp (deferredColumnIndex ix))-   status <- FFI.glp_get_status lp-   either (return . Left) (\typ -> Right . (,) typ <$> examine) $-      analyzeStatus status+   peekSimplexSolution (Array.shape obj) lp +{-# INLINE solveSuccessiveInPlace #-}+solveSuccessiveInPlace ::+   (Traversable f, Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   (Foreign.Ptr FFI.Problem -> IO FFI.GlpkSimplexStatus) ->+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) ->+   f ((SolutionType, (Double, Array sh Double)) -> Constraints ix,+      (Direction, Objective sh)) ->+   Either NoSolutionType+      (NonEmpty.T f (SolutionType, (Double, Array sh Double)))+solveSuccessiveInPlace solver bounds constrs0 (dir0,obj0) objs =+      unsafePerformIO $+   bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> ME.runExceptT $ do+   let shape = Array.shape obj0+   sol0 <- ME.ExceptT $ do+      storeProblem bounds constrs0 (dir0,obj0) lp+      void $ solver lp+      peekSimplexSolution shape lp+   NonEmpty.cons sol0 <$>+      MS.evalStateT+         (for objs $+            \(makeNewConstrs,(dir,obj)) -> MS.StateT $ \sol ->+                  fmap (\sol1 -> (sol1, sol1)) $ ME.ExceptT $ do+               setDirection lp dir+               when (shape /= Array.shape obj) $+                  error "GLPK.solveSuccessiveInplace: objective shapes mismatch"+               setObjective lp obj+               let newConstrs = makeNewConstrs sol+               newRow <- FFI.glp_add_rows lp $ fromIntegral $ length newConstrs+               for_ (zip [newRow..] (map prepareBounds newConstrs)) $+                     \(row, (terms,(bnd,lo,up))) -> do+                  FFI.glp_set_row_bnds lp row bnd lo up+                  let numTerms = length terms+                  allocaArray numTerms $ \indicesPtr ->+                     allocaArray numTerms $ \coeffsPtr -> do+                     for_ (zip [1..] terms) $+                        \(k, Term c x) -> do+                           pokeElem indicesPtr k (columnIndex shape x)+                           pokeElem coeffsPtr k (realToFrac c)+                     FFI.glp_set_mat_row lp row+                        (fromIntegral numTerms) indicesPtr coeffsPtr+               void $ solver lp+               peekSimplexSolution shape lp)+         sol0 ++ {- | >>> case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapPair (round3, Array.toTuple . Array.map round3)) <$> LP.interior [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double) Right (Optimal,(28.0,(5.0,0.0,4.0)))++prop> QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.interior bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False -} interior ::    (Shape.Indexed sh, Shape.Index sh ~ ix) =>@@ -312,27 +369,19 @@    Bounds ix -> Constraints ix ->    (Direction, Objective sh) -> Foreign.Ptr FFI.Problem -> IO () storeProblem bounds constrs (dir,obj) lp = do-   let shape = Array.shape obj    void $ FFI.glp_term_out FFI.glpkOff-   FFI.glp_set_obj_dir lp $-      case dir of-         Minimize -> FFI.glpkMin-         Maximize -> FFI.glpkMax+   let shape = Array.shape obj+   setDirection lp dir    firstRow <- FFI.glp_add_rows lp $ fromIntegral $ length constrs    for_ (zip [firstRow..] $       map prepareBounds constrs) $ \(row,(_x,(bnd,lo,up))) ->       FFI.glp_set_row_bnds lp row bnd lo up-   _firstCol <- FFI.glp_add_cols lp $ fromIntegral $ Shape.size shape-   for_ (Shape.indices $ Shape.Deferred shape) $ \x ->-      FFI.glp_set_col_bnds lp (deferredColumnIndex x) FFI.glpkGT 0 0-   for_ (map prepareBounds bounds) $ \(x,(bnd,lo,up)) ->-      FFI.glp_set_col_bnds lp (columnIndex shape x) bnd lo up-   for_ (Array.toAssociations obj) $ \(x,c) ->-      FFI.glp_set_obj_coef lp (columnIndex shape x) (realToFrac c)-   let numRows = length $ concatMap (fst . prepareBounds) constrs-   allocaArray numRows $ \ia ->-      allocaArray numRows $ \ja ->-      allocaArray numRows $ \ar -> do+   storeBounds lp shape bounds+   setObjective lp obj+   let numTerms = length $ concatMap (fst . prepareBounds) constrs+   allocaArray numTerms $ \ia ->+      allocaArray numTerms $ \ja ->+      allocaArray numTerms $ \ar -> do       for_ (zip [1..] $ concat $             zipWith (map . (,)) [firstRow..] $             map (fst . prepareBounds) constrs) $@@ -340,28 +389,78 @@             pokeElem ia k row             pokeElem ja k (columnIndex shape x)             pokeElem ar k (realToFrac c)-      FFI.glp_load_matrix lp (fromIntegral numRows) ia ja ar+      FFI.glp_load_matrix lp (fromIntegral numTerms) ia ja ar -{-# INLINE readGLPArray #-}-readGLPArray ::-   (Shape.C sh, Foreign.Storable a, Num a) =>-   sh ->-   (Mutable.Array IO (Shape.Deferred sh) a ->-    Shape.DeferredIndex sh -> IO ()) ->-   IO (Array sh a)-readGLPArray shape act = do-   let defShape = Shape.Deferred shape-   arr <- Mutable.new defShape 0-   traverse_ (act arr) (Shape.indices defShape)-   Array.reshape shape <$> Mutable.freeze arr -analyzeStatus :: FFI.GlpkSolutionStatus -> Either NoSolutionType SolutionType-analyzeStatus status =-   fromMaybe (error "glpk-solver: solution type unknown") $ lookup status $-      (FFI.glpkUndefined,  Left Undefined) :-      (FFI.glpkFeasible,   Right Feasible) :-      (FFI.glpkInfeasible, Right Infeasible) :-      (FFI.glpkNoFeasible, Left NoFeasible) :-      (FFI.glpkOptimal,    Right Optimal) :-      (FFI.glpkUnbounded,  Left Unbounded) :-      []+++class FormatIdentifier ix where+   formatIdentifier :: ix -> String++instance FormatIdentifier Char where+   formatIdentifier x = [x]++instance FormatIdentifier c => FormatIdentifier [c] where+   formatIdentifier = concatMap formatIdentifier++instance FormatIdentifier Int where+   formatIdentifier = printf "x%d"++instance FormatIdentifier Integer where+   formatIdentifier = printf "x%d"+++formatBound :: (FormatIdentifier ix) => Inequality ix -> String+formatBound (Inequality ix bnd) =+   printf "var %s%s;" (formatIdentifier ix) $+   case bnd of+      LessEqual up -> printf ", <=%f" up+      GreaterEqual lo -> printf ", >=%f" lo+      Between lo up -> printf ", >=%f, <=%f" lo up+      Equal x -> printf ", =%f" x+      Free -> ""+++formatSum :: (FormatIdentifier ix) => [Term ix] -> String+formatSum [] = "0"+formatSum xs =+   let formatTerm (Term c ix) = printf "%f*%s" c (formatIdentifier ix) in+   List.intercalate "+" $ map formatTerm xs++formatConstraint :: (FormatIdentifier ix) => Inequality [Term ix] -> String+formatConstraint (Inequality terms bnd) =+   let sumStr = formatSum terms in+   case bnd of+      LessEqual up -> printf "%s <= %f" sumStr up+      GreaterEqual lo -> printf "%f <= %s" lo sumStr+      Between lo up -> printf "%f <= %s <= %f" lo sumStr up+      Equal x -> printf "%s = %f" sumStr x+      Free -> sumStr++formatDirection :: Direction -> String+formatDirection Minimize = "minimize"+formatDirection Maximize = "maximize"++formatObjective ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, FormatIdentifier ix) =>+   Objective sh -> String+formatObjective =+   formatSum . map (\(ix,c) -> Term c ix) . Array.toAssociations++formatMathProg ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, FormatIdentifier ix) =>+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) -> [String]+formatMathProg bounds constrs (dir,obj) =+   map formatBound bounds +++   "" :+   formatDirection dir :+   printf "value: %s;" (formatObjective obj) :+   "" :+   "subject to" :+   zipWith+      (\k constr -> printf "constr%d: %s;" k $ formatConstraint constr)+      [(0::Int)..] constrs +++   "" :+   "end;" :+   []

@@ -0,0 +1,110 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{- |+The monadic interface to GLPK allows to optimize+with respect to multiple objectives, successively.+-}+module Numeric.GLPK.Monad (+   T,+   run,+   simplex,+   exact,+   Direction(..),+   ) where++import qualified Math.Programming.Glpk.Header as FFI+import Numeric.GLPK.Private+         (Term(Term), Constraints, Solution,+          allocaArray, pokeElem, columnIndex, prepareBounds, storeBounds,+          setDirection, setObjective, peekSimplexSolution)+import Numeric.GLPK (Bounds, Direction(..), Objective)++import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Foldable (for_)++import qualified Control.Monad.Trans.Reader as MR+import Control.Monad.IO.Class (liftIO)+import Control.Monad (void, when)+import Control.Exception (bracket)++import System.IO.Unsafe (unsafePerformIO)++import Foreign.Ptr (Ptr, nullPtr)+++{- $setup+>>> import qualified Numeric.GLPK.Monad as LP+>>> import Test.Numeric.GLPK (TripletShape, tripletShape)+>>> import Numeric.GLPK ((.*), (<=.))+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import qualified Data.Array.Comfort.Shape as Shape+>>>+>>> import Data.Tuple.HT (mapSnd)+-}+++newtype T sh a = Cons (MR.ReaderT (sh, Ptr FFI.Problem) IO a)+   deriving (Functor, Applicative, Monad)++run ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   sh -> Bounds ix -> T sh a -> a+run shape bounds (Cons act) =+   unsafePerformIO $+   bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> do+      void $ FFI.glp_term_out FFI.glpkOff+      storeBounds lp shape bounds+      liftIO $ MR.runReaderT act (shape, lp)+++{- |+Add new constraints to an existing problem+and solve with a new direction and objective.++>>> case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd (mapSnd Array.toTuple) <$> LP.run tripletShape [] (LP.simplex [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))+Right (Optimal,(28.0,(5.0,0.0,4.0)))+-}+simplex ::+   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Constraints ix ->+   (Direction, Objective sh) -> T sh (Solution sh)+simplex = solve (flip FFI.glp_simplex nullPtr)++exact ::+   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Constraints ix ->+   (Direction, Objective sh) -> T sh (Solution sh)+exact = solve (flip FFI.glp_exact nullPtr)++solve ::+   (Eq sh, Shape.Indexed sh, Shape.Index sh ~ ix) =>+   (Ptr FFI.Problem -> IO FFI.GlpkSimplexStatus) ->+   Constraints ix ->+   (Direction, Objective sh) ->+   T sh (Solution sh)+solve method constrs (dir,obj) = Cons $ do+   (shape, lp) <- MR.ask+   when (shape /= Array.shape obj) $+      error "GLPK.Monad.solve: objective shape mismatch"++   liftIO $ do+      setDirection lp dir+      setObjective lp obj+      newRow <- FFI.glp_add_rows lp $ fromIntegral $ length constrs+      for_ (zip [newRow..] (map prepareBounds constrs)) $+            \(row, (terms,(bnd,lo,up))) -> do+         FFI.glp_set_row_bnds lp row bnd lo up+         let numTerms = length terms+         allocaArray numTerms $ \indicesPtr ->+            allocaArray numTerms $ \coeffsPtr -> do+            for_ (zip [1..] terms) $+               \(k, Term c x) -> do+                  pokeElem indicesPtr k (columnIndex shape x)+                  pokeElem coeffsPtr k (realToFrac c)+            FFI.glp_set_mat_row lp row+               (fromIntegral numTerms) indicesPtr coeffsPtr+      void $ method lp+      peekSimplexSolution shape lp

@@ -0,0 +1,177 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+module Numeric.GLPK.Private where++import qualified Math.Programming.Glpk.Header as FFI++import qualified Data.Array.Comfort.Storable.Mutable as Mutable+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import Data.Array.Comfort.Storable (Array)+import Data.Maybe (fromMaybe)+import Data.Foldable (for_)++import Control.Applicative (liftA2)+import Control.DeepSeq (NFData, rnf)++import qualified Foreign+import Foreign.C.Types (CDouble)++++data Term ix = Term Double ix+   deriving (Show)+++data Inequality x = Inequality x Bound+   deriving Show++data Bound =+     LessEqual Double+   | GreaterEqual Double+   | Between Double Double+   | Equal Double+   | Free+   deriving Show++instance Functor Inequality where+   fmap f (Inequality x bnd)  =  Inequality (f x) bnd++type Bounds ix = [Inequality ix]++type Constraints ix = [Inequality [Term ix]]++data Direction = Minimize | Maximize+   deriving (Eq, Enum, Bounded, Show)++type Objective sh = Array sh Double+++data NoSolutionType =+     Undefined+   | NoFeasible+   | Unbounded+   deriving (Eq, Show)++data SolutionType =+     Feasible+   | Infeasible+   | Optimal+   deriving (Eq, Show)++instance NFData NoSolutionType where+    rnf NoFeasible = ()+    rnf _ = ()++instance NFData SolutionType where+    rnf Optimal = ()+    rnf _ = ()++type Solution sh =+      Either NoSolutionType (SolutionType, (Double, Array sh Double))+++{- |+@libglpk@ considers (Between x x) an error. @glpsol@ does not.+In handwritten problems, (Between x x) might indicate a mistake.+In automatically generated problems it will certainly not.+-}+canonicalizeBounds :: Inequality a -> Inequality a+canonicalizeBounds (Inequality x bnd) =+   Inequality x $+   case bnd of+      Between lo up -> if lo == up then Equal lo else bnd+      _ -> bnd++prepareBoundsFFI ::+   Inequality a -> (a, (FFI.GlpkConstraintType, CDouble, CDouble))+prepareBoundsFFI (Inequality x bnd) =+   (,) x $+   (\(bndType,lo,up) -> (bndType, realToFrac lo, realToFrac up)) $+   case bnd of+      LessEqual up    -> (FFI.glpkLT,      0,  up)+      GreaterEqual lo -> (FFI.glpkGT,      lo, 0)+      Between lo up   -> (FFI.glpkBounded, lo, up)+      Equal y         -> (FFI.glpkFixed,   y,  y)+      Free            -> (FFI.glpkFree,    0,  0)++prepareBounds ::+   Inequality a -> (a, (FFI.GlpkConstraintType, CDouble, CDouble))+prepareBounds = prepareBoundsFFI . canonicalizeBounds++storeBounds ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Foreign.Ptr FFI.Problem -> sh -> Bounds ix -> IO ()+storeBounds lp shape bounds = do+   _firstCol <- FFI.glp_add_cols lp $ fromIntegral $ Shape.size shape+   for_ (Shape.indices $ Shape.Deferred shape) $ \x ->+      FFI.glp_set_col_bnds lp (deferredColumnIndex x) FFI.glpkGT 0 0+   for_ (map prepareBounds bounds) $ \(x,(bnd,lo,up)) ->+      FFI.glp_set_col_bnds lp (columnIndex shape x) bnd lo up++++columnIndex :: (Shape.Indexed sh, Shape.Index sh ~ ix) => sh -> ix -> FFI.Column+columnIndex shape var = 1 + fromIntegral (Shape.offset shape var)++deferredColumnIndex :: Shape.DeferredIndex ix -> FFI.Column+deferredColumnIndex (Shape.DeferredIndex var) = 1 + fromIntegral var++allocaArray :: (Foreign.Storable a) => Int -> (FFI.GlpkArray a -> IO b) -> IO b+allocaArray n f = Foreign.allocaArray (n+1) $ f . FFI.GlpkArray++pokeElem :: (Foreign.Storable a) => FFI.GlpkArray a -> Int -> a -> IO ()+pokeElem (FFI.GlpkArray ptr) k a = Foreign.pokeElemOff ptr k a++++setDirection :: Foreign.Ptr FFI.Problem -> Direction -> IO ()+setDirection lp dir =+   FFI.glp_set_obj_dir lp $+      case dir of+         Minimize -> FFI.glpkMin+         Maximize -> FFI.glpkMax++setObjective ::+   (Shape.Indexed sh) => Foreign.Ptr FFI.Problem -> Objective sh -> IO ()+setObjective lp obj =+   for_ (Array.toAssociations obj) $ \(x,c) ->+      FFI.glp_set_obj_coef lp (columnIndex (Array.shape obj) x) (realToFrac c)++{-# INLINE readGLPArray #-}+readGLPArray ::+   (Shape.C sh, Foreign.Storable a, Num a) =>+   sh ->+   (Mutable.Array IO (Shape.Deferred sh) a ->+    Shape.DeferredIndex sh -> IO ()) ->+   IO (Array sh a)+readGLPArray shape act = do+   let defShape = Shape.Deferred shape+   arr <- Mutable.new defShape 0+   for_ (Shape.indices defShape) (act arr)+   Array.reshape shape <$> Mutable.freeze arr++analyzeStatus :: FFI.GlpkSolutionStatus -> Either NoSolutionType SolutionType+analyzeStatus status =+   fromMaybe (error "glpk-solver: solution type unknown") $ lookup status $+      (FFI.glpkUndefined,  Left Undefined) :+      (FFI.glpkFeasible,   Right Feasible) :+      (FFI.glpkInfeasible, Right Infeasible) :+      (FFI.glpkNoFeasible, Left NoFeasible) :+      (FFI.glpkOptimal,    Right Optimal) :+      (FFI.glpkUnbounded,  Left Unbounded) :+      []+++peekSimplexSolution ::+   (Shape.C sh) => sh -> Foreign.Ptr FFI.Problem -> IO (Solution sh)+peekSimplexSolution shape lp = do+   let examine =+         liftA2 (,)+            (realToFrac <$> FFI.glp_get_obj_val lp)+            (readGLPArray shape $ \arr ix ->+               Mutable.write arr ix . realToFrac+                  =<< FFI.glp_get_col_prim lp (deferredColumnIndex ix))+   status <- FFI.glp_get_status lp+   either (return . Left) (\typ -> Right . (,) typ <$> examine) $+      analyzeStatus status

@@ -2,9 +2,11 @@ module Main where  import qualified Test.Numeric.GLPK+import qualified Test.Numeric.GLPK.Monad  import qualified Test.DocTest.Driver as DocTest  main :: IO () main = DocTest.run $ do     Test.Numeric.GLPK.test+    Test.Numeric.GLPK.Monad.test

@@ -1,12 +1,13 @@ -- Do not edit! Automatically created with doctest-extract from src/Numeric/GLPK.hs-{-# LINE 68 "src/Numeric/GLPK.hs" #-}+{-# LINE 82 "src/Numeric/GLPK.hs" #-}  module Test.Numeric.GLPK where  import Test.DocTest.Base import qualified Test.DocTest.Driver as DocTest -{-# LINE 69 "src/Numeric/GLPK.hs" #-}+{-# LINE 83 "src/Numeric/GLPK.hs" #-}+import     qualified Test.Numeric.GLPK.Generator as TestLP import     qualified Numeric.GLPK as LP import     Numeric.GLPK ((.*), (<=.), (==.)) @@ -33,26 +34,71 @@  test :: DocTest.T () test = do- DocTest.printPrefix "Numeric.GLPK:198: "-{-# LINE 198 "src/Numeric/GLPK.hs" #-}+ DocTest.printPrefix "Numeric.GLPK:141: "+{-# LINE 141 "src/Numeric/GLPK.hs" #-}  DocTest.property-{-# LINE 198 "src/Numeric/GLPK.hs" #-}+{-# LINE 141 "src/Numeric/GLPK.hs" #-}+     (\target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (1,1) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (TestLP.approxReal 0.001 absol (abs target)) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False)+ DocTest.printPrefix "Numeric.GLPK:142: "+{-# LINE 142 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 142 "src/Numeric/GLPK.hs" #-}      (\(QC.Positive posWeight) (QC.Positive negWeight) target -> case Shape.indexTupleFromShape pairShape of (pos,neg) -> case mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[1.*pos, (-1).*neg] ==. target] (LP.Minimize, Array.fromTuple (posWeight,negWeight) :: Array.Array PairShape Double) of (Right (LP.Optimal,(absol,(posResult,negResult)))) -> QC.property (absol>=0) .&&. (posResult === 0 .||. negResult === 0); _ -> QC.property False)- DocTest.printPrefix "Numeric.GLPK:195: "-{-# LINE 195 "src/Numeric/GLPK.hs" #-}+ DocTest.printPrefix "Numeric.GLPK:143: "+{-# LINE 143 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 143 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case LP.simplex bounds constrs (dir,obj) of Right (LP.Optimal, _) -> True; _ -> False)+ DocTest.printPrefix "Numeric.GLPK:138: "+{-# LINE 138 "src/Numeric/GLPK.hs" #-}  DocTest.example-{-# LINE 195 "src/Numeric/GLPK.hs" #-}+{-# LINE 138 "src/Numeric/GLPK.hs" #-}    (case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.simplex [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]- DocTest.printPrefix "Numeric.GLPK:256: "-{-# LINE 256 "src/Numeric/GLPK.hs" #-}+ DocTest.printPrefix "Numeric.GLPK:180: "+{-# LINE 180 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 180 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case LP.simplexMulti bounds constrs (Array.shape origin) objs of (_, Right (LP.Optimal, _)) -> QC.property True; result -> QC.counterexample (show result) False)+ DocTest.printPrefix "Numeric.GLPK:219: "+{-# LINE 219 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 219 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.simplexSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False)+ DocTest.printPrefix "Numeric.GLPK:220: "+{-# LINE 220 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 220 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAllShrink (TestLP.genProblem origin) TestLP.shrinkProblem $ \(bounds, constrs) -> QC.forAllShrink (TestLP.genObjectives origin) TestLP.shrinkObjectives $ \objs -> case uncurry (LP.exactSuccessive bounds constrs) $ TestLP.successiveObjectives origin 0.01 objs of result -> QC.counterexample (show result) $ case result of Right results -> all (\r -> case r of (LP.Optimal, _) -> True; _ -> False) results; _ -> False)+ DocTest.printPrefix "Numeric.GLPK:240: "+{-# LINE 240 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 240 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.simplexSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.simplex bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB)+ DocTest.printPrefix "Numeric.GLPK:241: "+{-# LINE 241 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 241 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjectives origin) $ (. TestLP.successiveObjectives origin 0.01) $ \(obj,objs) -> case (LP.exactSuccessive bounds constrs obj objs, LP.solveSuccessive (LP.exact bounds) constrs obj objs) of (resultA,resultB) -> TestLP.approxSuccession 0.01 resultA resultB)+ DocTest.printPrefix "Numeric.GLPK:274: "+{-# LINE 274 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 274 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.exact bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False)+ DocTest.printPrefix "Numeric.GLPK:271: "+{-# LINE 271 "src/Numeric/GLPK.hs" #-}  DocTest.example-{-# LINE 256 "src/Numeric/GLPK.hs" #-}+{-# LINE 271 "src/Numeric/GLPK.hs" #-}    (case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapSnd Array.toTuple) <$> LP.exact [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]- DocTest.printPrefix "Numeric.GLPK:288: "-{-# LINE 288 "src/Numeric/GLPK.hs" #-}+ DocTest.printPrefix "Numeric.GLPK:346: "+{-# LINE 346 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 346 "src/Numeric/GLPK.hs" #-}+     (QC.forAllShrink TestLP.genOrigin TestLP.shrinkOrigin $ \origin -> QC.forAll (TestLP.genProblem origin) $ \(bounds, constrs) -> QC.forAll (TestLP.genObjective origin) $ \(dir,obj) -> case (LP.simplex bounds constrs (dir,obj), LP.interior bounds constrs (dir,obj)) of (Right (LP.Optimal, (optSimplex,_)), Right (LP.Optimal, (optExact,_))) -> TestLP.approx "optimum" 0.001 optSimplex optExact; _ -> QC.property False)+ DocTest.printPrefix "Numeric.GLPK:343: "+{-# LINE 343 "src/Numeric/GLPK.hs" #-}  DocTest.example-{-# LINE 288 "src/Numeric/GLPK.hs" #-}+{-# LINE 343 "src/Numeric/GLPK.hs" #-}    (case Shape.indexTupleFromShape tripletShape of (x1,x2,x3) -> mapSnd (mapPair (round3, Array.toTuple . Array.map round3)) <$> LP.interior [] [[2.*x1, 1.*x2] <=. 10, [1.*x2, 5.*x3] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double))   [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]

@@ -0,0 +1,208 @@+{-# LANGUAGE TypeFamilies #-}+module Test.Numeric.GLPK.Generator where++import qualified Numeric.GLPK as LP+import Numeric.GLPK ((<=.), (>=.))++import qualified Test.QuickCheck as QC+import System.Random (Random)++import qualified Data.Array.Comfort.Boxed as BoxedArray+import qualified Data.Array.Comfort.Storable as Array+import qualified Data.Array.Comfort.Shape as Shape+import qualified Data.NonEmpty as NonEmpty+import qualified Data.List.HT as ListHT+import qualified Data.Ix as Ix+import Data.Array.Comfort.Storable (Array, (!))+import Data.Traversable (sequenceA, for)+import Data.Tuple.HT (mapSnd)+import Data.Maybe (fromMaybe)+import Data.Int (Int64)++import Control.Applicative (liftA2)++import Text.Printf (PrintfArg, printf)++import Foreign.Storable (Storable)+++{- |+Generate constraints in the form of a polyhedron+which contains warrantedly the zero vector.+That is, there is an admissible solution.+In order to assert that the polyhedron is closed,+we bound all variables by a hypercube.+-}+genProblem ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>+   Array sh a -> QC.Gen (LP.Bounds ix, LP.Constraints ix)+genProblem origin =+   liftA2 (,)+      (for (Array.toAssociations origin) $ \(ix,x) ->+         LP.Inequality ix <$>+         liftA2 LP.Between+            (doubleFromElement . (x+) <$> QC.choose (-100,-50))+            (doubleFromElement . (x+) <$> QC.choose (50,100)))+      (do+         numConstraints <- QC.choose (1,20)+         QC.vectorOf numConstraints $ do+            ixs <- QC.sublistOf $ Shape.indices $ Array.shape origin+            terms <- for ixs $ \ix -> do+               coeff <- QC.choose (-10,10)+               return (coeff, ix)+            let offset = scalarProduct terms origin+            let deviation = 25+            LP.Inequality+               (map (uncurry (LP.Term . doubleFromElement)) terms)+               <$>+               QC.oneof (+                  (do bound <- QC.choose (offset-deviation, offset+deviation)+                      return $+                         if bound > offset+                            then LP.LessEqual    $ doubleFromElement bound+                            else LP.GreaterEqual $ doubleFromElement bound) :+                  (liftA2 LP.Between+                     (doubleFromElement <$>+                        QC.choose (offset-deviation, offset))+                     (doubleFromElement <$>+                        QC.choose (offset, offset+deviation))) :+                  []))++scalarProduct ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Storable a, Num a) =>+   [(a,ix)] -> Array sh a -> a+scalarProduct terms origin =+   sum $ map (\(coeff, ix) -> coeff * origin!ix) terms++genVarShape :: QC.Gen (Shape.Range Char)+genVarShape = Shape.Range 'a' <$> QC.choose ('a','j')++genOrigin :: QC.Gen (Array (Shape.Range Char) Int64)+genOrigin = genVector =<< genVarShape++_genOrigin :: QC.Gen (Array (Shape.Range Char) Double)+_genOrigin = genVector =<< genVarShape+++shrinkVarShape :: Shape.Range Char -> [Shape.Range Char]+shrinkVarShape (Shape.Range from to) =+   if from<to then [Shape.Range from (pred to)] else []++shrinkOrigin ::+   (Storable a) => Array (Shape.Range Char) a -> [Array (Shape.Range Char) a]+shrinkOrigin vec =+   case Array.shape vec of+      Shape.Range from to ->+         if from<to+            then [Array.sample (Shape.Range from (pred to)) (vec!)]+            else []+++class (Storable a, Random a, Num a, Ord a) => Element a where+   doubleFromElement :: a -> Double++instance Element Double where+   doubleFromElement = id++instance Element Int64 where+   doubleFromElement = fromIntegral++genObjective ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>+   Array sh a -> QC.Gen (LP.Direction, LP.Objective sh)+genObjective origin =+   liftA2 (,) QC.arbitraryBoundedEnum+      (fmap (Array.map doubleFromElement . flip asTypeOf origin) $+       genVector $ Array.shape origin)++genVector :: (Shape.Indexed sh, Element a) => sh -> QC.Gen (Array sh a)+genVector shape =+   fmap Array.fromBoxed $ sequenceA $+   BoxedArray.fromAssociations (QC.choose (-10,10)) shape []++shrinkProblem ::+   (LP.Bounds ix, LP.Constraints ix) ->+   [(LP.Bounds ix, LP.Constraints ix)]+shrinkProblem (bounds, constraints) =+   map (\shrinked -> (bounds, shrinked)) $+   filter (not . null) $ QC.shrinkList (const []) constraints++genObjectives ::+   (Shape.Indexed sh, Shape.Index sh ~ ix, Element a) =>+   Array sh a -> QC.Gen (NonEmpty.T [] (LP.Direction, [LP.Term ix]))+genObjectives origin = do+   let shape = Array.shape origin+   let stageRange :: (Int,Int)+       stageRange = (0,3)+   stages <- for (Shape.indices shape) $ \ix -> (,) ix <$> QC.choose stageRange+   let varSets =+         fromMaybe (error "there should be at least one stage") $+         NonEmpty.fetch $+         filter (not . null) $+         map (\k -> map fst $ filter ((k==) . snd) stages) $+         Ix.range stageRange+   let asTypeOfElement :: a -> f a -> a+       asTypeOfElement = const+   for varSets $ \varSet ->+      liftA2 (,)+         QC.arbitraryBoundedEnum+         (for varSet $ \ix ->+            flip LP.Term ix . doubleFromElement+               <$> QC.choose (-10, 10 `asTypeOfElement` origin))++shrinkObjectives ::+   NonEmpty.T [] (LP.Direction, [LP.Term ix]) ->+   [NonEmpty.T [] (LP.Direction, [LP.Term ix])]+shrinkObjectives (NonEmpty.Cons obj objs) =+   map (NonEmpty.Cons obj) $+   QC.shrinkList+      (\(dir,terms) ->+         map ((,) dir) $ filter (not . null) $+         QC.shrinkList (const []) terms)+      objs++successiveObjectives ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Array sh a -> Double ->+   NonEmpty.T [] (LP.Direction, [LP.Term ix]) ->+   ((LP.Direction, LP.Objective sh),+    [((LP.SolutionType, (Double, Array sh Double)) -> LP.Constraints ix,+      (LP.Direction, LP.Objective sh))])+successiveObjectives origin tol xs =+   let shape = Array.shape origin in+   (mapSnd (LP.objectiveFromTerms shape) $ NonEmpty.head xs,+    NonEmpty.mapAdjacent+      (\(dir0,obj0) y1 ->+         (\(_sol,(opt,_vec)) ->+            case dir0 of+               LP.Minimize -> [obj0 <=. opt + tol]+               LP.Maximize -> [obj0 >=. opt - tol],+          mapSnd (LP.objectiveFromTerms shape) y1))+      xs)+++approxReal :: (Ord a, Num a) => a -> a -> a -> Bool+approxReal tol x y = abs (x-y) <= tol++approx :: (PrintfArg a, Ord a, Num a) => String -> a -> a -> a -> QC.Property+approx name tol x y =+   QC.counterexample (printf "%s: %f - %f" name x y) (approxReal tol x y)++approxSuccession ::+   (Shape.C sh, Show sh, Show a, Ord a, Num a, Storable a) =>+   a ->+   Either LP.NoSolutionType+      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->+   Either LP.NoSolutionType+      (NonEmpty.T [] (LP.SolutionType, (a, Array sh a))) ->+   QC.Property+approxSuccession tol x y =+   QC.counterexample (show x) $+   QC.counterexample (show y) $+   case (x,y) of+      (Left sx, Left sy) -> sx==sy+      (Right (NonEmpty.Cons xh xs), Right (NonEmpty.Cons yh ys)) ->+         let equalSol (solX, (optX, _)) (solY, (optY, _)) =+               solX == solY && approxReal tol optX optY+         in equalSol xh yh  &&  ListHT.equalWith equalSol xs ys+      _ -> False

@@ -0,0 +1,26 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/GLPK/Monad.hs+{-# LINE 37 "src/Numeric/GLPK/Monad.hs" #-}++module Test.Numeric.GLPK.Monad where++import Test.DocTest.Base+import qualified Test.DocTest.Driver as DocTest++{-# LINE 38 "src/Numeric/GLPK/Monad.hs" #-}+import     qualified Numeric.GLPK.Monad as LP+import     Test.Numeric.GLPK (TripletShape, tripletShape)+import     Numeric.GLPK ((.*), (<=.))++import     qualified Data.Array.Comfort.Storable as Array+import     qualified Data.Array.Comfort.Shape as Shape++import     Data.Tuple.HT (mapSnd)++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.GLPK.Monad:67: "+{-# LINE 67 "src/Numeric/GLPK/Monad.hs" #-}+ DocTest.example+{-# LINE 67 "src/Numeric/GLPK/Monad.hs" #-}+   (case Shape.indexTupleFromShape tripletShape of (x,y,z) -> mapSnd (mapSnd Array.toTuple) <$> LP.run tripletShape [] (LP.simplex [[2.*x, 1.*y] <=. 10, [1.*y, (5::Double).*z] <=. 20] (LP.Maximize, Array.fromTuple (4,-3,2) :: Array.Array TripletShape Double)))+  [ExpectedLine [LineChunk "Right (Optimal,(28.0,(5.0,0.0,4.0)))"]]