{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.Simplex.Textbook.MIPSolver.Simple
-- Copyright   :  (c) Masahiro Sakai 2011
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- References:
--
-- * [Gomory1960]
--   Ralph E. Gomory.
--   An Algorithm for the Mixed Integer Problem, Technical Report
--   RM-2597, 1960, The Rand Corporation, Santa Monica, CA.
--   <http://www.rand.org/pubs/research_memoranda/RM2597.html>
--
-- * [Gomory1958]
--   Ralph E. Gomory.
--   Outline of an algorithm for integer solutions to linear programs.
--   Bull. Amer. Math. Soc., Vol. 64, No. 5. (1958), pp. 275-278.
--   <http://projecteuclid.org/euclid.bams/1183522679>
-----------------------------------------------------------------------------

module ToySolver.Arith.Simplex.Textbook.MIPSolver.Simple
  ( module Data.OptDir
  , OptResult (..)
  , minimize
  , maximize
  , optimize
  ) where

import Control.Exception
import Control.Monad
import Control.Monad.State
import Data.Default.Class
import Data.Ord
import Data.Maybe
import Data.List (maximumBy)
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.OptDir
import Data.VectorSpace

import ToySolver.Data.OrdRel
import ToySolver.Data.IntVar
import qualified ToySolver.Data.LA as LA
import qualified ToySolver.Arith.Simplex.Textbook as Simplex
import qualified ToySolver.Arith.Simplex.Textbook.LPSolver as LPSolver
import ToySolver.Arith.Simplex.Textbook.LPSolver hiding (OptResult (..))
import ToySolver.Arith.Simplex.Textbook.LPSolver.Simple (OptResult (..))
import qualified ToySolver.Arith.OmegaTest as OmegaTest
import ToySolver.Internal.Util (isInteger, fracPart)

-- ---------------------------------------------------------------------------

data Node r
  = Node
  { forall r. Node r -> Solver r
ndSolver :: LPSolver.Solver r
  , forall r. Node r -> Int
ndDepth  :: {-# UNPACK #-} !Int
--  , ndCutSlackVariables :: VarSet
  }

ndTableau :: Node r  -> Simplex.Tableau r
ndTableau :: forall r. Node r -> Tableau r
ndTableau Node r
node = State (Solver r) (Tableau r) -> Solver r -> Tableau r
forall s a. State s a -> s -> a
evalState State (Solver r) (Tableau r)
forall r. LP r (Tableau r)
getTableau (Node r -> Solver r
forall r. Node r -> Solver r
ndSolver Node r
node)

ndLowerBound :: Node r -> r
ndLowerBound :: forall r. Node r -> r
ndLowerBound Node r
node = State (Solver r) r -> Solver r -> r
forall s a. State s a -> s -> a
evalState ((Tableau r -> r)
-> StateT (Solver r) Identity (Tableau r) -> State (Solver r) r
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM Tableau r -> r
forall r. Tableau r -> r
Simplex.currentObjValue StateT (Solver r) Identity (Tableau r)
forall r. LP r (Tableau r)
getTableau) (Node r -> Solver r
forall r. Node r -> Solver r
ndSolver Node r
node)

data Err = ErrUnbounded | ErrUnsat deriving (Eq Err
Eq Err =>
(Err -> Err -> Ordering)
-> (Err -> Err -> Bool)
-> (Err -> Err -> Bool)
-> (Err -> Err -> Bool)
-> (Err -> Err -> Bool)
-> (Err -> Err -> Err)
-> (Err -> Err -> Err)
-> Ord Err
Err -> Err -> Bool
Err -> Err -> Ordering
Err -> Err -> Err
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Err -> Err -> Ordering
compare :: Err -> Err -> Ordering
$c< :: Err -> Err -> Bool
< :: Err -> Err -> Bool
$c<= :: Err -> Err -> Bool
<= :: Err -> Err -> Bool
$c> :: Err -> Err -> Bool
> :: Err -> Err -> Bool
$c>= :: Err -> Err -> Bool
>= :: Err -> Err -> Bool
$cmax :: Err -> Err -> Err
max :: Err -> Err -> Err
$cmin :: Err -> Err -> Err
min :: Err -> Err -> Err
Ord, Err -> Err -> Bool
(Err -> Err -> Bool) -> (Err -> Err -> Bool) -> Eq Err
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Err -> Err -> Bool
== :: Err -> Err -> Bool
$c/= :: Err -> Err -> Bool
/= :: Err -> Err -> Bool
Eq, Int -> Err -> ShowS
[Err] -> ShowS
Err -> String
(Int -> Err -> ShowS)
-> (Err -> String) -> ([Err] -> ShowS) -> Show Err
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Err -> ShowS
showsPrec :: Int -> Err -> ShowS
$cshow :: Err -> String
show :: Err -> String
$cshowList :: [Err] -> ShowS
showList :: [Err] -> ShowS
Show, Int -> Err
Err -> Int
Err -> [Err]
Err -> Err
Err -> Err -> [Err]
Err -> Err -> Err -> [Err]
(Err -> Err)
-> (Err -> Err)
-> (Int -> Err)
-> (Err -> Int)
-> (Err -> [Err])
-> (Err -> Err -> [Err])
-> (Err -> Err -> [Err])
-> (Err -> Err -> Err -> [Err])
-> Enum Err
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
$csucc :: Err -> Err
succ :: Err -> Err
$cpred :: Err -> Err
pred :: Err -> Err
$ctoEnum :: Int -> Err
toEnum :: Int -> Err
$cfromEnum :: Err -> Int
fromEnum :: Err -> Int
$cenumFrom :: Err -> [Err]
enumFrom :: Err -> [Err]
$cenumFromThen :: Err -> Err -> [Err]
enumFromThen :: Err -> Err -> [Err]
$cenumFromTo :: Err -> Err -> [Err]
enumFromTo :: Err -> Err -> [Err]
$cenumFromThenTo :: Err -> Err -> Err -> [Err]
enumFromThenTo :: Err -> Err -> Err -> [Err]
Enum, Err
Err -> Err -> Bounded Err
forall a. a -> a -> Bounded a
$cminBound :: Err
minBound :: Err
$cmaxBound :: Err
maxBound :: Err
Bounded)

maximize :: RealFrac r => LA.Expr r -> [LA.Atom r] -> VarSet -> OptResult r
maximize :: forall r. RealFrac r => Expr r -> [Atom r] -> VarSet -> OptResult r
maximize = OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
forall r.
RealFrac r =>
OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
optimize OptDir
OptMax

minimize :: RealFrac r => LA.Expr r -> [LA.Atom r] -> VarSet -> OptResult r
minimize :: forall r. RealFrac r => Expr r -> [Atom r] -> VarSet -> OptResult r
minimize = OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
forall r.
RealFrac r =>
OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
optimize OptDir
OptMin

optimize :: RealFrac r => OptDir -> LA.Expr r -> [LA.Atom r] -> VarSet -> OptResult r
optimize :: forall r.
RealFrac r =>
OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
optimize OptDir
optdir Expr r
obj [Atom r]
cs VarSet
ivs =
  case OptDir
-> Expr r -> [Atom r] -> VarSet -> Either Err (Node r, VarSet)
forall r.
RealFrac r =>
OptDir
-> Expr r -> [Atom r] -> VarSet -> Either Err (Node r, VarSet)
mkInitialNode OptDir
optdir Expr r
obj [Atom r]
cs VarSet
ivs of
    Left Err
err ->
      case Err
err of
        Err
ErrUnsat -> OptResult r
forall r. OptResult r
OptUnsat
        Err
ErrUnbounded ->
          if VarSet -> Bool
IS.null VarSet
ivs
          then OptResult r
forall r. OptResult r
Unbounded
          else
            {-
               Fallback to Fourier-Motzkin + OmegaTest
               * In general, original problem may have optimal
                 solution even though LP relaxiation is unbounded.
               * But if restricted to rational numbers, the
                 original problem is unbounded or unsatisfiable
                 when LP relaxation is unbounded.
            -}
            case Options
-> VarSet -> [Atom Rational] -> VarSet -> Maybe (Model Rational)
OmegaTest.solveQFLIRAConj Options
forall a. Default a => a
def ([Atom r] -> VarSet
forall a. Variables a => a -> VarSet
vars [Atom r]
cs VarSet -> VarSet -> VarSet
`IS.union` VarSet
ivs) ((Atom r -> Atom Rational) -> [Atom r] -> [Atom Rational]
forall a b. (a -> b) -> [a] -> [b]
map Atom r -> Atom Rational
forall r. RealFrac r => Atom r -> Atom Rational
conv [Atom r]
cs) VarSet
ivs of
              Maybe (Model Rational)
Nothing -> OptResult r
forall r. OptResult r
OptUnsat
              Just Model Rational
_ -> OptResult r
forall r. OptResult r
Unbounded
    Right (Node r
node0, VarSet
ivs2) ->
      case OptDir -> Expr r -> VarSet -> Node r -> Either Err (Node r)
forall r.
RealFrac r =>
OptDir -> Expr r -> VarSet -> Node r -> Either Err (Node r)
traverseBBTree OptDir
optdir Expr r
obj VarSet
ivs2 Node r
node0 of
        Left Err
ErrUnbounded -> String -> OptResult r
forall a. HasCallStack => String -> a
error String
"shoud not happen"
        Left Err
ErrUnsat -> OptResult r
forall r. OptResult r
OptUnsat
        Right Node r
node -> (State (Solver r) (OptResult r) -> Solver r -> OptResult r)
-> Solver r -> State (Solver r) (OptResult r) -> OptResult r
forall a b c. (a -> b -> c) -> b -> a -> c
flip State (Solver r) (OptResult r) -> Solver r -> OptResult r
forall s a. State s a -> s -> a
evalState (Node r -> Solver r
forall r. Node r -> Solver r
ndSolver Node r
node) (State (Solver r) (OptResult r) -> OptResult r)
-> State (Solver r) (OptResult r) -> OptResult r
forall a b. (a -> b) -> a -> b
$ do
          Tableau r
tbl <- LP r (Tableau r)
forall r. LP r (Tableau r)
getTableau
          Model r
model <- VarSet -> LP r (Model r)
forall r. Fractional r => VarSet -> LP r (Model r)
getModel VarSet
vs
          OptResult r -> State (Solver r) (OptResult r)
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (OptResult r -> State (Solver r) (OptResult r))
-> OptResult r -> State (Solver r) (OptResult r)
forall a b. (a -> b) -> a -> b
$ r -> Model r -> OptResult r
forall r. r -> Model r -> OptResult r
Optimum (Tableau r -> r
forall r. Tableau r -> r
Simplex.currentObjValue Tableau r
tbl) Model r
model
  where
    vs :: VarSet
vs = [Atom r] -> VarSet
forall a. Variables a => a -> VarSet
vars [Atom r]
cs VarSet -> VarSet -> VarSet
`IS.union` Expr r -> VarSet
forall a. Variables a => a -> VarSet
vars Expr r
obj

tableau' :: (RealFrac r) => [LA.Atom r] -> VarSet -> LP r VarSet
tableau' :: forall r. RealFrac r => [Atom r] -> VarSet -> LP r VarSet
tableau' [Atom r]
cs VarSet
ivs = do
  let (VarSet
nonnegVars, [Atom r]
cs') = [Atom r] -> VarSet -> (VarSet, [Atom r])
forall r. RealFrac r => [Atom r] -> VarSet -> (VarSet, [Atom r])
collectNonnegVars [Atom r]
cs VarSet
ivs
      fvs :: VarSet
fvs = [Atom r] -> VarSet
forall a. Variables a => a -> VarSet
vars [Atom r]
cs VarSet -> VarSet -> VarSet
`IS.difference` VarSet
nonnegVars
  VarSet
ivs2 <- ([VarSet] -> VarSet)
-> StateT (Solver r) Identity [VarSet] -> LP r VarSet
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM [VarSet] -> VarSet
forall (f :: * -> *). Foldable f => f VarSet -> VarSet
IS.unions (StateT (Solver r) Identity [VarSet] -> LP r VarSet)
-> StateT (Solver r) Identity [VarSet] -> LP r VarSet
forall a b. (a -> b) -> a -> b
$ [Int]
-> (Int -> LP r VarSet) -> StateT (Solver r) Identity [VarSet]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM (VarSet -> [Int]
IS.toList VarSet
fvs) ((Int -> LP r VarSet) -> StateT (Solver r) Identity [VarSet])
-> (Int -> LP r VarSet) -> StateT (Solver r) Identity [VarSet]
forall a b. (a -> b) -> a -> b
$ \Int
v -> do
    Int
v1 <- LP r Int
forall r. LP r Int
newVar
    Int
v2 <- LP r Int
forall r. LP r Int
newVar
    Int -> Expr r -> LP r ()
forall r. Int -> Expr r -> LP r ()
define Int
v (Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
v1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
v2)
    VarSet -> LP r VarSet
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (VarSet -> LP r VarSet) -> VarSet -> LP r VarSet
forall a b. (a -> b) -> a -> b
$ if Int
v Int -> VarSet -> Bool
`IS.member` VarSet
ivs then [Int] -> VarSet
IS.fromList [Int
v1,Int
v2] else VarSet
IS.empty
  (Atom r -> LP r ()) -> [Atom r] -> LP r ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ Atom r -> LP r ()
forall r. Real r => Atom r -> LP r ()
addConstraint [Atom r]
cs'
  VarSet -> LP r VarSet
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return VarSet
ivs2

conv :: RealFrac r => LA.Atom r -> LA.Atom Rational
conv :: forall r. RealFrac r => Atom r -> Atom Rational
conv = (Expr r -> Expr Rational) -> OrdRel (Expr r) -> Atom Rational
forall a b. (a -> b) -> OrdRel a -> OrdRel b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((r -> Rational) -> Expr r -> Expr Rational
forall b a. (Num b, Eq b) => (a -> b) -> Expr a -> Expr b
LA.mapCoeff r -> Rational
forall a. Real a => a -> Rational
toRational)

mkInitialNode :: RealFrac r => OptDir -> LA.Expr r -> [LA.Atom r] -> VarSet -> Either Err (Node r, VarSet)
mkInitialNode :: forall r.
RealFrac r =>
OptDir
-> Expr r -> [Atom r] -> VarSet -> Either Err (Node r, VarSet)
mkInitialNode OptDir
optdir Expr r
obj [Atom r]
cs VarSet
ivs =
  (State (Solver r) (Either Err (Node r, VarSet))
 -> Solver r -> Either Err (Node r, VarSet))
-> Solver r
-> State (Solver r) (Either Err (Node r, VarSet))
-> Either Err (Node r, VarSet)
forall a b c. (a -> b -> c) -> b -> a -> c
flip State (Solver r) (Either Err (Node r, VarSet))
-> Solver r -> Either Err (Node r, VarSet)
forall s a. State s a -> s -> a
evalState (VarSet -> Solver r
forall r. VarSet -> Solver r
emptySolver VarSet
vs) (State (Solver r) (Either Err (Node r, VarSet))
 -> Either Err (Node r, VarSet))
-> State (Solver r) (Either Err (Node r, VarSet))
-> Either Err (Node r, VarSet)
forall a b. (a -> b) -> a -> b
$ do
    VarSet
ivs2 <- [Atom r] -> VarSet -> LP r VarSet
forall r. RealFrac r => [Atom r] -> VarSet -> LP r VarSet
tableau' [Atom r]
cs VarSet
ivs
    OptResult
ret <- OptDir -> Expr r -> LP r OptResult
forall r.
(Fractional r, Real r) =>
OptDir -> Expr r -> LP r OptResult
LPSolver.twoPhaseSimplex OptDir
optdir Expr r
obj
    case OptResult
ret of
      OptResult
LPSolver.Unsat -> Either Err (Node r, VarSet)
-> State (Solver r) (Either Err (Node r, VarSet))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Err -> Either Err (Node r, VarSet)
forall a b. a -> Either a b
Left Err
ErrUnsat)
      OptResult
LPSolver.Unbounded -> Either Err (Node r, VarSet)
-> State (Solver r) (Either Err (Node r, VarSet))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Err -> Either Err (Node r, VarSet)
forall a b. a -> Either a b
Left Err
ErrUnbounded)
      OptResult
LPSolver.Optimum -> do
        Solver r
solver <- StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
        Either Err (Node r, VarSet)
-> State (Solver r) (Either Err (Node r, VarSet))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return (Either Err (Node r, VarSet)
 -> State (Solver r) (Either Err (Node r, VarSet)))
-> Either Err (Node r, VarSet)
-> State (Solver r) (Either Err (Node r, VarSet))
forall a b. (a -> b) -> a -> b
$ (Node r, VarSet) -> Either Err (Node r, VarSet)
forall a b. b -> Either a b
Right ((Node r, VarSet) -> Either Err (Node r, VarSet))
-> (Node r, VarSet) -> Either Err (Node r, VarSet)
forall a b. (a -> b) -> a -> b
$
          ( Node{ ndSolver :: Solver r
ndSolver = Solver r
solver
                , ndDepth :: Int
ndDepth = Int
0
--                , ndCutSlackVariables = IS.empty
                }
          , VarSet
ivs VarSet -> VarSet -> VarSet
`IS.union` VarSet
ivs2
          )
  where
    vs :: VarSet
vs = [Atom r] -> VarSet
forall a. Variables a => a -> VarSet
vars [Atom r]
cs VarSet -> VarSet -> VarSet
`IS.union` Expr r -> VarSet
forall a. Variables a => a -> VarSet
vars Expr r
obj

isStrictlyBetter :: RealFrac r => OptDir -> r -> r -> Bool
isStrictlyBetter :: forall r. RealFrac r => OptDir -> r -> r -> Bool
isStrictlyBetter OptDir
optdir = if OptDir
optdirOptDir -> OptDir -> Bool
forall a. Eq a => a -> a -> Bool
==OptDir
OptMin then r -> r -> Bool
forall a. Ord a => a -> a -> Bool
(<) else r -> r -> Bool
forall a. Ord a => a -> a -> Bool
(>)

traverseBBTree :: forall r. RealFrac r => OptDir -> LA.Expr r -> VarSet -> Node r -> Either Err (Node r)
traverseBBTree :: forall r.
RealFrac r =>
OptDir -> Expr r -> VarSet -> Node r -> Either Err (Node r)
traverseBBTree OptDir
optdir Expr r
obj VarSet
ivs Node r
node0 = [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop [Node r
node0] Maybe (Node r)
forall a. Maybe a
Nothing
  where
    loop :: [Node r] -> Maybe (Node r) -> Either Err (Node r)
    loop :: [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop [] (Just Node r
best) = Node r -> Either Err (Node r)
forall a b. b -> Either a b
Right Node r
best
    loop [] Maybe (Node r)
Nothing = Err -> Either Err (Node r)
forall a b. a -> Either a b
Left Err
ErrUnsat
    loop (Node r
n:[Node r]
ns) Maybe (Node r)
Nothing =
      case Node r -> Maybe [Node r]
children Node r
n of
        Maybe [Node r]
Nothing -> [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop [Node r]
ns (Node r -> Maybe (Node r)
forall a. a -> Maybe a
Just Node r
n)
        Just [Node r]
cs -> [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop ([Node r]
cs[Node r] -> [Node r] -> [Node r]
forall a. [a] -> [a] -> [a]
++[Node r]
ns) Maybe (Node r)
forall a. Maybe a
Nothing
    loop (Node r
n:[Node r]
ns) (Just Node r
best)
      | OptDir -> r -> r -> Bool
forall r. RealFrac r => OptDir -> r -> r -> Bool
isStrictlyBetter OptDir
optdir (Node r -> r
forall r. Node r -> r
ndLowerBound Node r
n) (Node r -> r
forall r. Node r -> r
ndLowerBound Node r
best) =
          case Node r -> Maybe [Node r]
children Node r
n of
            Maybe [Node r]
Nothing -> [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop [Node r]
ns (Node r -> Maybe (Node r)
forall a. a -> Maybe a
Just Node r
n)
            Just [Node r]
cs -> [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop ([Node r]
cs[Node r] -> [Node r] -> [Node r]
forall a. [a] -> [a] -> [a]
++[Node r]
ns) (Node r -> Maybe (Node r)
forall a. a -> Maybe a
Just Node r
best)
      | Bool
otherwise = [Node r] -> Maybe (Node r) -> Either Err (Node r)
loop [Node r]
ns (Node r -> Maybe (Node r)
forall a. a -> Maybe a
Just Node r
best)

    reopt :: Solver r -> Maybe (Solver r)
    reopt :: Solver r -> Maybe (Solver r)
reopt Solver r
s = (State (Solver r) (Maybe (Solver r))
 -> Solver r -> Maybe (Solver r))
-> Solver r
-> State (Solver r) (Maybe (Solver r))
-> Maybe (Solver r)
forall a b c. (a -> b -> c) -> b -> a -> c
flip State (Solver r) (Maybe (Solver r)) -> Solver r -> Maybe (Solver r)
forall s a. State s a -> s -> a
evalState Solver r
s (State (Solver r) (Maybe (Solver r)) -> Maybe (Solver r))
-> State (Solver r) (Maybe (Solver r)) -> Maybe (Solver r)
forall a b. (a -> b) -> a -> b
$ do
      Bool
ret <- OptDir -> Expr r -> LP r Bool
forall r. (Fractional r, Real r) => OptDir -> Expr r -> LP r Bool
dualSimplex OptDir
optdir Expr r
obj
      if Bool
ret
        then (Solver r -> Maybe (Solver r))
-> StateT (Solver r) Identity (Solver r)
-> State (Solver r) (Maybe (Solver r))
forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM Solver r -> Maybe (Solver r)
forall a. a -> Maybe a
Just StateT (Solver r) Identity (Solver r)
forall s (m :: * -> *). MonadState s m => m s
get
        else Maybe (Solver r) -> State (Solver r) (Maybe (Solver r))
forall a. a -> StateT (Solver r) Identity a
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Solver r)
forall a. Maybe a
Nothing

    children :: Node r -> Maybe [Node r]
    children :: Node r -> Maybe [Node r]
children Node r
node
      | [(Int, VarMap r, r)] -> Bool
forall a. [a] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Int, VarMap r, r)]
xs = Maybe [Node r]
forall a. Maybe a
Nothing -- no violation
      | Node r -> Int
forall r. Node r -> Int
ndDepth Node r
node Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
100 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 = -- cut
          let
            (r
f0, VarMap r
m0) = ((r, VarMap r) -> (r, VarMap r) -> Ordering)
-> [(r, VarMap r)] -> (r, VarMap r)
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (((r, VarMap r) -> r) -> (r, VarMap r) -> (r, VarMap r) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (r, VarMap r) -> r
forall a b. (a, b) -> a
fst) [(r -> r
forall a. RealFrac a => a -> a
fracPart r
val, VarMap r
m) | (Int
_,VarMap r
m,r
val) <- [(Int, VarMap r, r)]
xs]
            sv :: Solver r
sv = (State (Solver r) () -> Solver r -> Solver r)
-> Solver r -> State (Solver r) () -> Solver r
forall a b c. (a -> b -> c) -> b -> a -> c
flip State (Solver r) () -> Solver r -> Solver r
forall s a. State s a -> s -> s
execState (Node r -> Solver r
forall r. Node r -> Solver r
ndSolver Node r
node) (State (Solver r) () -> Solver r)
-> State (Solver r) () -> Solver r
forall a b. (a -> b) -> a -> b
$ do
                   Int
s <- LP r Int
forall r. LP r Int
newVar
                   let g :: Int -> r -> r
g Int
j r
x = Bool -> r -> r
forall a. HasCallStack => Bool -> a -> a
assert (r
a r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>= r
0) r
a
                         where
                           a :: r
a | Int
j Int -> VarSet -> Bool
`IS.member` VarSet
ivs =
                                if r -> r
forall a. RealFrac a => a -> a
fracPart r
x r -> r -> Bool
forall a. Ord a => a -> a -> Bool
<= r
f0
                                then r -> r
forall a. RealFrac a => a -> a
fracPart r
x
                                else (r
f0 r -> r -> r
forall a. Fractional a => a -> a -> a
/ (r
f0 r -> r -> r
forall a. Num a => a -> a -> a
- r
1)) r -> r -> r
forall a. Num a => a -> a -> a
* (r -> r
forall a. RealFrac a => a -> a
fracPart r
x r -> r -> r
forall a. Num a => a -> a -> a
- r
1)
                                     -- [Gomory1960] では (f0 / (1 - f0)) * (fracPart x - 1) としているが、
                                     -- これは誤り
                             | Bool
otherwise =
                                if r
x r -> r -> Bool
forall a. Ord a => a -> a -> Bool
>= r
0
                                then r
x
                                else (r
f0 r -> r -> r
forall a. Fractional a => a -> a -> a
/ (r
f0 r -> r -> r
forall a. Num a => a -> a -> a
- r
1)) r -> r -> r
forall a. Num a => a -> a -> a
* r
x
                   Tableau r -> State (Solver r) ()
forall r. Tableau r -> LP r ()
putTableau (Tableau r -> State (Solver r) ())
-> Tableau r -> State (Solver r) ()
forall a b. (a -> b) -> a -> b
$ Int -> Row r -> Tableau r -> Tableau r
forall a. Int -> a -> IntMap a -> IntMap a
IM.insert Int
s ((Int -> r -> r) -> VarMap r -> VarMap r
forall a b. (Int -> a -> b) -> IntMap a -> IntMap b
IM.mapWithKey (\Int
j r
x -> r -> r
forall a. Num a => a -> a
negate (Int -> r -> r
g Int
j r
x)) VarMap r
m0, r -> r
forall a. Num a => a -> a
negate r
f0) Tableau r
tbl
          in [Node r] -> Maybe [Node r]
forall a. a -> Maybe a
Just ([Node r] -> Maybe [Node r]) -> [Node r] -> Maybe [Node r]
forall a b. (a -> b) -> a -> b
$ [Node r
node{ ndSolver = sv2, ndDepth = ndDepth node + 1 } | Solver r
sv2 <- Maybe (Solver r) -> [Solver r]
forall a. Maybe a -> [a]
maybeToList (Solver r -> Maybe (Solver r)
reopt Solver r
sv)]
      | Bool
otherwise = -- branch
          let (Int
v0, r
val0) = (r, (Int, r)) -> (Int, r)
forall a b. (a, b) -> b
snd ((r, (Int, r)) -> (Int, r)) -> (r, (Int, r)) -> (Int, r)
forall a b. (a -> b) -> a -> b
$ ((r, (Int, r)) -> (r, (Int, r)) -> Ordering)
-> [(r, (Int, r))] -> (r, (Int, r))
forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (((r, (Int, r)) -> r) -> (r, (Int, r)) -> (r, (Int, r)) -> Ordering
forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing (r, (Int, r)) -> r
forall a b. (a, b) -> a
fst) [(r -> r
forall a. RealFrac a => a -> a
fracPart r
val, (Int
v, r
val)) | (Int
v,VarMap r
_,r
val) <- [(Int, VarMap r, r)]
xs]
              cs :: [OrdRel (Expr r)]
cs = [ Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
v0 Expr r -> Expr r -> OrdRel (Expr r)
forall e r. IsOrdRel e r => e -> e -> r
.>=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant (Integer -> r
forall a b. (Integral a, Num b) => a -> b
fromIntegral (r -> Integer
forall b. Integral b => r -> b
forall a b. (RealFrac a, Integral b) => a -> b
ceiling r
val0 :: Integer))
                   , Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
v0 Expr r -> Expr r -> OrdRel (Expr r)
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant (Integer -> r
forall a b. (Integral a, Num b) => a -> b
fromIntegral (r -> Integer
forall b. Integral b => r -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor r
val0 :: Integer))
                   ]
              svs :: [Solver r]
svs = [State (Solver r) () -> Solver r -> Solver r
forall s a. State s a -> s -> s
execState (OrdRel (Expr r) -> State (Solver r) ()
forall r. Real r => Atom r -> LP r ()
addConstraint OrdRel (Expr r)
c) (Node r -> Solver r
forall r. Node r -> Solver r
ndSolver Node r
node) | OrdRel (Expr r)
c <- [OrdRel (Expr r)]
cs]
          in [Node r] -> Maybe [Node r]
forall a. a -> Maybe a
Just ([Node r] -> Maybe [Node r]) -> [Node r] -> Maybe [Node r]
forall a b. (a -> b) -> a -> b
$ [Node r
node{ ndSolver = sv, ndDepth = ndDepth node + 1 } | Just Solver r
sv <- (Solver r -> Maybe (Solver r)) -> [Solver r] -> [Maybe (Solver r)]
forall a b. (a -> b) -> [a] -> [b]
map Solver r -> Maybe (Solver r)
reopt [Solver r]
svs]

      where
        tbl :: Simplex.Tableau r
        tbl :: Tableau r
tbl = Node r -> Tableau r
forall r. Node r -> Tableau r
ndTableau Node r
node

        xs :: [(Var, VarMap r, r)]
        xs :: [(Int, VarMap r, r)]
xs = [ (Int
v, VarMap r
m, r
val)
             | Int
v <- VarSet -> [Int]
IS.toList VarSet
ivs
             , Just (VarMap r
m, r
val) <- Maybe (Row r) -> [Maybe (Row r)]
forall a. a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> Tableau r -> Maybe (Row r)
forall a. Int -> IntMap a -> Maybe a
IM.lookup Int
v Tableau r
tbl)
             , Bool -> Bool
not (r -> Bool
forall a. RealFrac a => a -> Bool
isInteger r
val)
             ]

-- ---------------------------------------------------------------------------

example1 :: (Fractional r, Eq r) => (OptDir, LA.Expr r, [LA.Atom r], VarSet)
example1 :: forall r.
(Fractional r, Eq r) =>
(OptDir, Expr r, [Atom r], VarSet)
example1 = (OptDir
optdir, Expr r
obj, [Atom r]
cs, VarSet
ivs)
  where
    optdir :: OptDir
optdir = OptDir
OptMax
    x1 :: Expr r
x1 = Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
1
    x2 :: Expr r
x2 = Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
2
    x3 :: Expr r
x3 = Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
3
    x4 :: Expr r
x4 = Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var Int
4
    obj :: Expr r
obj = Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ r
Scalar (Expr r)
2 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ r
Scalar (Expr r)
3 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x3 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x4
    cs :: [Atom r]
cs =
      [ (-r
1) Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x3 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ r
Scalar (Expr r)
10Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^Expr r
x4 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
20
      , Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ r
Scalar (Expr r)
3 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x3 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
30
      , Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ r
Scalar (Expr r)
3.5 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x4 Expr r -> Expr r -> Atom r
forall e r. IsEqRel e r => e -> e -> r
.==. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x1
      , Expr r
x1 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
40
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x2
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x3
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
2 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x4
      , Expr r
x4 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
3
      ]
    ivs :: VarSet
ivs = Int -> VarSet
IS.singleton Int
4

_test1 :: Bool
_test1 :: Bool
_test1 = OptResult Rational
resultOptResult Rational -> OptResult Rational -> Bool
forall a. Eq a => a -> a -> Bool
==OptResult Rational
expected
  where
    (OptDir
optdir, Expr Rational
obj, [Atom Rational]
cs, VarSet
ivs) = (OptDir, Expr Rational, [Atom Rational], VarSet)
forall r.
(Fractional r, Eq r) =>
(OptDir, Expr r, [Atom r], VarSet)
example1
    result, expected :: OptResult Rational
    result :: OptResult Rational
result = OptDir
-> Expr Rational -> [Atom Rational] -> VarSet -> OptResult Rational
forall r.
RealFrac r =>
OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
optimize OptDir
optdir Expr Rational
obj [Atom Rational]
cs VarSet
ivs
    expected :: OptResult Rational
expected = Rational -> Model Rational -> OptResult Rational
forall r. r -> Model r -> OptResult r
Optimum (Rational
245Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2) ([(Int, Rational)] -> Model Rational
forall a. [(Int, a)] -> IntMap a
IM.fromList [(Int
1,Rational
40),(Int
2,Rational
21Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2),(Int
3,Rational
39Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2),(Int
4,Rational
3)])

_test1' :: Bool
_test1' :: Bool
_test1' = OptResult Rational
resultOptResult Rational -> OptResult Rational -> Bool
forall a. Eq a => a -> a -> Bool
==OptResult Rational
expected
  where
    (OptDir
optdir, Expr Rational
obj, [Atom Rational]
cs, VarSet
ivs) = (OptDir, Expr Rational, [Atom Rational], VarSet)
forall r.
(Fractional r, Eq r) =>
(OptDir, Expr r, [Atom r], VarSet)
example1
    f :: OptDir -> OptDir
f OptDir
OptMin = OptDir
OptMax
    f OptDir
OptMax = OptDir
OptMin
    result, expected :: OptResult Rational
    result :: OptResult Rational
result = OptDir
-> Expr Rational -> [Atom Rational] -> VarSet -> OptResult Rational
forall r.
RealFrac r =>
OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
optimize (OptDir -> OptDir
f OptDir
optdir) (Expr Rational -> Expr Rational
forall v. AdditiveGroup v => v -> v
negateV Expr Rational
obj) [Atom Rational]
cs VarSet
ivs
    expected :: OptResult Rational
expected = Rational -> Model Rational -> OptResult Rational
forall r. r -> Model r -> OptResult r
Optimum (-Rational
245Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2) ([(Int, Rational)] -> Model Rational
forall a. [(Int, a)] -> IntMap a
IM.fromList [(Int
1,Rational
40),(Int
2,Rational
21Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2),(Int
3,Rational
39Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2),(Int
4,Rational
3)])

-- 『数理計画法の基礎』(坂和 正敏) p.109 例 3.8
example2 :: (Fractional r, Eq r) => (OptDir, LA.Expr r, [LA.Atom r], VarSet)
example2 :: forall r.
(Fractional r, Eq r) =>
(OptDir, Expr r, [Atom r], VarSet)
example2 = (OptDir
optdir, Expr r
obj, [Atom r]
cs, VarSet
ivs)
  where
    optdir :: OptDir
optdir = OptDir
OptMin
    [Expr r
x1,Expr r
x2,Expr r
x3] = (Int -> Expr r) -> [Int] -> [Expr r]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Expr r
forall r. Num r => Int -> Expr r
LA.var [Int
1..Int
3]
    obj :: Expr r
obj = (-r
1) Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ r
Scalar (Expr r)
3 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^-^ r
Scalar (Expr r)
5 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x3
    cs :: [Atom r]
cs =
      [ r
Scalar (Expr r)
3 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ r
Scalar (Expr r)
4 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x2 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
10
      , r
Scalar (Expr r)
2 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x3 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
7
      , r
Scalar (Expr r)
3Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^Expr r
x1 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ Expr r
x2 Expr r -> Expr r -> Expr r
forall v. AdditiveGroup v => v -> v -> v
^+^ r
Scalar (Expr r)
4 Scalar (Expr r) -> Expr r -> Expr r
forall v. VectorSpace v => Scalar v -> v -> v
*^ Expr r
x3 Expr r -> Expr r -> Atom r
forall e r. IsEqRel e r => e -> e -> r
.==. r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
12
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x1
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x2
      , r -> Expr r
forall r. (Num r, Eq r) => r -> Expr r
LA.constant r
0 Expr r -> Expr r -> Atom r
forall e r. IsOrdRel e r => e -> e -> r
.<=. Expr r
x3
      ]
    ivs :: VarSet
ivs = [Int] -> VarSet
IS.fromList [Int
1,Int
2]

_test2 :: Bool
_test2 :: Bool
_test2 = OptResult Rational
result OptResult Rational -> OptResult Rational -> Bool
forall a. Eq a => a -> a -> Bool
== OptResult Rational
expected
  where
    result, expected :: OptResult Rational
    result :: OptResult Rational
result = OptDir
-> Expr Rational -> [Atom Rational] -> VarSet -> OptResult Rational
forall r.
RealFrac r =>
OptDir -> Expr r -> [Atom r] -> VarSet -> OptResult r
optimize OptDir
optdir Expr Rational
obj [Atom Rational]
cs VarSet
ivs
    expected :: OptResult Rational
expected = Rational -> Model Rational -> OptResult Rational
forall r. r -> Model r -> OptResult r
Optimum (-Rational
37Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2) ([(Int, Rational)] -> Model Rational
forall a. [(Int, a)] -> IntMap a
IM.fromList [(Int
1,Rational
0),(Int
2,Rational
2),(Int
3,Rational
5Rational -> Rational -> Rational
forall a. Fractional a => a -> a -> a
/Rational
2)])
    (OptDir
optdir, Expr Rational
obj, [Atom Rational]
cs, VarSet
ivs) = (OptDir, Expr Rational, [Atom Rational], VarSet)
forall r.
(Fractional r, Eq r) =>
(OptDir, Expr r, [Atom r], VarSet)
example2

-- ---------------------------------------------------------------------------