diff --git a/comfort-glpk.cabal b/comfort-glpk.cabal --- a/comfort-glpk.cabal +++ b/comfort-glpk.cabal @@ -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 @@ -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 diff --git a/src/Numeric/GLPK.hs b/src/Numeric/GLPK.hs --- a/src/Numeric/GLPK.hs +++ b/src/Numeric/GLPK.hs @@ -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;" : + [] diff --git a/src/Numeric/GLPK/Monad.hs b/src/Numeric/GLPK/Monad.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/GLPK/Monad.hs @@ -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 diff --git a/src/Numeric/GLPK/Private.hs b/src/Numeric/GLPK/Private.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/GLPK/Private.hs @@ -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 diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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 diff --git a/test/Test/Numeric/GLPK.hs b/test/Test/Numeric/GLPK.hs --- a/test/Test/Numeric/GLPK.hs +++ b/test/Test/Numeric/GLPK.hs @@ -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)))"]] diff --git a/test/Test/Numeric/GLPK/Generator.hs b/test/Test/Numeric/GLPK/Generator.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/GLPK/Generator.hs @@ -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 Array (Shape.Range Char) a -> [Array (Shape.Range Char) a] +shrinkOrigin vec = + case Array.shape vec of + Shape.Range from to -> + if from 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 diff --git a/test/Test/Numeric/GLPK/Monad.hs b/test/Test/Numeric/GLPK/Monad.hs new file mode 100644 --- /dev/null +++ b/test/Test/Numeric/GLPK/Monad.hs @@ -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)))"]]