{-# OPTIONS_GHC -Wall #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.FourierMotzkin.Optimization
-- Copyright   :  (c) Masahiro Sakai 2014-2015
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- Naïve implementation of Fourier-Motzkin Variable Elimination
--
-- Reference:
--
-- * <http://users.cecs.anu.edu.au/~michaeln/pubs/arithmetic-dps.pdf>
--
-----------------------------------------------------------------------------
module ToySolver.Arith.FourierMotzkin.Optimization
  ( optimize
  ) where

import Control.Exception (assert)
import Control.Monad
import Data.Maybe
import qualified Data.IntMap as IntMap
import qualified Data.IntSet as IntSet
import Data.ExtendedReal
import qualified Data.Interval as Interval
import Data.Interval (Interval, (<=..<), (<..<=))
import Data.OptDir

import ToySolver.Data.DNF
import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar
import ToySolver.Arith.FourierMotzkin.Base

-- | @optimize dir obj φ@ returns @(I, lift)@ where
--
-- * @I@ is convex hull of feasible region, and
--
-- * @lift@ is a function, that takes @x ∈ I@ and returns the feasible solution with objective value /better than or equal to/ @x@.
--
-- Note:
--
-- * @'Interval.lowerBound' i@ (resp. @'Interval.upperBound' i@) is the optimal value in case of minimization (resp. maximization).
--
-- * If @I@ is empty, then the problem is INFEASIBLE.
--
optimize :: VarSet -> OptDir -> LA.Expr Rational -> [LA.Atom Rational] -> (Interval Rational, Rational -> Model Rational)
optimize :: VarSet
-> OptDir
-> Expr Rational
-> [Atom Rational]
-> (Interval Rational, Rational -> Model Rational)
optimize VarSet
vs OptDir
dir Expr Rational
obj [Atom Rational]
cs = (Interval Rational
ival, Rational -> Model Rational
m)
  where
    rs :: [(Interval Rational, Rational -> Model Rational)]
rs = [VarSet
-> Expr Rational
-> [Constr]
-> (Interval Rational, Rational -> Model Rational)
projectToObj' VarSet
vs Expr Rational
obj [Constr]
cs' | [Constr]
cs' <- forall lit. DNF lit -> [[lit]]
unDNF ([Atom Rational] -> DNF Constr
constraintsToDNF [Atom Rational]
cs)]
    ival :: Interval Rational
ival = forall r. Ord r => [Interval r] -> Interval r
Interval.hulls (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> a
fst [(Interval Rational, Rational -> Model Rational)]
rs)

    m :: Rational -> Model Rational
    m :: Rational -> Model Rational
m Rational
x = forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) (m :: * -> *) a.
(Foldable t, MonadPlus m) =>
t (m a) -> m a
msum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (Interval Rational, Rational -> Model Rational)
-> Maybe (Model Rational)
f [(Interval Rational, Rational -> Model Rational)]
rs
      where
        f :: (Interval Rational, Rational -> Model Rational) -> Maybe (Model Rational)
        f :: (Interval Rational, Rational -> Model Rational)
-> Maybe (Model Rational)
f (Interval Rational
i1,Rational -> Model Rational
m1) = do
          Rational
x' <- forall r. RealFrac r => Interval r -> Maybe Rational
Interval.simplestRationalWithin (forall r. Ord r => Interval r -> Interval r -> Interval r
Interval.intersection Interval Rational
i1 Interval Rational
ib)
          forall (m :: * -> *) a. Monad m => a -> m a
return (Rational -> Model Rational
m1 Rational
x')

        ib :: Interval Rational
ib = case OptDir
dir of
               OptDir
OptMin -> forall r. Extended r
NegInf forall r. Ord r => Extended r -> Extended r -> Interval r
<..<= forall r. r -> Extended r
Finite Rational
x
               OptDir
OptMax -> forall r. r -> Extended r
Finite Rational
x forall r. Ord r => Extended r -> Extended r -> Interval r
<=..< forall r. Extended r
PosInf

projectToObj' :: VarSet -> LA.Expr Rational -> [Constr] -> (Interval Rational, Rational -> Model Rational)
projectToObj' :: VarSet
-> Expr Rational
-> [Constr]
-> (Interval Rational, Rational -> Model Rational)
projectToObj' VarSet
vs Expr Rational
obj [Constr]
cs = VarSet
-> [Constr]
-> Var
-> (Interval Rational, Rational -> Model Rational)
projectTo' VarSet
vs (Rat -> Rat -> Constr
eqR (Expr Rational -> Rat
toRat (forall r. Num r => Var -> Expr r
LA.var Var
z)) (Expr Rational -> Rat
toRat Expr Rational
obj) forall a. a -> [a] -> [a]
: [Constr]
cs) Var
z
  where
    z :: Var
z = forall a. a -> Maybe a -> a
fromMaybe Var
0 (forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((forall a. Num a => a -> a -> a
+Var
1) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fst) (VarSet -> Maybe (Var, VarSet)
IntSet.maxView VarSet
vs))

projectTo' :: VarSet -> [Constr] -> Var -> (Interval Rational, Rational -> Model Rational)
projectTo' :: VarSet
-> [Constr]
-> Var
-> (Interval Rational, Rational -> Model Rational)
projectTo' VarSet
vs [Constr]
cs Var
z = forall a. a -> Maybe a -> a
fromMaybe (forall r. Ord r => Interval r
Interval.empty, \Rational
_ -> forall a. HasCallStack => [Char] -> a
error [Char]
"ToySolver.FourierMotzkin.projectTo': should not be called") forall a b. (a -> b) -> a -> b
$ do
  ([Constr]
ys,Model Rational -> Model Rational
mt) <- VarSet
-> [Constr] -> Maybe ([Constr], Model Rational -> Model Rational)
projectN' VarSet
vs forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< [Constr] -> Maybe [Constr]
simplify [Constr]
cs
  let (Bounds
bs,[Constr]
ws) = Var -> [Constr] -> (Bounds, [Constr])
collectBounds Var
z [Constr]
ys
  forall a. HasCallStack => Bool -> a -> a
assert (forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Constr]
ws) forall a b. (a -> b) -> a -> b
$ forall (m :: * -> *) a. Monad m => a -> m a
return ()
  let ival :: Interval Rational
ival = Model Rational -> Bounds -> Interval Rational
evalBounds forall a. IntMap a
IntMap.empty Bounds
bs
  forall (m :: * -> *) a. Monad m => a -> m a
return (Interval Rational
ival, \Rational
v -> Model Rational -> Model Rational
mt (forall a. Var -> a -> IntMap a
IntMap.singleton Var
z Rational
v))