{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
module Numeric.GLPK (
Term(..),
Bound(..),
Inequality(..),
free, (<=.), (>=.), (==.), (>=<.),
NoSolutionType(..),
SolutionType(..),
Solution,
Constraints,
Direction(..),
Objective,
Bounds,
(.*),
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.Traversable (for)
import Data.Foldable (for_)
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 System.IO.Unsafe (unsafePerformIO)
import qualified Foreign
import Foreign.Ptr (nullPtr)
infix 7 .*
(.*) :: Double -> ix -> Term ix
(.*) = Term
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
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))
simplex ::
(Shape.Indexed sh, Shape.Index sh ~ ix) =>
Bounds ix -> Constraints ix ->
(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" #-}
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 ::
(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
{-# DEPRECATED simplexSuccessive "use GLPK.Monad instead" #-}
{-# DEPRECATED exactSuccessive "use GLPK.Monad instead" #-}
{-# DEPRECATED interiorSuccessive "run 'interior' in Either monad instead" #-}
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" #-}
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
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
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
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
void $ FFI.glp_term_out FFI.glpkOff
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
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) $
\(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 numTerms) ia ja ar
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;" :
[]