packages feed

comfort-glpk (empty) → 0.0

raw patch · 6 files changed

+541/−0 lines, 6 filesdep +QuickCheckdep +basedep +comfort-arraysetup-changed

Dependencies added: QuickCheck, base, comfort-array, comfort-glpk, deepseq, doctest-exitcode-stdio, doctest-lib, glpk-headers, non-empty, utility-ht

Files

+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) Henning Thielemann 2022++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+   notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+   notice, this list of conditions and the following disclaimer in the+   documentation and/or other materials provided with the distribution.+3. Neither the name of the author nor the names of his contributors+   may be used to endorse or promote products derived from this software+   without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ comfort-glpk.cabal view
@@ -0,0 +1,76 @@+Cabal-Version:    2.2+Name:             comfort-glpk+Version:          0.0+License:          BSD-3-Clause+License-File:     LICENSE+Author:           Henning Thielemann <haskell@henning-thielemann.de>+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>+Category:         Math+Tested-With:      GHC ==8.6.5+Build-Type:       Simple+Synopsis:         Linear Programming using GLPK and comfort-array+Description:+  Simple interface to linear programming functions provided by GLPK+  using the flexible Array shape framework from @comfort-array@.+  .+  E.g. you can use @Shape.Tuple@ to convert safely+  between nested tuples and arrays with the same number of elements.+  .+  > type X = Shape.Element+  > type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)+  >+  > case Shape.indexTupleFromShape (Shape.static :: PairShape) of+  >   (posIx,negIx) ->+  >     case mapSnd (mapSnd Array.toTuple) <$>+  >          LP.simplex [] [[1.*posIx, (-1).*negIx] ==. 314]+  >            (LP.Minimize,+  >             Array.fromTuple (23,42) :: Array PairShape Double)+  >       of+  >         (Right (LP.Optimal, (absol, (pos, neg)))) ->+  >           printf "absol %f,  pos %f, neg %f\n" absol pos neg+  >         _ -> fail "GLPK solver failed"+  .+  Alternatives: @hmatrix-glpk@, @glpk-hs@++Source-Repository this+  Tag:         0.0+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/comfort-glpk/++Source-Repository head+  Type:        darcs+  Location:    https://hub.darcs.net/thielema/comfort-glpk/++Library+  Build-Depends:+    glpk-headers >=0.4.1 && <0.5,+    comfort-array >=0.4 && <0.6,+    deepseq >=1.3 && <1.5,+    non-empty >=0.3.2 && <0.4,+    utility-ht >=0.0.16 && <0.1,+    base >=4.5 && <5++  GHC-Options:      -Wall+  Hs-Source-Dirs:   src+  Default-Language: Haskell98+  Extra-Libraries: glpk+  Exposed-modules:+    Numeric.GLPK++Test-Suite comfort-glpk-test+  Type:             exitcode-stdio-1.0+  Build-Depends:+    comfort-glpk,+    comfort-array >=0.5.2,+    utility-ht,+    doctest-exitcode-stdio >=0.0 && <0.1,+    doctest-lib >=0.1 && <0.2,+    QuickCheck >=2.1 && <3,+    base >=4.5 && <5++  GHC-Options:      -Wall+  Hs-Source-Dirs:   test+  Default-Language: Haskell98+  Main-Is: Main.hs+  Other-Modules:+    Test.Numeric.GLPK
+ src/Numeric/GLPK.hs view
@@ -0,0 +1,367 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{- |+The following LP problem++maximize @4 x_1 - 3 x_2 + 2 x_3@ subject to++@2 x_1 + x_2 <= 10@++@x_2 + 5 x_3 <= 20@++and++@x_i >= 0@++is used as an example in the doctest comments.+++By default all indeterminates are non-negative.+A given bound for a variable completely replaces the default,+so @0 <= x_i <= b@ must be explicitly given as @i >=<. (0,b)@.+Multiple bounds for a variable are not allowed,+instead of @[i >=. a, i <=. b]@ use @i >=<. (a,b)@.+-}+module Numeric.GLPK (+   Term(..),+   Bound(..),+   Inequality(..),+   free, (<=.), (>=.), (==.), (>=<.),+   NoSolutionType(..),+   SolutionType(..),+   Solution,+   Constraints,+   Direction(..),+   Objective,+   Bounds,+   (.*),+   objectiveFromTerms,+   simplex,+   simplexMulti,+   exact,+   interior,+   interiorMulti,+   ) 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 qualified Data.NonEmpty as NonEmpty+import Data.Array.Comfort.Storable (Array)+import Data.Tuple.HT (mapFst, mapSnd)+import Data.Maybe (fromMaybe)+import Data.Foldable (traverse_, for_)++import Control.Monad (void)+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 Numeric.GLPK as LP+>>> import Numeric.GLPK ((.*), (<=.), (==.))+>>>+>>> import qualified Data.Array.Comfort.Storable as Array+>>> import qualified Data.Array.Comfort.Shape as Shape+>>>+>>> import Data.Tuple.HT (mapPair, mapSnd)+>>>+>>> import qualified Test.QuickCheck as QC+>>> import Test.QuickCheck ((===), (.&&.), (.||.))+>>>+>>> type X = Shape.Element+>>> type PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)+>>> type TripletShape = Shape.NestedTuple Shape.TupleIndex (X,X,X)+>>>+>>> pairShape :: PairShape+>>> pairShape = Shape.static+>>>+>>> tripletShape :: TripletShape+>>> tripletShape = Shape.static+>>>+>>> round3 :: Double -> Double+>>> round3 x = fromInteger (round (1000*x)) / 1000+-}+++infix 7 .*++(.*) :: Double -> ix -> Term ix+(.*) = Term++data Term ix = Term Double ix+   deriving (Show)+++infix 4 <=., >=., >=<., ==.++(<=.), (>=.), (==.) :: x -> Double -> Inequality x+x <=. bnd = Inequality x $ LessEqual bnd+x >=. bnd = Inequality x $ GreaterEqual bnd+x ==. bnd = Inequality x $ Equal bnd++(>=<.) :: x -> (Double,Double) -> Inequality x+x >=<. bnd = Inequality x $ uncurry Between bnd++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> \(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+-}+simplex ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) -> Solution sh+simplex = solve (flip FFI.glp_simplex nullPtr)++{- |+Optimize for one objective after another.+That is, if the first optimization succeeds+then the optimum is fixed as constraint+and the optimization is continued with respect to the second objective and so on.+The iteration fails if one optimization fails.+The obtained objective values are returned as well.+Their number equals the number of attempted optimizations.++The last objective value is included in the Solution value.+This is a bit inconsistent,+but this way you have a warranty that there is an objective value+if the optimization is successful.++The objectives are expected as 'Term's+because after successful optimization step+they are used as (sparse) constraints.+It's also easy to assert that the same array shape+is used for all objectives.++The function does not work reliably,+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.+-}+simplexMulti, interiorMulti ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Bounds ix -> Constraints ix ->+   sh -> NonEmpty.T [] (Direction, [Term ix]) -> ([Double], Solution sh)+simplexMulti = solveMulti . simplex+interiorMulti = solveMulti . interior++solveMulti ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   (Constraints ix -> (Direction, Objective sh) -> Solution sh) ->+   Constraints ix ->+   sh -> NonEmpty.T [] (Direction, [Term ix]) -> ([Double], Solution sh)+solveMulti solver constrs sh (NonEmpty.Cons obj0 objs0) =+   let go curObj ((dir,obj):objs) (Right (Optimal, (opt,_))) =+         mapFst (opt:) $+         go obj objs $+         solver+            ((curObj==.opt) : constrs)+            (dir, objectiveFromTerms sh obj)+       go _ _ sol = ([], sol)+   in go (snd obj0) objs0 $+      solver constrs $ mapSnd (objectiveFromTerms sh) obj0+++{- |+>>> 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)))+-}+exact ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) -> Solution sh+exact = solve (flip FFI.glp_exact nullPtr)+++{-# INLINE solve #-}+solve ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   (Foreign.Ptr FFI.Problem -> IO FFI.GlpkSimplexStatus) ->+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) -> Solution sh+solve solver bounds constrs (dir,obj) = unsafePerformIO $+   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+++{- |+>>> 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)))+-}+interior ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   Bounds ix -> Constraints ix ->+   (Direction, Objective sh) -> Solution sh+interior bounds constrs (dir,obj) = unsafePerformIO $+   bracket FFI.glp_create_prob FFI.glp_delete_prob $ \lp -> do+   storeProblem bounds constrs (dir,obj) lp+   void $ FFI.glp_interior lp nullPtr+   let examine =+         liftA2 (,)+            (realToFrac <$> FFI.glp_ipt_obj_val lp)+            (readGLPArray (Array.shape obj) $ \arr ix ->+               Mutable.write arr ix . realToFrac+                  =<< FFI.glp_ipt_col_prim lp (deferredColumnIndex ix))+   status <- FFI.glp_ipt_status lp+   either (return . Left) (\typ -> Right . (,) typ <$> examine) $+      analyzeStatus status+++storeProblem ::+   (Shape.Indexed sh, Shape.Index sh ~ ix) =>+   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+   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+      for_ (zip [1..] $ concat $+            zipWith (map . (,)) [firstRow..] $+            map (fst . prepareBounds) constrs) $+         \(k, (row, Term c x)) -> do+            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++{-# 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) :+      []
+ test/Main.hs view
@@ -0,0 +1,10 @@+-- Do not edit! Automatically created with doctest-extract.+module Main where++import qualified Test.Numeric.GLPK++import qualified Test.DocTest.Driver as DocTest++main :: IO ()+main = DocTest.run $ do+    Test.Numeric.GLPK.test
+ test/Test/Numeric/GLPK.hs view
@@ -0,0 +1,58 @@+-- Do not edit! Automatically created with doctest-extract from src/Numeric/GLPK.hs+{-# LINE 68 "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" #-}+import     qualified Numeric.GLPK as LP+import     Numeric.GLPK ((.*), (<=.), (==.))++import     qualified Data.Array.Comfort.Storable as Array+import     qualified Data.Array.Comfort.Shape as Shape++import     Data.Tuple.HT (mapPair, mapSnd)++import     qualified Test.QuickCheck as QC+import     Test.QuickCheck ((===), (.&&.), (.||.))++type     X = Shape.Element+type     PairShape = Shape.NestedTuple Shape.TupleIndex (X,X)+type     TripletShape = Shape.NestedTuple Shape.TupleIndex (X,X,X)++pairShape     :: PairShape+pairShape     = Shape.static++tripletShape     :: TripletShape+tripletShape     = Shape.static++round3     :: Double -> Double+round3     x = fromInteger (round (1000*x)) / 1000++test :: DocTest.T ()+test = do+ DocTest.printPrefix "Numeric.GLPK:198: "+{-# LINE 198 "src/Numeric/GLPK.hs" #-}+ DocTest.property+{-# LINE 198 "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.example+{-# LINE 195 "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.example+{-# LINE 256 "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.example+{-# LINE 288 "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)))"]]