{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.MIP
-- Copyright   :  (c) Masahiro Sakai 2012
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-portable
--
-- Naïve implementation of MIP solver based on Simplex module
--
-- Reference:
--
-- * <http://www.math.cuhk.edu.hk/~wei/lpch3.pdf>
--
-- * 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>
--
-- * 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>
--
-- * R. C. Daniel and Martyn Jeffreys.
--   \"Unboundedness in Integer and Discrete Programming L.P. Relaxations\"
--   The Journal of the Operational Research Society, Vol. 30, No. 12. (1979)
--   <http://www.jstor.org/stable/3009435>
--
-----------------------------------------------------------------------------
module ToySolver.Arith.MIP
  (
  -- * The @Solver@ type
    Solver
  , newSolver

  -- * Solving
  , optimize

  -- * Extract results
  , getBestSolution
  , getBestValue
  , getBestModel

  -- * Configulation
  , setNThread
  , setLogger
  , setOnUpdateBestSolution
  , setShowRational
  ) where

import Prelude hiding (log)

import Control.Monad
import Control.Exception
import Control.Concurrent
import Control.Concurrent.STM
import Data.Default.Class
import Data.List
import Data.OptDir
import Data.Ord
import Data.IORef
import Data.Maybe
import qualified Data.IntSet as IS
import qualified Data.IntMap as IM
import qualified Data.Map as Map
import qualified Data.Sequence as Seq
import qualified Data.Foldable as F
import Data.VectorSpace
import System.Clock
import System.Timeout
import Text.Printf

import qualified ToySolver.Data.LA as LA
import ToySolver.Data.OrdRel ((.<=.), (.>=.))
import qualified ToySolver.Arith.Simplex as Simplex
import ToySolver.Arith.Simplex (OptResult (..), Var, Model)
import ToySolver.Internal.Util (isInteger, fracPart)

data Solver
  = MIP
  { Solver -> Solver
mipRootLP :: Simplex.Solver
  , Solver -> IntSet
mipIVs    :: IS.IntSet
  , Solver -> TVar (Maybe Node)
mipBest   :: TVar (Maybe Node)

  , Solver -> IORef Int
mipNThread :: IORef Int
  , Solver -> IORef (Maybe (String -> IO ()))
mipLogger  :: IORef (Maybe (String -> IO ()))
  , Solver -> IORef (Model -> Rational -> IO ())
mipOnUpdateBestSolution :: IORef (Model -> Rational -> IO ())
  , Solver -> IORef Bool
mipShowRational :: IORef Bool
  }

data Node =
  Node
  { Node -> Solver
ndLP    :: Simplex.Solver
  , Node -> Int
ndDepth :: {-# UNPACK #-} !Int
  , Node -> Rational
ndValue :: Rational
  }

newSolver :: Simplex.Solver -> IS.IntSet -> IO Solver
newSolver :: Solver -> IntSet -> IO Solver
newSolver Solver
lp IntSet
ivs = do
  Solver
lp2 <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> m (GenericSolverM m v)
Simplex.cloneSolver Solver
lp

  forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ (IntSet -> [Int]
IS.toList IntSet
ivs) forall a b. (a -> b) -> a -> b
$ \Int
v -> do
    Bound Rational
lb <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getLB Solver
lp2 Int
v
    case Bound Rational
lb of
      Just (Rational
l,IntSet
_) | Bool -> Bool
not (forall a. RealFrac a => a -> Bool
isInteger Rational
l) ->
        forall (m :: * -> *) v.
(PrimMonad m, SolverValue v) =>
GenericSolverM m v -> Int -> v -> m ()
Simplex.assertLower Solver
lp2 Int
v (forall a. Num a => Integer -> a
fromInteger (forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
l))
      Bound Rational
_ -> forall (m :: * -> *) a. Monad m => a -> m a
return ()
    Bound Rational
ub <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getUB Solver
lp2 Int
v
    case Bound Rational
ub of
      Just (Rational
u,IntSet
_) | Bool -> Bool
not (forall a. RealFrac a => a -> Bool
isInteger Rational
u) ->
        forall (m :: * -> *) v.
(PrimMonad m, SolverValue v) =>
GenericSolverM m v -> Int -> v -> m ()
Simplex.assertLower Solver
lp2 Int
v (forall a. Num a => Integer -> a
fromInteger (forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
u))
      Bound Rational
_ -> forall (m :: * -> *) a. Monad m => a -> m a
return ()

  TVar (Maybe Node)
bestRef <- forall a. a -> IO (TVar a)
newTVarIO forall a. Maybe a
Nothing

  IORef Int
nthreadRef <- forall a. a -> IO (IORef a)
newIORef Int
0
  IORef (Maybe (String -> IO ()))
logRef  <- forall a. a -> IO (IORef a)
newIORef forall a. Maybe a
Nothing
  IORef Bool
showRef <- forall a. a -> IO (IORef a)
newIORef Bool
False
  IORef (Model -> Rational -> IO ())
updateRef <- forall a. a -> IO (IORef a)
newIORef (\Model
_ Rational
_ -> forall (m :: * -> *) a. Monad m => a -> m a
return ())

  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$
    MIP
    { mipRootLP :: Solver
mipRootLP = Solver
lp2
    , mipIVs :: IntSet
mipIVs    = IntSet
ivs
    , mipBest :: TVar (Maybe Node)
mipBest   = TVar (Maybe Node)
bestRef

    , mipNThread :: IORef Int
mipNThread = IORef Int
nthreadRef
    , mipLogger :: IORef (Maybe (String -> IO ()))
mipLogger = IORef (Maybe (String -> IO ()))
logRef
    , mipOnUpdateBestSolution :: IORef (Model -> Rational -> IO ())
mipOnUpdateBestSolution = IORef (Model -> Rational -> IO ())
updateRef
    , mipShowRational :: IORef Bool
mipShowRational = IORef Bool
showRef
    }

optimize :: Solver -> IO OptResult
optimize :: Solver -> IO OptResult
optimize Solver
solver = do
  let lp :: Solver
lp = Solver -> Solver
mipRootLP Solver
solver
  Model -> Rational -> IO ()
update <- forall a. IORef a -> IO a
readIORef (Solver -> IORef (Model -> Rational -> IO ())
mipOnUpdateBestSolution Solver
solver)
  Solver -> String -> IO ()
log Solver
solver String
"MIP: Solving LP relaxation..."
  Bool
ret <- forall (m :: * -> *) v.
(PrimMonad m, SolverValue v) =>
GenericSolverM m v -> m Bool
Simplex.check Solver
lp
  if Bool -> Bool
not Bool
ret
  then forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unsat
  else do
    String
s0 <- Solver -> Rational -> IO String
showValue Solver
solver forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall (m :: * -> *) v. PrimMonad m => GenericSolverM m v -> m v
Simplex.getObjValue Solver
lp
    Solver -> String -> IO ()
log Solver
solver (forall r. PrintfType r => String -> r
printf String
"MIP: LP relaxation is satisfiable with obj = %s" String
s0)
    Solver -> String -> IO ()
log Solver
solver String
"MIP: Optimizing LP relaxation"
    OptResult
ret2 <- forall (m :: * -> *).
PrimMonad m =>
GenericSolverM m Rational -> Options -> m OptResult
Simplex.optimize Solver
lp forall a. Default a => a
def
    case OptResult
ret2 of
      OptResult
Unsat    -> forall a. (?callStack::CallStack) => String -> a
error String
"should not happen"
      OptResult
ObjLimit -> forall a. (?callStack::CallStack) => String -> a
error String
"should not happen"
      OptResult
Unbounded -> do
        Solver -> String -> IO ()
log Solver
solver String
"MIP: LP relaxation is unbounded"
        let ivs :: IntSet
ivs = Solver -> IntSet
mipIVs Solver
solver
        if IntSet -> Bool
IS.null IntSet
ivs
          then forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unbounded
          else do
            {-
              * 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.
            -}
            Expr Rational
origObj <- forall (m :: * -> *) v.
(PrimMonad m, SolverValue v) =>
GenericSolverM m v -> m (Expr Rational)
Simplex.getObj Solver
lp
            Solver
lp2 <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> m (GenericSolverM m v)
Simplex.cloneSolver Solver
lp
            forall (m :: * -> *) v. PrimMonad m => GenericSolverM m v -> m ()
Simplex.clearLogger Solver
lp2
            forall (m :: * -> *) v.
(PrimMonad m, SolverValue v) =>
GenericSolverM m v -> Expr Rational -> m ()
Simplex.setObj Solver
lp2 (forall r. (Num r, Eq r) => r -> Expr r
LA.constant Rational
0)
            Solver -> Solver -> (Model -> Rational -> IO ()) -> IO ()
branchAndBound Solver
solver Solver
lp2 forall a b. (a -> b) -> a -> b
$ \Model
m Rational
_ -> do
              Model -> Rational -> IO ()
update Model
m (forall m e v. Eval m e v => m -> e -> v
LA.eval Model
m Expr Rational
origObj)
            Maybe Node
best <- forall a. TVar a -> IO a
readTVarIO (Solver -> TVar (Maybe Node)
mipBest Solver
solver)
            case Maybe Node
best of
              Just Node
nd -> do
                Model
m <- forall v (m :: * -> *).
(SolverValue v, PrimMonad m) =>
GenericSolverM m v -> m Model
Simplex.getModel (Node -> Solver
ndLP Node
nd)
                forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ forall a. TVar a -> a -> STM ()
writeTVar (Solver -> TVar (Maybe Node)
mipBest Solver
solver) forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just Node
nd{ ndValue :: Rational
ndValue = forall m e v. Eval m e v => m -> e -> v
LA.eval Model
m Expr Rational
origObj }
                forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unbounded
              Maybe Node
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unsat
      OptResult
Optimum -> do
        String
s1 <- Solver -> Rational -> IO String
showValue Solver
solver forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< forall (m :: * -> *) v. PrimMonad m => GenericSolverM m v -> m v
Simplex.getObjValue Solver
lp
        Solver -> String -> IO ()
log Solver
solver forall a b. (a -> b) -> a -> b
$ String
"MIP: LP relaxation optimum is " forall a. [a] -> [a] -> [a]
++ String
s1
        Solver -> String -> IO ()
log Solver
solver String
"MIP: Integer optimization begins..."
        forall (m :: * -> *) v. PrimMonad m => GenericSolverM m v -> m ()
Simplex.clearLogger Solver
lp
        Solver -> Solver -> (Model -> Rational -> IO ()) -> IO ()
branchAndBound Solver
solver Solver
lp Model -> Rational -> IO ()
update
        Maybe Node
m <- forall a. TVar a -> IO a
readTVarIO (Solver -> TVar (Maybe Node)
mipBest Solver
solver)
        case Maybe Node
m of
          Maybe Node
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Unsat
          Just Node
_ -> forall (m :: * -> *) a. Monad m => a -> m a
return OptResult
Optimum

branchAndBound :: Solver -> Simplex.Solver -> (Model -> Rational -> IO ()) -> IO ()
branchAndBound :: Solver -> Solver -> (Model -> Rational -> IO ()) -> IO ()
branchAndBound Solver
solver Solver
rootLP Model -> Rational -> IO ()
update = do
  OptDir
dir <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> m OptDir
Simplex.getOptDir Solver
rootLP
  Rational
rootVal <- forall (m :: * -> *) v. PrimMonad m => GenericSolverM m v -> m v
Simplex.getObjValue Solver
rootLP
  let root :: Node
root = Node{ ndLP :: Solver
ndLP = Solver
rootLP, ndDepth :: Int
ndDepth = Int
0, ndValue :: Rational
ndValue = Rational
rootVal }

  TVar (Seq Node)
pool <- forall a. a -> IO (TVar a)
newTVarIO (forall a. a -> Seq a
Seq.singleton Node
root)
  TVar (Map ThreadId Node)
activeThreads <- forall a. a -> IO (TVar a)
newTVarIO (forall k a. Map k a
Map.empty)
  TVar Int
visitedNodes <- forall a. a -> IO (TVar a)
newTVarIO Int
0
  TChan Node
solchan <- forall a. IO (TChan a)
newTChanIO

  let addNode :: Node -> STM ()
      addNode :: Node -> STM ()
addNode Node
nd = do
        forall a. TVar a -> (a -> a) -> STM ()
modifyTVar TVar (Seq Node)
pool (forall a. Seq a -> a -> Seq a
Seq.|> Node
nd)

      pickNode :: IO (Maybe Node)
      pickNode :: IO (Maybe Node)
pickNode = do
        ThreadId
self <- IO ThreadId
myThreadId
        forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ forall a. TVar a -> (a -> a) -> STM ()
modifyTVar TVar (Map ThreadId Node)
activeThreads (forall k a. Ord k => k -> Map k a -> Map k a
Map.delete ThreadId
self)
        forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ do
          Seq Node
s <- forall a. TVar a -> STM a
readTVar TVar (Seq Node)
pool
          case forall a. Seq a -> ViewL a
Seq.viewl Seq Node
s of
            Node
nd Seq.:< Seq Node
s2 -> do
              forall a. TVar a -> a -> STM ()
writeTVar TVar (Seq Node)
pool Seq Node
s2
              forall a. TVar a -> (a -> a) -> STM ()
modifyTVar TVar (Map ThreadId Node)
activeThreads (forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert ThreadId
self Node
nd)
              forall (m :: * -> *) a. Monad m => a -> m a
return (forall a. a -> Maybe a
Just Node
nd)
            ViewL Node
Seq.EmptyL -> do
              Map ThreadId Node
ths <- forall a. TVar a -> STM a
readTVar TVar (Map ThreadId Node)
activeThreads
              if forall k a. Map k a -> Bool
Map.null Map ThreadId Node
ths
                then forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
                else forall a. STM a
retry

      processNode :: Node -> IO ()
      processNode :: Node -> IO ()
processNode Node
node = do
        let lp :: Solver
lp = Node -> Solver
ndLP Node
node
        Maybe Rational
lim <- forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Node -> Rational
ndValue) forall a b. (a -> b) -> a -> b
$ forall a. TVar a -> IO a
readTVarIO (Solver -> TVar (Maybe Node)
mipBest Solver
solver)
        OptResult
ret <- forall (m :: * -> *).
PrimMonad m =>
GenericSolverM m Rational -> Options -> m OptResult
Simplex.dualSimplex Solver
lp forall a. Default a => a
def{ objLimit :: Maybe Rational
Simplex.objLimit = Maybe Rational
lim }

        case OptResult
ret of
          OptResult
Unbounded -> forall a. (?callStack::CallStack) => String -> a
error String
"should not happen"
          OptResult
Unsat ->  forall (m :: * -> *) a. Monad m => a -> m a
return ()
          OptResult
ObjLimit -> forall (m :: * -> *) a. Monad m => a -> m a
return ()
          OptResult
Optimum -> do
            Rational
val <- forall (m :: * -> *) v. PrimMonad m => GenericSolverM m v -> m v
Simplex.getObjValue Solver
lp
            Bool
p <- Solver -> Rational -> IO Bool
prune Solver
solver Rational
val
            forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
unless Bool
p forall a b. (a -> b) -> a -> b
$ do
              [(Int, Rational)]
xs <- Node -> IntSet -> IO [(Int, Rational)]
violated Node
node (Solver -> IntSet
mipIVs Solver
solver)
              case [(Int, Rational)]
xs of
                [] -> forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ forall a. TChan a -> a -> STM ()
writeTChan TChan Node
solchan forall a b. (a -> b) -> a -> b
$ Node
node { ndValue :: Rational
ndValue = Rational
val }
                [(Int, Rational)]
_ -> do
                  Maybe Int
r <- if Node -> Int
ndDepth Node
node forall a. Integral a => a -> a -> a
`mod` Int
100 forall a. Eq a => a -> a -> Bool
/= Int
0
                       then forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
                       else forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM forall a. [a] -> Maybe a
listToMaybe forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a.
Applicative m =>
(a -> m Bool) -> [a] -> m [a]
filterM (Solver -> Int -> IO Bool
canDeriveGomoryCut Solver
lp) forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> a
fst [(Int, Rational)]
xs
                  case Maybe Int
r of
                    Maybe Int
Nothing -> do -- branch
                      let (Int
v0,Rational
val0) = forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a.
Foldable t =>
(a -> a -> Ordering) -> t a -> a
maximumBy (forall a b. Ord a => (b -> a) -> b -> b -> Ordering
comparing forall a b. (a, b) -> b
snd)
                                      [((Int
v,Rational
vval), forall a. Num a => a -> a
abs (forall a. Num a => Integer -> a
fromInteger (forall a b. (RealFrac a, Integral b) => a -> b
round Rational
vval) forall a. Num a => a -> a -> a
- Rational
vval)) | (Int
v,Rational
vval) <- [(Int, Rational)]
xs]
                      let lp1 :: Solver
lp1 = Solver
lp
                      Solver
lp2 <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> m (GenericSolverM m v)
Simplex.cloneSolver Solver
lp
                      forall (m :: * -> *).
PrimMonad m =>
GenericSolverM m Rational -> Atom Rational -> m ()
Simplex.assertAtom Solver
lp1 (forall r. Num r => Int -> Expr r
LA.var Int
v0 forall e r. IsOrdRel e r => e -> e -> r
.<=. forall r. (Num r, Eq r) => r -> Expr r
LA.constant (forall a. Num a => Integer -> a
fromInteger (forall a b. (RealFrac a, Integral b) => a -> b
floor Rational
val0)))
                      forall (m :: * -> *).
PrimMonad m =>
GenericSolverM m Rational -> Atom Rational -> m ()
Simplex.assertAtom Solver
lp2 (forall r. Num r => Int -> Expr r
LA.var Int
v0 forall e r. IsOrdRel e r => e -> e -> r
.>=. forall r. (Num r, Eq r) => r -> Expr r
LA.constant (forall a. Num a => Integer -> a
fromInteger (forall a b. (RealFrac a, Integral b) => a -> b
ceiling Rational
val0)))
                      forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ do
                        Node -> STM ()
addNode forall a b. (a -> b) -> a -> b
$ Solver -> Int -> Rational -> Node
Node Solver
lp1 (Node -> Int
ndDepth Node
node forall a. Num a => a -> a -> a
+ Int
1) Rational
val
                        Node -> STM ()
addNode forall a b. (a -> b) -> a -> b
$ Solver -> Int -> Rational -> Node
Node Solver
lp2 (Node -> Int
ndDepth Node
node forall a. Num a => a -> a -> a
+ Int
1) Rational
val
                        forall a. TVar a -> (a -> a) -> STM ()
modifyTVar TVar Int
visitedNodes (forall a. Num a => a -> a -> a
+Int
1)
                    Just Int
v -> do -- cut
                      Atom Rational
atom <- Solver -> IntSet -> Int -> IO (Atom Rational)
deriveGomoryCut Solver
lp (Solver -> IntSet
mipIVs Solver
solver) Int
v
                      forall (m :: * -> *).
PrimMonad m =>
GenericSolverM m Rational -> Atom Rational -> m ()
Simplex.assertAtom Solver
lp Atom Rational
atom
                      forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ do
                        Node -> STM ()
addNode forall a b. (a -> b) -> a -> b
$ Solver -> Int -> Rational -> Node
Node Solver
lp (Node -> Int
ndDepth Node
node forall a. Num a => a -> a -> a
+ Int
1) Rational
val

  let isCompleted :: STM Bool
isCompleted = do
        Seq Node
nodes <- forall a. TVar a -> STM a
readTVar TVar (Seq Node)
pool
        Map ThreadId Node
threads <- forall a. TVar a -> STM a
readTVar TVar (Map ThreadId Node)
activeThreads
        forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Seq a -> Bool
Seq.null Seq Node
nodes Bool -> Bool -> Bool
&& forall k a. Map k a -> Bool
Map.null Map ThreadId Node
threads

  -- fork worker threads
  Int
nthreads <- forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM (forall a. Ord a => a -> a -> a
max Int
1) forall a b. (a -> b) -> a -> b
$ forall a. IORef a -> IO a
readIORef (Solver -> IORef Int
mipNThread Solver
solver)

  Solver -> String -> IO ()
log Solver
solver forall a b. (a -> b) -> a -> b
$ forall r. PrintfType r => String -> r
printf String
"MIP: forking %d worker threads..." Int
nthreads

  TimeSpec
startCPU <- Clock -> IO TimeSpec
getTime Clock
ProcessCPUTime
  TimeSpec
startWC  <- Clock -> IO TimeSpec
getTime Clock
Monotonic
  TMVar SomeException
ex <- forall a. IO (TMVar a)
newEmptyTMVarIO

  let printStatus :: Seq.Seq Node -> Int -> IO ()
      printStatus :: Seq Node -> Int -> IO ()
printStatus Seq Node
nodes Int
visited
        | forall a. Seq a -> Bool
Seq.null Seq Node
nodes = forall (m :: * -> *) a. Monad m => a -> m a
return () -- should not happen
        | Bool
otherwise = do
            TimeSpec
nowCPU <- Clock -> IO TimeSpec
getTime Clock
ProcessCPUTime
            TimeSpec
nowWC  <- Clock -> IO TimeSpec
getTime Clock
Monotonic
            let spentCPU :: Int64
spentCPU = TimeSpec -> Int64
sec (TimeSpec
nowCPU TimeSpec -> TimeSpec -> TimeSpec
`diffTimeSpec` TimeSpec
startCPU)
            let spentWC :: Int64
spentWC  = TimeSpec -> Int64
sec (TimeSpec
nowWC  TimeSpec -> TimeSpec -> TimeSpec
`diffTimeSpec` TimeSpec
startWC )

            let vs :: [Rational]
vs = forall a b. (a -> b) -> [a] -> [b]
map Node -> Rational
ndValue (forall (t :: * -> *) a. Foldable t => t a -> [a]
F.toList Seq Node
nodes)
                dualBound :: Rational
dualBound =
                  case OptDir
dir of
                    OptDir
OptMin -> forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum [Rational]
vs
                    OptDir
OptMax -> forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum [Rational]
vs

            Maybe Rational
primalBound <- do
              Maybe Node
x <- forall a. TVar a -> IO a
readTVarIO (Solver -> TVar (Maybe Node)
mipBest Solver
solver) -- TODO: 引数にするようにした方が良い?
              forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ case Maybe Node
x of
                Maybe Node
Nothing -> forall a. Maybe a
Nothing
                Just Node
node -> forall a. a -> Maybe a
Just (Node -> Rational
ndValue Node
node)

            (String
p,String
g) <- case Maybe Rational
primalBound of
                   Maybe Rational
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return (String
"not yet found", String
"--")
                   Just Rational
val -> do
                     String
p <- Solver -> Rational -> IO String
showValue Solver
solver Rational
val
                     let g :: String
g = if Rational
val forall a. Eq a => a -> a -> Bool
== Rational
0
                             then String
"inf"
                             else forall r. PrintfType r => String -> r
printf String
"%.2f%%" (forall a. Fractional a => Rational -> a
fromRational (forall a. Num a => a -> a
abs (Rational
dualBound forall a. Num a => a -> a -> a
- Rational
val) forall a. Num a => a -> a -> a
* Rational
100 forall a. Fractional a => a -> a -> a
/ forall a. Num a => a -> a
abs Rational
val) :: Double)
                     forall (m :: * -> *) a. Monad m => a -> m a
return (String
p, String
g)
            String
d <- Solver -> Rational -> IO String
showValue Solver
solver Rational
dualBound

            let range :: String
range =
                  case OptDir
dir of
                    OptDir
OptMin -> String
p forall a. [a] -> [a] -> [a]
++ String
" >= " forall a. [a] -> [a] -> [a]
++ String
d
                    OptDir
OptMax -> String
p forall a. [a] -> [a] -> [a]
++ String
" <= " forall a. [a] -> [a] -> [a]
++ String
d

            Solver -> String -> IO ()
log Solver
solver forall a b. (a -> b) -> a -> b
$ forall r. PrintfType r => String -> r
printf String
"cpu time = %d sec; wc time = %d sec; active nodes = %d; visited nodes = %d; %s; gap = %s"
              Int64
spentCPU Int64
spentWC (forall a. Seq a -> Int
Seq.length Seq Node
nodes) Int
visited String
range String
g

  forall b. ((forall a. IO a -> IO a) -> IO b) -> IO b
mask forall a b. (a -> b) -> a -> b
$ \(forall a. IO a -> IO a
restore :: forall a. IO a -> IO a) -> do
    [ThreadId]
threads <- forall (m :: * -> *) a. Applicative m => Int -> m a -> m [a]
replicateM Int
nthreads forall a b. (a -> b) -> a -> b
$ do
      IO () -> IO ThreadId
forkIO forall a b. (a -> b) -> a -> b
$ do
        let loop :: IO ()
loop = do
              Maybe Node
m <- IO (Maybe Node)
pickNode
              case Maybe Node
m of
                Maybe Node
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return ()
                Just Node
node -> Node -> IO ()
processNode Node
node forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> IO ()
loop
        Either SomeException ()
ret <- forall e a. Exception e => IO a -> IO (Either e a)
try forall a b. (a -> b) -> a -> b
$ forall a. IO a -> IO a
restore IO ()
loop
        case Either SomeException ()
ret of
          Left SomeException
e -> forall a. STM a -> IO a
atomically (forall a. TMVar a -> a -> STM ()
putTMVar TMVar SomeException
ex SomeException
e)
          Right ()
_ -> forall (m :: * -> *) a. Monad m => a -> m a
return ()

    let propagateException :: SomeException -> IO ()
        propagateException :: SomeException -> IO ()
propagateException SomeException
e = do
          forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
(a -> m b) -> t a -> m ()
mapM_ (\ThreadId
t -> forall e. Exception e => ThreadId -> e -> IO ()
throwTo ThreadId
t SomeException
e) [ThreadId]
threads
          forall e a. Exception e => e -> IO a
throwIO SomeException
e

    let loop :: IO ()
loop = do
          Either SomeException (Maybe (IO ()))
ret <- forall e a. Exception e => IO a -> IO (Either e a)
try forall a b. (a -> b) -> a -> b
$ forall a. Int -> IO a -> IO (Maybe a)
timeout (Int
2forall a. Num a => a -> a -> a
*Int
1000forall a. Num a => a -> a -> a
*Int
1000) forall a b. (a -> b) -> a -> b
$ forall a. IO a -> IO a
restore forall a b. (a -> b) -> a -> b
$ forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) (m :: * -> *) a.
(Foldable t, MonadPlus m) =>
t (m a) -> m a
msum
            [ do Node
node <- forall a. TChan a -> STM a
readTChan TChan Node
solchan
                 Bool
ret <- do
                   Maybe Node
old <- forall a. TVar a -> STM a
readTVar (Solver -> TVar (Maybe Node)
mipBest Solver
solver)
                   case Maybe Node
old of
                     Maybe Node
Nothing -> do
                       forall a. TVar a -> a -> STM ()
writeTVar (Solver -> TVar (Maybe Node)
mipBest Solver
solver) (forall a. a -> Maybe a
Just Node
node)
                       forall (m :: * -> *) a. Monad m => a -> m a
return Bool
True
                     Just Node
best -> do
                       let isBetter :: Bool
isBetter = if OptDir
dirforall a. Eq a => a -> a -> Bool
==OptDir
OptMin then Node -> Rational
ndValue Node
node forall a. Ord a => a -> a -> Bool
< Node -> Rational
ndValue Node
best else Node -> Rational
ndValue Node
node forall a. Ord a => a -> a -> Bool
> Node -> Rational
ndValue Node
best
                       forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
isBetter forall a b. (a -> b) -> a -> b
$ forall a. TVar a -> a -> STM ()
writeTVar (Solver -> TVar (Maybe Node)
mipBest Solver
solver) (forall a. a -> Maybe a
Just Node
node)
                       forall (m :: * -> *) a. Monad m => a -> m a
return Bool
isBetter
                 forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ do
                   forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when Bool
ret forall a b. (a -> b) -> a -> b
$ do
                     let lp :: Solver
lp = Node -> Solver
ndLP Node
node
                     Model
m <- forall v (m :: * -> *).
(SolverValue v, PrimMonad m) =>
GenericSolverM m v -> m Model
Simplex.getModel Solver
lp
                     Model -> Rational -> IO ()
update Model
m (Node -> Rational
ndValue Node
node)
                   IO ()
loop
            , do Bool
b <- STM Bool
isCompleted
                 forall (f :: * -> *). Alternative f => Bool -> f ()
guard Bool
b
                 forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. Monad m => a -> m a
return ()
            , do SomeException
e <- forall a. TMVar a -> STM a
readTMVar TMVar SomeException
ex
                 forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ SomeException -> IO ()
propagateException SomeException
e
            ]

          case Either SomeException (Maybe (IO ()))
ret of
            Left (SomeException
e::SomeException) -> SomeException -> IO ()
propagateException SomeException
e
            Right (Just IO ()
m) -> IO ()
m
            Right Maybe (IO ())
Nothing -> do -- timeout
              (Seq Node
nodes, Int
visited) <- forall a. STM a -> IO a
atomically forall a b. (a -> b) -> a -> b
$ do
                Seq Node
nodes    <- forall a. TVar a -> STM a
readTVar TVar (Seq Node)
pool
                Map ThreadId Node
athreads <- forall a. TVar a -> STM a
readTVar TVar (Map ThreadId Node)
activeThreads
                Int
visited  <- forall a. TVar a -> STM a
readTVar TVar Int
visitedNodes
                forall (m :: * -> *) a. Monad m => a -> m a
return (forall a. [a] -> Seq a
Seq.fromList (forall k a. Map k a -> [a]
Map.elems Map ThreadId Node
athreads) forall a. Seq a -> Seq a -> Seq a
Seq.>< Seq Node
nodes, Int
visited)
              Seq Node -> Int -> IO ()
printStatus Seq Node
nodes Int
visited
              IO ()
loop

    IO ()
loop

getBestSolution :: Solver -> IO (Maybe (Model, Rational))
getBestSolution :: Solver -> IO (Maybe (Model, Rational))
getBestSolution Solver
solver = do
  Maybe Node
ret <- forall a. TVar a -> IO a
readTVarIO (Solver -> TVar (Maybe Node)
mipBest Solver
solver)
  case Maybe Node
ret of
    Maybe Node
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
    Just Node
node -> do
      Model
m <- forall v (m :: * -> *).
(SolverValue v, PrimMonad m) =>
GenericSolverM m v -> m Model
Simplex.getModel (Node -> Solver
ndLP Node
node)
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just (Model
m, Node -> Rational
ndValue Node
node)

getBestModel :: Solver -> IO (Maybe Model)
getBestModel :: Solver -> IO (Maybe Model)
getBestModel Solver
solver = forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a b. (a, b) -> a
fst) forall a b. (a -> b) -> a -> b
$ Solver -> IO (Maybe (Model, Rational))
getBestSolution Solver
solver

getBestValue :: Solver -> IO (Maybe Rational)
getBestValue :: Solver -> IO (Maybe Rational)
getBestValue Solver
solver = forall (m :: * -> *) a1 r. Monad m => (a1 -> r) -> m a1 -> m r
liftM (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a b. (a, b) -> b
snd) forall a b. (a -> b) -> a -> b
$ Solver -> IO (Maybe (Model, Rational))
getBestSolution Solver
solver

violated :: Node -> IS.IntSet -> IO [(Var, Rational)]
violated :: Node -> IntSet -> IO [(Int, Rational)]
violated Node
node IntSet
ivs = do
  Model
m <- forall v (m :: * -> *).
(SolverValue v, PrimMonad m) =>
GenericSolverM m v -> m Model
Simplex.getModel (Node -> Solver
ndLP Node
node)
  let p :: (Int, a) -> Bool
p (Int
v,a
val) = Int
v Int -> IntSet -> Bool
`IS.member` IntSet
ivs Bool -> Bool -> Bool
&& Bool -> Bool
not (forall a. RealFrac a => a -> Bool
isInteger a
val)
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
filter forall {a}. RealFrac a => (Int, a) -> Bool
p (forall a. IntMap a -> [(Int, a)]
IM.toList Model
m)

prune :: Solver -> Rational -> IO Bool
prune :: Solver -> Rational -> IO Bool
prune Solver
solver Rational
lb = do
  Maybe Node
b <- forall a. TVar a -> IO a
readTVarIO (Solver -> TVar (Maybe Node)
mipBest Solver
solver)
  case Maybe Node
b of
    Maybe Node
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return Bool
False
    Just Node
node -> do
      OptDir
dir <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> m OptDir
Simplex.getOptDir (Solver -> Solver
mipRootLP Solver
solver)
      forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ if OptDir
dirforall a. Eq a => a -> a -> Bool
==OptDir
OptMin then Node -> Rational
ndValue Node
node forall a. Ord a => a -> a -> Bool
<= Rational
lb else Node -> Rational
ndValue Node
node forall a. Ord a => a -> a -> Bool
>= Rational
lb

showValue :: Solver -> Rational -> IO String
showValue :: Solver -> Rational -> IO String
showValue Solver
solver Rational
v = do
  Bool
printRat <- forall a. IORef a -> IO a
readIORef (Solver -> IORef Bool
mipShowRational Solver
solver)
  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall v. SolverValue v => Bool -> v -> String
Simplex.showValue Bool
printRat Rational
v

setShowRational :: Solver -> Bool -> IO ()
setShowRational :: Solver -> Bool -> IO ()
setShowRational Solver
solver = forall a. IORef a -> a -> IO ()
writeIORef (Solver -> IORef Bool
mipShowRational Solver
solver)

setNThread :: Solver -> Int -> IO ()
setNThread :: Solver -> Int -> IO ()
setNThread Solver
solver = forall a. IORef a -> a -> IO ()
writeIORef (Solver -> IORef Int
mipNThread Solver
solver)

{--------------------------------------------------------------------
  Logging
--------------------------------------------------------------------}

-- | set callback function for receiving messages.
setLogger :: Solver -> (String -> IO ()) -> IO ()
setLogger :: Solver -> (String -> IO ()) -> IO ()
setLogger Solver
solver String -> IO ()
logger = do
  forall a. IORef a -> a -> IO ()
writeIORef (Solver -> IORef (Maybe (String -> IO ()))
mipLogger Solver
solver) (forall a. a -> Maybe a
Just String -> IO ()
logger)

setOnUpdateBestSolution :: Solver -> (Model -> Rational -> IO ()) -> IO ()
setOnUpdateBestSolution :: Solver -> (Model -> Rational -> IO ()) -> IO ()
setOnUpdateBestSolution Solver
solver Model -> Rational -> IO ()
cb = do
  forall a. IORef a -> a -> IO ()
writeIORef (Solver -> IORef (Model -> Rational -> IO ())
mipOnUpdateBestSolution Solver
solver) Model -> Rational -> IO ()
cb

log :: Solver -> String -> IO ()
log :: Solver -> String -> IO ()
log Solver
solver String
msg = Solver -> IO String -> IO ()
logIO Solver
solver (forall (m :: * -> *) a. Monad m => a -> m a
return String
msg)

logIO :: Solver -> IO String -> IO ()
logIO :: Solver -> IO String -> IO ()
logIO Solver
solver IO String
action = do
  Maybe (String -> IO ())
m <- forall a. IORef a -> IO a
readIORef (Solver -> IORef (Maybe (String -> IO ()))
mipLogger Solver
solver)
  case Maybe (String -> IO ())
m of
    Maybe (String -> IO ())
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return ()
    Just String -> IO ()
logger -> IO String
action forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= String -> IO ()
logger

{--------------------------------------------------------------------
  GomoryCut
--------------------------------------------------------------------}

deriveGomoryCut :: Simplex.Solver -> IS.IntSet -> Var -> IO (LA.Atom Rational)
deriveGomoryCut :: Solver -> IntSet -> Int -> IO (Atom Rational)
deriveGomoryCut Solver
lp IntSet
ivs Int
xi = do
  Rational
v0 <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m v
Simplex.getValue Solver
lp Int
xi
  let f0 :: Rational
f0 = forall a. RealFrac a => a -> a
fracPart Rational
v0
  forall a. (?callStack::CallStack) => Bool -> a -> a
assert (Rational
0 forall a. Ord a => a -> a -> Bool
< Rational
f0 Bool -> Bool -> Bool
&& Rational
f0 forall a. Ord a => a -> a -> Bool
< Rational
1) forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. Monad m => a -> m a
return ()

  Expr Rational
row <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Expr Rational)
Simplex.getRow Solver
lp Int
xi

  -- remove fixed variables
  let p :: (a, Int) -> IO Bool
p (a
_,Int
xj) = do
        Bound Rational
lb <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getLB Solver
lp Int
xj
        Bound Rational
ub <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getUB Solver
lp Int
xj
        case (Bound Rational
lb,Bound Rational
ub) of
          (Just (Rational, IntSet)
l, Just (Rational, IntSet)
u) -> forall (m :: * -> *) a. Monad m => a -> m a
return ((Rational, IntSet)
l forall a. Ord a => a -> a -> Bool
< (Rational, IntSet)
u)
          (Bound Rational, Bound Rational)
_ -> forall (m :: * -> *) a. Monad m => a -> m a
return Bool
True
  [(Rational, Int)]
ns <- forall (m :: * -> *) a.
Applicative m =>
(a -> m Bool) -> [a] -> m [a]
filterM forall {a}. (a, Int) -> IO Bool
p forall a b. (a -> b) -> a -> b
$ forall r. Expr r -> [(r, Int)]
LA.terms Expr Rational
row

  [(Rational, Int)]
js <- forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (m :: * -> *) a.
Applicative m =>
(a -> m Bool) -> [a] -> m [a]
filterM [(Rational, Int)]
ns forall a b. (a -> b) -> a -> b
$ \(Rational
_, Int
xj) -> do
    Rational
vj <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m v
Simplex.getValue Solver
lp Int
xj
    Bound Rational
lb <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getLB Solver
lp Int
xj
    forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just Rational
vj forall a. Eq a => a -> a -> Bool
== forall v. SolverValue v => Bound v -> Maybe v
Simplex.boundValue Bound Rational
lb
  [(Rational, Int)]
ks <- forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (m :: * -> *) a.
Applicative m =>
(a -> m Bool) -> [a] -> m [a]
filterM [(Rational, Int)]
ns forall a b. (a -> b) -> a -> b
$ \(Rational
_, Int
xj) -> do
    Rational
vj <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m v
Simplex.getValue Solver
lp Int
xj
    Bound Rational
ub <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getUB Solver
lp Int
xj
    forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just Rational
vj forall a. Eq a => a -> a -> Bool
== forall v. SolverValue v => Bound v -> Maybe v
Simplex.boundValue Bound Rational
ub

  [Expr Rational]
xs1 <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [(Rational, Int)]
js forall a b. (a -> b) -> a -> b
$ \(Rational
aij, Int
xj) -> do
    let fj :: Rational
fj = forall a. RealFrac a => a -> a
fracPart Rational
aij
    Just (Rational
lj,IntSet
_) <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getLB Solver
lp Int
xj
    let c :: Rational
c = if Int
xj Int -> IntSet -> Bool
`IS.member` IntSet
ivs
            then (if Rational
fj forall a. Ord a => a -> a -> Bool
<= Rational
1 forall a. Num a => a -> a -> a
- Rational
f0 then Rational
fj  forall a. Fractional a => a -> a -> a
/ (Rational
1 forall a. Num a => a -> a -> a
- Rational
f0) else ((Rational
1 forall a. Num a => a -> a -> a
- Rational
fj) forall a. Fractional a => a -> a -> a
/ Rational
f0))
            else (if Rational
aij forall a. Ord a => a -> a -> Bool
> Rational
0      then Rational
aij forall a. Fractional a => a -> a -> a
/ (Rational
1 forall a. Num a => a -> a -> a
- Rational
f0) else (-Rational
aij     forall a. Fractional a => a -> a -> a
/ Rational
f0))
    forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Rational
c forall v. VectorSpace v => Scalar v -> v -> v
*^ (forall r. Num r => Int -> Expr r
LA.var Int
xj forall v. AdditiveGroup v => v -> v -> v
^-^ forall r. (Num r, Eq r) => r -> Expr r
LA.constant Rational
lj)
  [Expr Rational]
xs2 <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [(Rational, Int)]
ks forall a b. (a -> b) -> a -> b
$ \(Rational
aij, Int
xj) -> do
    let fj :: Rational
fj = forall a. RealFrac a => a -> a
fracPart Rational
aij
    Just (Rational
uj, IntSet
_) <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getUB Solver
lp Int
xj
    let c :: Rational
c = if Int
xj Int -> IntSet -> Bool
`IS.member` IntSet
ivs
            then (if Rational
fj forall a. Ord a => a -> a -> Bool
<= Rational
f0 then Rational
fj  forall a. Fractional a => a -> a -> a
/ Rational
f0 else ((Rational
1 forall a. Num a => a -> a -> a
- Rational
fj) forall a. Fractional a => a -> a -> a
/ (Rational
1 forall a. Num a => a -> a -> a
- Rational
f0)))
            else (if Rational
aij forall a. Ord a => a -> a -> Bool
> Rational
0  then Rational
aij forall a. Fractional a => a -> a -> a
/ Rational
f0 else (-Rational
aij     forall a. Fractional a => a -> a -> a
/ (Rational
1 forall a. Num a => a -> a -> a
- Rational
f0)))
    forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Rational
c forall v. VectorSpace v => Scalar v -> v -> v
*^ (forall r. (Num r, Eq r) => r -> Expr r
LA.constant Rational
uj forall v. AdditiveGroup v => v -> v -> v
^-^ forall r. Num r => Int -> Expr r
LA.var Int
xj)

  forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [Expr Rational]
xs1 forall v. AdditiveGroup v => v -> v -> v
^+^ forall (f :: * -> *) v. (Foldable f, AdditiveGroup v) => f v -> v
sumV [Expr Rational]
xs2 forall e r. IsOrdRel e r => e -> e -> r
.>=. forall r. (Num r, Eq r) => r -> Expr r
LA.constant Rational
1

-- TODO: Simplexをδに対応させたら、xi, xj がδを含まない有理数であるという条件も必要
canDeriveGomoryCut :: Simplex.Solver -> Var -> IO Bool
canDeriveGomoryCut :: Solver -> Int -> IO Bool
canDeriveGomoryCut Solver
lp Int
xi = do
  Bool
b <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m Bool
Simplex.isBasicVariable Solver
lp Int
xi
  if Bool -> Bool
not Bool
b
    then forall (m :: * -> *) a. Monad m => a -> m a
return Bool
False
    else do
      Rational
val <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m v
Simplex.getValue Solver
lp Int
xi
      if forall a. RealFrac a => a -> Bool
isInteger Rational
val
        then forall (m :: * -> *) a. Monad m => a -> m a
return Bool
False
        else do
          Expr Rational
row <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Expr Rational)
Simplex.getRow Solver
lp Int
xi
          [Bool]
ys <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM (forall r. Expr r -> [(r, Int)]
LA.terms Expr Rational
row) forall a b. (a -> b) -> a -> b
$ \(Rational
_,Int
xj) -> do
            Rational
vj <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m v
Simplex.getValue Solver
lp Int
xj
            Bound Rational
lb <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getLB Solver
lp Int
xj
            Bound Rational
ub <- forall (m :: * -> *) v.
PrimMonad m =>
GenericSolverM m v -> Int -> m (Bound v)
Simplex.getUB Solver
lp Int
xj
            forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just Rational
vj forall a. Eq a => a -> a -> Bool
== forall v. SolverValue v => Bound v -> Maybe v
Simplex.boundValue Bound Rational
lb Bool -> Bool -> Bool
|| forall a. a -> Maybe a
Just Rational
vj forall a. Eq a => a -> a -> Bool
== forall v. SolverValue v => Bound v -> Maybe v
Simplex.boundValue Bound Rational
ub
          forall (m :: * -> *) a. Monad m => a -> m a
return (forall (t :: * -> *). Foldable t => t Bool -> Bool
and [Bool]
ys)