{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
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)
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
simplex ::
(Shape.Indexed sh, Shape.Index sh ~ ix) =>
Bounds ix -> Constraints ix ->
(Direction, Objective sh) -> Solution sh
simplex = solve (flip FFI.glp_simplex nullPtr)
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 constrs0 sh (NonEmpty.Cons obj0 objs0) =
let go constrs curObj ((dir,obj):objs) (Right (Optimal, (opt,_))) =
mapFst (opt:) $
let extConstrs = (curObj==.opt) : constrs in
go extConstrs obj objs $
solver extConstrs
(dir, objectiveFromTerms sh obj)
go _ _ _ sol = ([], sol)
in go constrs0 (snd obj0) objs0 $
solver constrs0 $ mapSnd (objectiveFromTerms sh) obj0
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
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) :
[]