-----------------------------------------------------------------------------
-- |
-- 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
  forall (m :: * -> *) a. Monad m => a -> m a
return forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map forall a b. (RealFrac a, Integral b) => a -> b
round forall b c a. (b -> c) -> (a -> b) -> a -> c
. IntMap Rational -> IntMap Rational
mt forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> IntMap a -> IntMap b
IM.map forall a. Num a => Integer -> a
fromInteger 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 forall a. Eq a => a -> a -> Bool
== OptDir
OptMin then Expr Rational
obj else forall v. AdditiveGroup v => v -> v
negateV Expr Rational
obj, [Atom Rational]
cs)
    obj3 :: Expr Integer
obj3 = 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 = forall a b. (RealFrac a, Integral b) => a -> b
round forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational
cforall a. Num a => a -> a -> a
*)
        c :: Rational
c = forall a. Num a => Integer -> a
fromInteger forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' forall a. Integral a => a -> a -> a
lcm Integer
1 [forall a. Ratio a -> a
denominator Rational
c | (Rational
c,Var
_) <- forall r. Expr r -> [(r, Var)]
LA.terms Expr Rational
obj]
    cs3 :: [(Expr Integer, Integer)]
cs3 = forall a b. (a -> b) -> [a] -> [b]
map 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) = (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 = forall a b. (RealFrac a, Integral b) => a -> b
round forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational
cforall a. Num a => a -> a -> a
*)
        c :: Rational
c = forall a. Num a => Integer -> a
fromInteger forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' forall a. Integral a => a -> a -> a
lcm Integer
1 [forall a. Ratio a -> a
denominator Rational
c | (Rational
c,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
  | forall (t :: * -> *). Foldable t => t Bool -> Bool
or [Integer
c forall a. Ord a => a -> a -> Bool
< Integer
0 | (Integer
c,Var
x) <- forall r. Expr r -> [(r, Var)]
LA.terms Expr Integer
obj, Var
x forall a. Eq a => a -> a -> Bool
/= Var
LA.unitVar] = forall a. HasCallStack => [Char] -> a
error [Char]
"all coefficient of cost function should be non-negative"
  | Bool
otherwise =
  if forall a. IntMap a -> VarSet
IM.keysSet (forall a. (a -> Bool) -> IntMap a -> IntMap a
IM.filter (forall a. Eq a => a -> a -> Bool
/= Integer
0) Model Integer
m) VarSet -> VarSet -> Bool
`IS.isSubsetOf` VarSet
vs'
    then forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ 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 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' forall a. Num a => a -> a -> a
+ Var
1

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

    t :: Var
    t :: Var
t = Var
v2 forall a. Num a => a -> a -> a
+ 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) forall a. Monoid a => a -> a -> a
`mappend` VarSet -> MonomialOrder Var
elimOrdering (Var -> VarSet
IS.singleton Var
t) forall a. Monoid a => a -> a -> a
`mappend` Expr Integer -> MonomialOrder Var
costOrdering Expr Integer
obj forall a. Monoid a => a -> a -> a
`mappend` MonomialOrder Var
cmp

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

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

    [(Rational
_,Monomial Var
z)] = forall k v. Polynomial k v -> [Term k v]
P.terms (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]
vsforall a. [a] -> [a] -> [a]
++[Var]
vs2forall 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 = forall a. [(Var, a)] -> IntMap a
IM.fromDistinctAscList (forall k a. Map k a -> [(k, a)]
Map.toAscList (forall v. Monomial v -> Map v Integer
P.mindicesMap Monomial Var
xs)) forall a. IntMap a -> IntMap a -> IntMap a
`IM.union` 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 = forall a. Ord a => a -> a -> Ordering
compare forall b c a. (b -> b -> c) -> (a -> b) -> a -> a -> c
`on` Monomial Var -> Integer
f
  where
    vs :: VarSet
vs = forall a. Variables a => a -> VarSet
vars Expr Integer
obj
    f :: Monomial Var -> Integer
f Monomial Var
xs = 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 = forall a. Ord a => a -> a -> Ordering
compare 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 forall a b. (a -> b) -> a -> b
$ VarSet -> Bool
IS.null 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
_) <- forall k a. Map k a -> [(k, a)]
Map.toAscList forall a b. (a -> b) -> a -> b
$ forall v. Monomial v -> Map v Integer
P.mindicesMap Monomial Var
ys]