-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.ContiTraverso
-- Copyright   :  (c) Masahiro Sakai 2012
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  portable
--
-- References:
--
-- * P. Conti and C. Traverso, "Buchberger algorithm and integer programming,"
--   Applied Algebra, Algebraic Algorithms and Error-Correcting Codes,
--   Lecture Notes in Computer Science Volume 539, 1991, pp 130-139
--   <https://doi.org/10.1007/3-540-54522-0_102>
--   <http://posso.dm.unipi.it/users/traverso/conti-traverso-ip.ps>
--
-- * IKEGAMI Daisuke, "数列と多項式の愛しい関係," 2011,
--   <http://madscientist.jp/~ikegami/articles/IntroSequencePolynomial.html>
--
-- * 伊藤雅史, , 平林 隆一, "整数計画問題のための b-Gröbner 基底変換アルゴリズム,"
--   <http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1295-27.pdf>
--
--
-----------------------------------------------------------------------------
module ToySolver.Arith.ContiTraverso
  ( solve
  , solve'
  ) where

import Data.Default.Class
import Data.Function
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import qualified Data.Map as Map
import Data.List
import Data.Monoid
import Data.Ratio
import Data.VectorSpace

import Data.OptDir

import ToySolver.Data.OrdRel
import qualified ToySolver.Data.LA as LA
import ToySolver.Data.Polynomial (Polynomial, UPolynomial, Monomial, MonomialOrder)
import qualified ToySolver.Data.Polynomial as P
import ToySolver.Data.Polynomial.GroebnerBasis as GB
import ToySolver.Data.IntVar
import qualified ToySolver.Arith.LPUtil as LPUtil

solve :: MonomialOrder Var -> VarSet -> OptDir -> LA.Expr Rational -> [LA.Atom Rational] -> Maybe (Model Integer)
solve :: MonomialOrder Var
-> VarSet
-> OptDir
-> Expr Rational
-> [Atom Rational]
-> Maybe (Model Integer)
solve MonomialOrder Var
cmp VarSet
vs OptDir
dir Expr Rational
obj [Atom Rational]
cs = do
  Model Integer
m <- MonomialOrder Var
-> VarSet
-> Expr Integer
-> [(Expr Integer, Integer)]
-> Maybe (Model Integer)
solve' MonomialOrder Var
cmp VarSet
vs Expr Integer
obj3 [(Expr Integer, Integer)]
cs3
  Model Integer -> Maybe (Model Integer)
forall a. a -> Maybe a
forall (m :: * -> *) a. Monad m => a -> m a
return (Model Integer -> Maybe (Model Integer))
-> (Model Integer -> Model Integer)
-> Model Integer
-> Maybe (Model Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Integer) -> IntMap Rational -> Model Integer
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (IntMap Rational -> Model Integer)
-> (Model Integer -> IntMap Rational)
-> Model Integer
-> Model Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IntMap Rational -> IntMap Rational
mt (IntMap Rational -> IntMap Rational)
-> (Model Integer -> IntMap Rational)
-> Model Integer
-> IntMap Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Rational) -> Model Integer -> IntMap Rational
forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Model Integer -> Maybe (Model Integer))
-> Model Integer -> Maybe (Model Integer)
forall a b. (a -> b) -> a -> b
$ Model Integer
m
  where
    ((Expr Rational
obj2,[(Expr Rational, Rational)]
cs2), IntMap Rational -> IntMap Rational
mt) = (Expr Rational, [Atom Rational])
-> ((Expr Rational, [(Expr Rational, Rational)]),
    IntMap Rational -> IntMap Rational)
LPUtil.toStandardForm (if OptDir
dir OptDir -> OptDir -> Bool
forall a. Eq a => a -> a -> Bool
== OptDir
OptMin then Expr Rational
obj else Expr Rational -> Expr Rational
forall v. AdditiveGroup v => v -> v
negateV Expr Rational
obj, [Atom Rational]
cs)
    obj3 :: Expr Integer
obj3 = (Rational -> Integer) -> Expr Rational -> Expr Integer
forall b a. (Num b, Eq b) => (a -> b) -> Expr a -> Expr b
LA.mapCoeff Rational -> Integer
g Expr Rational
obj2
      where
        g :: Rational -> Integer
g = Rational -> Integer
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Rational -> Rational) -> Rational -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational
cRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*)
        c :: Rational
c = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer -> Rational) -> Integer -> Rational
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
lcm Integer
1 [Rational -> Integer
forall a. Ratio a -> a
denominator Rational
c | (Rational
c,Var
_) <- Expr Rational -> [(Rational, Var)]
forall r. Expr r -> [(r, Var)]
LA.terms Expr Rational
obj]
    cs3 :: [(Expr Integer, Integer)]
cs3 = ((Expr Rational, Rational) -> (Expr Integer, Integer))
-> [(Expr Rational, Rational)] -> [(Expr Integer, Integer)]
forall a b. (a -> b) -> [a] -> [b]
map (Expr Rational, Rational) -> (Expr Integer, Integer)
forall {b}. Integral b => (Expr Rational, Rational) -> (Expr b, b)
f [(Expr Rational, Rational)]
cs2
    f :: (Expr Rational, Rational) -> (Expr b, b)
f (Expr Rational
lhs,Rational
rhs) = ((Rational -> b) -> Expr Rational -> Expr b
forall b a. (Num b, Eq b) => (a -> b) -> Expr a -> Expr b
LA.mapCoeff Rational -> b
g Expr Rational
lhs, Rational -> b
g Rational
rhs)
      where
        g :: Rational -> b
g = Rational -> b
forall b. Integral b => Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> b) -> (Rational -> Rational) -> Rational -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational
cRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*)
        c :: Rational
c = Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Integer -> Rational) -> Integer -> Rational
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> Integer) -> Integer -> [Integer] -> Integer
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
lcm Integer
1 [Rational -> Integer
forall a. Ratio a -> a
denominator Rational
c | (Rational
c,Var
_) <- Expr Rational -> [(Rational, Var)]
forall r. Expr r -> [(r, Var)]
LA.terms Expr Rational
lhs]

solve' :: MonomialOrder Var -> VarSet -> LA.Expr Integer -> [(LA.Expr Integer, Integer)] -> Maybe (Model Integer)
solve' :: MonomialOrder Var
-> VarSet
-> Expr Integer
-> [(Expr Integer, Integer)]
-> Maybe (Model Integer)
solve' MonomialOrder Var
cmp VarSet
vs' Expr Integer
obj [(Expr Integer, Integer)]
cs
  | [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
or [Integer
c Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 | (Integer
c,Var
x) <- Expr Integer -> [(Integer, Var)]
forall r. Expr r -> [(r, Var)]
LA.terms Expr Integer
obj, Var
x Var -> Var -> Bool
forall a. Eq a => a -> a -> Bool
/= Var
LA.unitVar] = [Char] -> Maybe (Model Integer)
forall a. HasCallStack => [Char] -> a
error [Char]
"all coefficient of cost function should be non-negative"
  | Bool
otherwise =
  if Model Integer -> VarSet
forall a. IntMap a -> VarSet
IM.keysSet ((Integer -> Bool) -> Model Integer -> Model Integer
forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0) Model Integer
m) VarSet -> VarSet -> Bool
`IS.isSubsetOf` VarSet
vs'
    then Model Integer -> Maybe (Model Integer)
forall a. a -> Maybe a
Just (Model Integer -> Maybe (Model Integer))
-> Model Integer -> Maybe (Model Integer)
forall a b. (a -> b) -> a -> b
$ (Var -> Integer -> Bool) -> Model Integer -> Model Integer
forall a. (Var -> a -> Bool) -> IntMap a -> IntMap a
IM.filterWithKey (\Var
y Integer
_ -> Var
y Var -> VarSet -> Bool
`IS.member` VarSet
vs') Model Integer
m
    else Maybe (Model Integer)
forall a. Maybe a
Nothing

  where
    vs :: [Var]
    vs :: [Var]
vs = VarSet -> [Var]
IS.toList VarSet
vs'

    v2 :: Var
    v2 :: Var
v2 = if VarSet -> Bool
IS.null VarSet
vs' then Var
0 else VarSet -> Var
IS.findMax VarSet
vs' Var -> Var -> Var
forall a. Num a => a -> a -> a
+ Var
1

    vs2 :: [Var]
    vs2 :: [Var]
vs2 = [Var
v2 .. Var
v2 Var -> Var -> Var
forall a. Num a => a -> a -> a
+ [(Expr Integer, Integer)] -> Var
forall a. [a] -> Var
forall (t :: * -> *) a. Foldable t => t a -> Var
length [(Expr Integer, Integer)]
cs Var -> Var -> Var
forall a. Num a => a -> a -> a
- Var
1]

    t :: Var
    t :: Var
t = Var
v2 Var -> Var -> Var
forall a. Num a => a -> a -> a
+ [(Expr Integer, Integer)] -> Var
forall a. [a] -> Var
forall (t :: * -> *) a. Foldable t => t a -> Var
length [(Expr Integer, Integer)]
cs

    cmp2 :: MonomialOrder Var
    cmp2 :: MonomialOrder Var
cmp2 = VarSet -> MonomialOrder Var
elimOrdering ([Var] -> VarSet
IS.fromList [Var]
vs2) MonomialOrder Var -> MonomialOrder Var -> MonomialOrder Var
forall a. Monoid a => a -> a -> a
`mappend` VarSet -> MonomialOrder Var
elimOrdering (Var -> VarSet
IS.singleton Var
t) MonomialOrder Var -> MonomialOrder Var -> MonomialOrder Var
forall a. Monoid a => a -> a -> a
`mappend` Expr Integer -> MonomialOrder Var
costOrdering Expr Integer
obj MonomialOrder Var -> MonomialOrder Var -> MonomialOrder Var
forall a. Monoid a => a -> a -> a
`mappend` MonomialOrder Var
cmp

    gb :: [Polynomial Rational Var]
    gb :: [Polynomial Rational Var]
gb = Options
-> MonomialOrder Var
-> [Polynomial Rational Var]
-> [Polynomial Rational Var]
forall k v.
(Eq k, Fractional k, Ord k, Ord v) =>
Options -> MonomialOrder v -> [Polynomial k v] -> [Polynomial k v]
GB.basis' Options
forall a. Default a => a
def MonomialOrder Var
cmp2 ([Polynomial Rational Var] -> Polynomial Rational Var
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ((Var -> Polynomial Rational Var)
-> [Var] -> [Polynomial Rational Var]
forall a b. (a -> b) -> [a] -> [b]
map Var -> Polynomial Rational Var
forall a v. Var a v => v -> a
P.var (Var
tVar -> [Var] -> [Var]
forall a. a -> [a] -> [a]
:[Var]
vs2)) Polynomial Rational Var
-> Polynomial Rational Var -> Polynomial Rational Var
forall a. Num a => a -> a -> a
- Polynomial Rational Var
1 Polynomial Rational Var
-> [Polynomial Rational Var] -> [Polynomial Rational Var]
forall a. a -> [a] -> [a]
: [Polynomial Rational Var]
phi)
      where
        phi :: [Polynomial Rational Var]
phi = do
          Var
xj <- [Var]
vs
          let aj :: [(Var, Integer)]
aj = [(Var
yi, Integer
aij) | (Var
yi,(Expr Integer
ai,Integer
_)) <- [Var]
-> [(Expr Integer, Integer)] -> [(Var, (Expr Integer, Integer))]
forall a b. [a] -> [b] -> [(a, b)]
zip [Var]
vs2 [(Expr Integer, Integer)]
cs, let aij :: Integer
aij = Var -> Expr Integer -> Integer
forall r. Num r => Var -> Expr r -> r
LA.coeff Var
xj Expr Integer
ai]
          Polynomial Rational Var -> [Polynomial Rational Var]
forall a. a -> [a]
forall (m :: * -> *) a. Monad m => a -> m a
return (Polynomial Rational Var -> [Polynomial Rational Var])
-> Polynomial Rational Var -> [Polynomial Rational Var]
forall a b. (a -> b) -> a -> b
$  [Polynomial Rational Var] -> Polynomial Rational Var
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Var -> Polynomial Rational Var
forall a v. Var a v => v -> a
P.var Var
yi Polynomial Rational Var -> Integer -> Polynomial Rational Var
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
aij    | (Var
yi, Integer
aij) <- [(Var, Integer)]
aj, Integer
aij Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0]
                  Polynomial Rational Var
-> Polynomial Rational Var -> Polynomial Rational Var
forall a. Num a => a -> a -> a
- [Polynomial Rational Var] -> Polynomial Rational Var
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Var -> Polynomial Rational Var
forall a v. Var a v => v -> a
P.var Var
yi Polynomial Rational Var -> Integer -> Polynomial Rational Var
forall a b. (Num a, Integral b) => a -> b -> a
^ (-Integer
aij) | (Var
yi, Integer
aij) <- [(Var, Integer)]
aj, Integer
aij Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0] Polynomial Rational Var
-> Polynomial Rational Var -> Polynomial Rational Var
forall a. Num a => a -> a -> a
* Var -> Polynomial Rational Var
forall a v. Var a v => v -> a
P.var Var
xj

    yb :: Polynomial Rational Var
yb = [Polynomial Rational Var] -> Polynomial Rational Var
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Var -> Polynomial Rational Var
forall a v. Var a v => v -> a
P.var Var
yi Polynomial Rational Var -> Integer -> Polynomial Rational Var
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
bi | ((Expr Integer
_,Integer
bi),Var
yi) <- [(Expr Integer, Integer)]
-> [Var] -> [((Expr Integer, Integer), Var)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(Expr Integer, Integer)]
cs [Var]
vs2]

    [(Rational
_,Monomial Var
z)] = Polynomial Rational Var -> [Term Rational Var]
forall k v. Polynomial k v -> [Term k v]
P.terms (MonomialOrder Var
-> Polynomial Rational Var
-> [Polynomial Rational Var]
-> Polynomial Rational Var
forall k v.
(Eq k, Fractional k, Ord v) =>
MonomialOrder v
-> Polynomial k v -> [Polynomial k v] -> Polynomial k v
P.reduce MonomialOrder Var
cmp2 Polynomial Rational Var
yb [Polynomial Rational Var]
gb)

    m :: Model Integer
m = [Var] -> Monomial Var -> Model Integer
mkModel ([Var]
vs[Var] -> [Var] -> [Var]
forall a. [a] -> [a] -> [a]
++[Var]
vs2[Var] -> [Var] -> [Var]
forall a. [a] -> [a] -> [a]
++[Var
t]) Monomial Var
z

mkModel :: [Var] -> Monomial Var -> Model Integer
mkModel :: [Var] -> Monomial Var -> Model Integer
mkModel [Var]
vs Monomial Var
xs = [(Var, Integer)] -> Model Integer
forall a. [(Var, a)] -> IntMap a
IM.fromDistinctAscList (Map Var Integer -> [(Var, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toAscList (Monomial Var -> Map Var Integer
forall v. Monomial v -> Map v Integer
P.mindicesMap Monomial Var
xs)) Model Integer -> Model Integer -> Model Integer
forall a. IntMap a -> IntMap a -> IntMap a
`IM.union` [(Var, Integer)] -> Model Integer
forall a. [(Var, a)] -> IntMap a
IM.fromList [(Var
x, Integer
0) | Var
x <- [Var]
vs]
-- IM.union is left-biased

costOrdering :: LA.Expr Integer -> MonomialOrder Var
costOrdering :: Expr Integer -> MonomialOrder Var
costOrdering Expr Integer
obj = Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Integer -> Ordering)
-> (Monomial Var -> Integer) -> MonomialOrder Var
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial Var -> Integer
f
  where
    vs :: VarSet
vs = Expr Integer -> VarSet
forall a. Variables a => a -> VarSet
vars Expr Integer
obj
    f :: Monomial Var -> Integer
f Monomial Var
xs = Model Integer -> Expr Integer -> Integer
forall m e v. Eval m e v => m -> e -> v
LA.eval ([Var] -> Monomial Var -> Model Integer
mkModel (VarSet -> [Var]
IS.toList VarSet
vs) Monomial Var
xs) Expr Integer
obj

elimOrdering :: IS.IntSet -> MonomialOrder Var
elimOrdering :: VarSet -> MonomialOrder Var
elimOrdering VarSet
xs = Bool -> Bool -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Bool -> Bool -> Ordering)
-> (Monomial Var -> Bool) -> MonomialOrder Var
forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial Var -> Bool
f
  where
    f :: Monomial Var -> Bool
f Monomial Var
ys = Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ VarSet -> Bool
IS.null (VarSet -> Bool) -> VarSet -> Bool
forall a b. (a -> b) -> a -> b
$ VarSet
xs VarSet -> VarSet -> VarSet
`IS.intersection` VarSet
ys'
      where
        ys' :: VarSet
ys' = [Var] -> VarSet
IS.fromDistinctAscList [Var
y | (Var
y,Integer
_) <- Map Var Integer -> [(Var, Integer)]
forall k a. Map k a -> [(k, a)]
Map.toAscList (Map Var Integer -> [(Var, Integer)])
-> Map Var Integer -> [(Var, Integer)]
forall a b. (a -> b) -> a -> b
$ Monomial Var -> Map Var Integer
forall v. Monomial v -> Map v Integer
P.mindicesMap Monomial Var
ys]