msakai / toysolver / 767

Committed 22 May 2025 12:57PM UTC coverage: 71.847% (-0.06%) from 71.906%

Build # 767

Build Type

push

github

Committed by

web-flow

Commit Message

Merge f4d92f6d1 into c10d256d2

Run Details

11104 of 15455 relevant lines covered (71.85%)

0.72 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

80.49

/src/ToySolver/Arith/FourierMotzkin/Base.hs

{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_HADDOCK show-extensions #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  ToySolver.Arith.FourierMotzkin.Base
-- Copyright   :  (c) Masahiro Sakai 2011-2013
-- License     :  BSD-style
--
-- Maintainer  :  masahiro.sakai@gmail.com
-- Stability   :  provisional
-- Portability :  non-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.Base
    (
    -- * Primitive constraints
      Constr (..)
    , eqR
    , ExprZ
    , fromLAAtom
    , toLAAtom
    , constraintsToDNF
    , simplify

    -- * Bounds
    , Bounds
    , evalBounds
    , boundsToConstrs
    , collectBounds

    -- * Projection
    , project
    , project'
    , projectN
    , projectN'

    -- * Solving
    , solve
    , solve'

    -- * Utilities used by other modules
    , Rat
    , toRat
    ) where

import Control.Monad
import Data.List (foldl', foldl1')
import Data.Maybe
import Data.Ratio
import qualified Data.IntMap as IM
import qualified Data.IntSet as IS
import Data.VectorSpace hiding (project)
import qualified Data.Interval as Interval
import Data.Interval (Interval, Extended (..), (<=..<), (<..<=), (<..<))

import ToySolver.Data.OrdRel
import ToySolver.Data.Boolean
import ToySolver.Data.DNF
import qualified ToySolver.Data.LA as LA
import ToySolver.Data.IntVar

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

type ExprZ = LA.Expr Integer

-- | (t,c) represents t/c, and c must be >0.
type Rat = (ExprZ, Integer)

toRat :: LA.Expr Rational -> Rat
toRat e = seq m $ (LA.mapCoeff f e, m)
  where
    f x = numerator (fromInteger m * x)
    m = foldl' lcm 1 [denominator c | (c,_) <- LA.terms e]

fromRat :: Rat -> LA.Expr Rational
fromRat (e,c) = LA.mapCoeff (% c) e

evalRat :: Model Rational -> Rat -> Rational
evalRat model (e, d) = LA.lift1 1 (model IM.!) (LA.mapCoeff fromIntegral e) / (fromIntegral d)

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

-- | Atomic constraint
data Constr
  = IsNonneg ExprZ
  -- ^ e ≥ 0
  | IsPos ExprZ
  -- ^ e > 0
  | IsZero ExprZ
  -- ^ e = 0
  deriving (Show, Eq, Ord)

instance Variables Constr where
  vars (IsPos t) = vars t
  vars (IsNonneg t) = vars t
  vars (IsZero t) = vars t

leR, ltR, geR, gtR, eqR :: Rat -> Rat -> Constr
leR (e1,c) (e2,d) = IsNonneg $ normalizeExpr $ c *^ e2 ^-^ d *^ e1
ltR (e1,c) (e2,d) = IsPos $ normalizeExpr $ c *^ e2 ^-^ d *^ e1
geR = flip leR
gtR = flip ltR
eqR (e1,c) (e2,d) = IsZero $ normalizeExpr $ c *^ e2 ^-^ d *^ e1

normalizeExpr :: ExprZ -> ExprZ
normalizeExpr e = LA.mapCoeff (`div` d) e
  where d = abs $ gcd' $ map fst $ LA.terms e

-- "subst1Constr x t c" computes c[t/x]
subst1Constr :: Var -> LA.Expr Rational -> Constr -> Constr
subst1Constr x t c =
  case c of
    IsPos e    -> IsPos (f e)
    IsNonneg e -> IsNonneg (f e)
    IsZero e   -> IsZero (f e)
  where
    f :: ExprZ -> ExprZ
    f = normalizeExpr . fst . toRat . LA.applySubst1 x t . LA.mapCoeff fromInteger

-- | Simplify conjunction of 'Constr's.
-- It returns 'Nothing' when a inconsistency is detected.
simplify :: [Constr] -> Maybe [Constr]
simplify = fmap concat . mapM f
  where
    f :: Constr -> Maybe [Constr]
    f c@(IsPos e) =
      case LA.asConst e of
        Just x -> guard (x > 0) >> return []
        Nothing -> return [c]
    f c@(IsNonneg e) =
      case LA.asConst e of
        Just x -> guard (x >= 0) >> return []
        Nothing -> return [c]
    f c@(IsZero e) =
      case LA.asConst e of
        Just x -> guard (x == 0) >> return []
        Nothing -> return [c]

instance Eval (Model Rational) Constr Bool where
  eval m (IsPos t)    = LA.eval m (LA.mapCoeff fromInteger t :: LA.Expr Rational) > 0
  eval m (IsNonneg t) = LA.eval m (LA.mapCoeff fromInteger t :: LA.Expr Rational) >= 0
  eval m (IsZero t)   = LA.eval m (LA.mapCoeff fromInteger t :: LA.Expr Rational) == 0

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

fromLAAtom :: LA.Atom Rational -> DNF Constr
fromLAAtom (OrdRel a op b) = atomR' op (toRat a) (toRat b)
  where
    atomR' :: RelOp -> Rat -> Rat -> DNF Constr
    atomR' op a b =
      case op of
        Le -> DNF [[a `leR` b]]
        Lt -> DNF [[a `ltR` b]]
        Ge -> DNF [[a `geR` b]]
        Gt -> DNF [[a `gtR` b]]
        Eql -> DNF [[a `eqR` b]]
        NEq -> DNF [[a `ltR` b], [a `gtR` b]]

toLAAtom :: Constr -> LA.Atom Rational
toLAAtom (IsNonneg e) = LA.mapCoeff fromInteger e .>=. LA.constant 0
toLAAtom (IsPos e)    = LA.mapCoeff fromInteger e .>.  LA.constant 0
toLAAtom (IsZero e)   = LA.mapCoeff fromInteger e .==. LA.constant 0

constraintsToDNF :: [LA.Atom Rational] -> DNF Constr
constraintsToDNF = andB . map fromLAAtom

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

{-
(ls1,ls2,us1,us2) represents
{ x | ∀(M,c)∈ls1. M/c≤x, ∀(M,c)∈ls2. M/c<x, ∀(M,c)∈us1. x≤M/c, ∀(M,c)∈us2. x<M/c }
-}
type Bounds = ([Rat], [Rat], [Rat], [Rat])

evalBounds :: Model Rational -> Bounds -> Interval Rational
evalBounds model (ls1,ls2,us1,us2) =
  Interval.intersections $
    [ Finite (evalRat model x) <=..< PosInf | x <- ls1 ] ++
    [ Finite (evalRat model x) <..<  PosInf | x <- ls2 ] ++
    [ NegInf <..<= Finite (evalRat model x) | x <- us1 ] ++
    [ NegInf <..<  Finite (evalRat model x) | x <- us2 ]

boundsToConstrs :: Bounds -> Maybe [Constr]
boundsToConstrs  (ls1, ls2, us1, us2) = simplify $
  [ x `leR` y | x <- ls1, y <- us1 ] ++
  [ x `ltR` y | x <- ls1, y <- us2 ] ++
  [ x `ltR` y | x <- ls2, y <- us1 ] ++
  [ x `ltR` y | x <- ls2, y <- us2 ]

collectBounds :: Var -> [Constr] -> (Bounds, [Constr])
collectBounds v = foldr phi (([],[],[],[]),[])
  where
    phi :: Constr -> (Bounds, [Constr]) -> (Bounds, [Constr])
    phi constr@(IsNonneg t) x = f False constr t x
    phi constr@(IsPos t) x = f True constr t x
    phi constr@(IsZero t) (bnd@(ls1,ls2,us1,us2), xs) =
      case LA.extractMaybe v t of
        Nothing -> (bnd, constr : xs)
        Just (c,t') -> ((t'' : ls1, ls2, t'' : us1, us2), xs)
          where
            t'' = (signum c *^ negateV t', abs c)

    f :: Bool -> Constr -> ExprZ -> (Bounds, [Constr]) -> (Bounds, [Constr])
    f strict constr t (bnd@(ls1,ls2,us1,us2), xs) =
      case LA.extract v t of
        (c,t') ->
          case c `compare` 0 of
            EQ -> (bnd, constr : xs)
            GT ->
              if strict
              then ((ls1, (negateV t', c) : ls2, us1, us2), xs) -- 0 < cx + M ⇔ -M/c <  x
              else (((negateV t', c) : ls1, ls2, us1, us2), xs) -- 0 ≤ cx + M ⇔ -M/c ≤ x
            LT ->
              if strict
              then ((ls1, ls2, us1, (t', negate c) : us2), xs) -- 0 < cx + M ⇔ x < M/-c
              else ((ls1, ls2, (t', negate c) : us1, us2), xs) -- 0 ≤ cx + M ⇔ x ≤ M/-c

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

{-| @'project' x φ@ returns @[(ψ_1, lift_1), …, (ψ_m, lift_m)]@ such that:

* @⊢_LRA ∀y1, …, yn. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_m)@ where @{y1, …, yn} = FV(φ) \\ {x}@, and

* if @M ⊧_LRA ψ_i@ then @lift_i M ⊧_LRA φ@.
-}
project :: Var -> [LA.Atom Rational] -> [([LA.Atom Rational], Model Rational -> Model Rational)]
project v xs = do
  ys <- unDNF $ constraintsToDNF xs
  (zs, mt) <- maybeToList $ project' v ys
  return (map toLAAtom zs, mt)

project' :: Var -> [Constr] -> Maybe ([Constr], Model Rational -> Model Rational)
project' v cs = projectN' (IS.singleton v) cs

{-| @'projectN' {x1,…,xm} φ@ returns @[(ψ_1, lift_1), …, (ψ_n, lift_n)]@ such that:

* @⊢_LRA ∀y1, …, yp. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_n)@ where @{y1, …, yp} = FV(φ) \\ {x1,…,xm}@, and

* if @M ⊧_LRA ψ_i@ then @lift_i M ⊧_LRA φ@.
-}
projectN :: VarSet -> [LA.Atom Rational] -> [([LA.Atom Rational], Model Rational -> Model Rational)]
projectN vs xs = do
  ys <- unDNF $ constraintsToDNF xs
  (zs, mt) <- maybeToList $ projectN' vs ys
  return (map toLAAtom zs, mt)

projectN' :: VarSet -> [Constr] -> Maybe ([Constr], Model Rational -> Model Rational)
projectN' = f
  where
    f vs xs
      | IS.null vs = return (xs, id)
      | Just (v,vdef,ys) <- findEq vs xs = do
          let mt1 m = IM.insert v (evalRat m vdef) m
          (zs, mt2) <- f (IS.delete v vs) [subst1Constr v (fromRat vdef) c | c <- ys]
          return (zs, mt1 . mt2)
      | otherwise =
          case IS.minView vs of
            Nothing -> return (xs, id) -- should not happen
            Just (v,vs') ->
              case collectBounds v xs of
                (bnd, rest) -> do
                  cond <- boundsToConstrs bnd
                  let mt1 m =
                       case Interval.simplestRationalWithin (evalBounds m bnd) of
                         Nothing  -> error "ToySolver.FourierMotzkin.project': should not happen"
                         Just val -> IM.insert v val m
                  (ys, mt2) <- f vs' (rest ++ cond)
                  return (ys, mt1 . mt2)

findEq :: VarSet -> [Constr] -> Maybe (Var, Rat, [Constr])
findEq vs = msum . map f . pickup
  where
    pickup :: [a] -> [(a,[a])]
    pickup [] = []
    pickup (x:xs) = (x,xs) : [(y,x:ys) | (y,ys) <- pickup xs]

    f :: (Constr, [Constr]) -> Maybe (Var, Rat, [Constr])
    f (IsZero e, cs) = do
      let vs2 = IS.intersection vs (vars e)
      guard $ not $ IS.null vs2
      let v = IS.findMin vs2
          (c, e') = LA.extract v e
      return (v, (negateV e', c), cs)
    f _ = Nothing

-- | @'solve' {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LRA φ@ when
-- such @M@ exists, returns @Nothing@ otherwise.
--
-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
--
solve :: VarSet -> [LA.Atom Rational] -> Maybe (Model Rational)
solve vs cs = msum [solve' vs cs2 | cs2 <- unDNF (constraintsToDNF cs)]

solve' :: VarSet -> [Constr] -> Maybe (Model Rational)
solve' vs cs = do
  (cs2,mt) <- projectN' vs =<< simplify cs
  let m :: Model Rational
      m = IM.empty
  guard $ all (eval m) cs2
  return $ mt m

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

gcd' :: [Integer] -> Integer
gcd' [] = 1
gcd' xs = foldl1' gcd xs

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

1	{-# LANGUAGE FlexibleInstances #-}
2	{-# LANGUAGE MultiParamTypeClasses #-}
3	{-# LANGUAGE TypeSynonymInstances #-}
4	{-# OPTIONS_GHC -Wall #-}
5	{-# OPTIONS_HADDOCK show-extensions #-}
6	-----------------------------------------------------------------------------
7	-- \|
8	-- Module : ToySolver.Arith.FourierMotzkin.Base
9	-- Copyright : (c) Masahiro Sakai 2011-2013
10	-- License : BSD-style
11	--
12	-- Maintainer : masahiro.sakai@gmail.com
13	-- Stability : provisional
14	-- Portability : non-portable
15	--
16	-- Naïve implementation of Fourier-Motzkin Variable Elimination
17	--
18	-- Reference:
19	--
20	-- * <http://users.cecs.anu.edu.au/~michaeln/pubs/arithmetic-dps.pdf>
21	--
22	-----------------------------------------------------------------------------
23	module ToySolver.Arith.FourierMotzkin.Base
24	(
25	-- * Primitive constraints
26	Constr (..)
27	, eqR
28	, ExprZ
29	, fromLAAtom
30	, toLAAtom
31	, constraintsToDNF
32	, simplify
33
34	-- * Bounds
35	, Bounds
36	, evalBounds
37	, boundsToConstrs
38	, collectBounds
39
40	-- * Projection
41	, project
42	, project'
43	, projectN
44	, projectN'
45
46	-- * Solving
47	, solve
48	, solve'
49
50	-- * Utilities used by other modules
51	, Rat
52	, toRat
53	) where
54
55	import Control.Monad
56	import Data.List (foldl', foldl1')
57	import Data.Maybe
58	import Data.Ratio
59	import qualified Data.IntMap as IM
60	import qualified Data.IntSet as IS
61	import Data.VectorSpace hiding (project)
62	import qualified Data.Interval as Interval
63	import Data.Interval (Interval, Extended (..), (<=..<), (<..<=), (<..<))
64
65	import ToySolver.Data.OrdRel
66	import ToySolver.Data.Boolean
67	import ToySolver.Data.DNF
68	import qualified ToySolver.Data.LA as LA
69	import ToySolver.Data.IntVar
70
71	-- ---------------------------------------------------------------------------
72
73	type ExprZ = LA.Expr Integer
74
75	-- \| (t,c) represents t/c, and c must be >0.
76	type Rat = (ExprZ, Integer)
77
78	toRat :: LA.Expr Rational -> Rat
79	toRat e = seq m $ (LA.mapCoeff f e, m)	1✔
80	where
81	f x = numerator (fromInteger m * x)	1✔
82	m = foldl' lcm 1 [denominator c \| (c,_) <- LA.terms e]	1✔
83
84	fromRat :: Rat -> LA.Expr Rational
85	fromRat (e,c) = LA.mapCoeff (% c) e	1✔
86
87	evalRat :: Model Rational -> Rat -> Rational
88	evalRat model (e, d) = LA.lift1 1 (model IM.!) (LA.mapCoeff fromIntegral e) / (fromIntegral d)	1✔
89
90	-- ---------------------------------------------------------------------------
91
92	-- \| Atomic constraint
93	data Constr
94	= IsNonneg ExprZ
95	-- ^ e ≥ 0
96	\| IsPos ExprZ
97	-- ^ e > 0
98	\| IsZero ExprZ
99	-- ^ e = 0
100	deriving (Show, Eq, Ord)	×
101
102	instance Variables Constr where
103	vars (IsPos t) = vars t	×
104	vars (IsNonneg t) = vars t	×
105	vars (IsZero t) = vars t	×
106
107	leR, ltR, geR, gtR, eqR :: Rat -> Rat -> Constr
108	leR (e1,c) (e2,d) = IsNonneg $ normalizeExpr $ c ^ e2 ^-^ d ^ e1	1✔
109	ltR (e1,c) (e2,d) = IsPos $ normalizeExpr $ c ^ e2 ^-^ d ^ e1	1✔
110	geR = flip leR	1✔
111	gtR = flip ltR	1✔
112	eqR (e1,c) (e2,d) = IsZero $ normalizeExpr $ c ^ e2 ^-^ d ^ e1	1✔
113
114	normalizeExpr :: ExprZ -> ExprZ
115	normalizeExpr e = LA.mapCoeff (`div` d) e	1✔
116	where d = abs $ gcd' $ map fst $ LA.terms e	1✔
117
118	-- "subst1Constr x t c" computes c[t/x]
119	subst1Constr :: Var -> LA.Expr Rational -> Constr -> Constr
120	subst1Constr x t c =	1✔
121	case c of	1✔
122	IsPos e -> IsPos (f e)	1✔
123	IsNonneg e -> IsNonneg (f e)	1✔
124	IsZero e -> IsZero (f e)	1✔
125	where
126	f :: ExprZ -> ExprZ
127	f = normalizeExpr . fst . toRat . LA.applySubst1 x t . LA.mapCoeff fromInteger	1✔
128
129	-- \| Simplify conjunction of 'Constr's.
130	-- It returns 'Nothing' when a inconsistency is detected.
131	simplify :: [Constr] -> Maybe [Constr]
132	simplify = fmap concat . mapM f	1✔
133	where
134	f :: Constr -> Maybe [Constr]
135	f c@(IsPos e) =	1✔
136	case LA.asConst e of	1✔
137	Just x -> guard (x > 0) >> return []	1✔
138	Nothing -> return [c]	1✔
139	f c@(IsNonneg e) =
140	case LA.asConst e of	1✔
141	Just x -> guard (x >= 0) >> return []	1✔
142	Nothing -> return [c]	1✔
143	f c@(IsZero e) =
144	case LA.asConst e of	1✔
145	Just x -> guard (x == 0) >> return []	1✔
146	Nothing -> return [c]	1✔
147
148	instance Eval (Model Rational) Constr Bool where
149	eval m (IsPos t) = LA.eval m (LA.mapCoeff fromInteger t :: LA.Expr Rational) > 0	1✔
150	eval m (IsNonneg t) = LA.eval m (LA.mapCoeff fromInteger t :: LA.Expr Rational) >= 0	1✔
151	eval m (IsZero t) = LA.eval m (LA.mapCoeff fromInteger t :: LA.Expr Rational) == 0	×
152
153	-- ---------------------------------------------------------------------------
154
155	fromLAAtom :: LA.Atom Rational -> DNF Constr
156	fromLAAtom (OrdRel a op b) = atomR' op (toRat a) (toRat b)	1✔
157	where
158	atomR' :: RelOp -> Rat -> Rat -> DNF Constr
159	atomR' op a b =	1✔
160	case op of	1✔
161	Le -> DNF [[a `leR` b]]	1✔
162	Lt -> DNF [[a `ltR` b]]	1✔
163	Ge -> DNF [[a `geR` b]]	1✔
164	Gt -> DNF [[a `gtR` b]]	1✔
165	Eql -> DNF [[a `eqR` b]]	1✔
166	NEq -> DNF [[a `ltR` b], [a `gtR` b]]	×
167
168	toLAAtom :: Constr -> LA.Atom Rational
169	toLAAtom (IsNonneg e) = LA.mapCoeff fromInteger e .>=. LA.constant 0	×
170	toLAAtom (IsPos e) = LA.mapCoeff fromInteger e .>. LA.constant 0	×
171	toLAAtom (IsZero e) = LA.mapCoeff fromInteger e .==. LA.constant 0	×
172
173	constraintsToDNF :: [LA.Atom Rational] -> DNF Constr
174	constraintsToDNF = andB . map fromLAAtom	1✔
175
176	-- ---------------------------------------------------------------------------
177
178	{-
179	(ls1,ls2,us1,us2) represents
180	{ x \| ∀(M,c)∈ls1. M/c≤x, ∀(M,c)∈ls2. M/c<x, ∀(M,c)∈us1. x≤M/c, ∀(M,c)∈us2. x<M/c }
181	-}
182	type Bounds = ([Rat], [Rat], [Rat], [Rat])
183
184	evalBounds :: Model Rational -> Bounds -> Interval Rational
185	evalBounds model (ls1,ls2,us1,us2) =	1✔
186	Interval.intersections $	1✔
187	[ Finite (evalRat model x) <=..< PosInf \| x <- ls1 ] ++	1✔
188	[ Finite (evalRat model x) <..< PosInf \| x <- ls2 ] ++	1✔
189	[ NegInf <..<= Finite (evalRat model x) \| x <- us1 ] ++	1✔
190	[ NegInf <..< Finite (evalRat model x) \| x <- us2 ]	1✔
191
192	boundsToConstrs :: Bounds -> Maybe [Constr]
193	boundsToConstrs (ls1, ls2, us1, us2) = simplify $	1✔
194	[ x `leR` y \| x <- ls1, y <- us1 ] ++	1✔
195	[ x `ltR` y \| x <- ls1, y <- us2 ] ++	1✔
196	[ x `ltR` y \| x <- ls2, y <- us1 ] ++	1✔
197	[ x `ltR` y \| x <- ls2, y <- us2 ]	1✔
198
199	collectBounds :: Var -> [Constr] -> (Bounds, [Constr])
200	collectBounds v = foldr phi (([],[],[],[]),[])	1✔
201	where
202	phi :: Constr -> (Bounds, [Constr]) -> (Bounds, [Constr])
203	phi constr@(IsNonneg t) x = f False constr t x	1✔
204	phi constr@(IsPos t) x = f True constr t x	1✔
205	phi constr@(IsZero t) (bnd@(ls1,ls2,us1,us2), xs) =
206	case LA.extractMaybe v t of	1✔
207	Nothing -> (bnd, constr : xs)	1✔
208	Just (c,t') -> ((t'' : ls1, ls2, t'' : us1, us2), xs)	×
209	where
210	t'' = (signum c *^ negateV t', abs c)	×
211
212	f :: Bool -> Constr -> ExprZ -> (Bounds, [Constr]) -> (Bounds, [Constr])
213	f strict constr t (bnd@(ls1,ls2,us1,us2), xs) =	1✔
214	case LA.extract v t of	1✔
215	(c,t') ->
216	case c `compare` 0 of	1✔
217	EQ -> (bnd, constr : xs)	1✔
218	GT ->
219	if strict	1✔
220	then ((ls1, (negateV t', c) : ls2, us1, us2), xs) -- 0 < cx + M ⇔ -M/c < x	1✔
221	else (((negateV t', c) : ls1, ls2, us1, us2), xs) -- 0 ≤ cx + M ⇔ -M/c ≤ x	1✔
222	LT ->
223	if strict	1✔
224	then ((ls1, ls2, us1, (t', negate c) : us2), xs) -- 0 < cx + M ⇔ x < M/-c	1✔
225	else ((ls1, ls2, (t', negate c) : us1, us2), xs) -- 0 ≤ cx + M ⇔ x ≤ M/-c	1✔
226
227	-- ---------------------------------------------------------------------------
228
229	{-\| @'project' x φ@ returns @[(ψ_1, lift_1), …, (ψ_m, lift_m)]@ such that:
230
231	* @⊢_LRA ∀y1, …, yn. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_m)@ where @{y1, …, yn} = FV(φ) \\ {x}@, and
232
233	* if @M ⊧_LRA ψ_i@ then @lift_i M ⊧_LRA φ@.
234	-}
235	project :: Var -> [LA.Atom Rational] -> [([LA.Atom Rational], Model Rational -> Model Rational)]
236	project v xs = do	×
237	ys <- unDNF $ constraintsToDNF xs	×
238	(zs, mt) <- maybeToList $ project' v ys	×
239	return (map toLAAtom zs, mt)	×
240
241	project' :: Var -> [Constr] -> Maybe ([Constr], Model Rational -> Model Rational)
242	project' v cs = projectN' (IS.singleton v) cs	×
243
244	{-\| @'projectN' {x1,…,xm} φ@ returns @[(ψ_1, lift_1), …, (ψ_n, lift_n)]@ such that:
245
246	* @⊢_LRA ∀y1, …, yp. (∃x. φ) ↔ (ψ_1 ∨ … ∨ φ_n)@ where @{y1, …, yp} = FV(φ) \\ {x1,…,xm}@, and
247
248	* if @M ⊧_LRA ψ_i@ then @lift_i M ⊧_LRA φ@.
249	-}
250	projectN :: VarSet -> [LA.Atom Rational] -> [([LA.Atom Rational], Model Rational -> Model Rational)]
251	projectN vs xs = do	×
252	ys <- unDNF $ constraintsToDNF xs	×
253	(zs, mt) <- maybeToList $ projectN' vs ys	×
254	return (map toLAAtom zs, mt)	×
255
256	projectN' :: VarSet -> [Constr] -> Maybe ([Constr], Model Rational -> Model Rational)
257	projectN' = f	1✔
258	where
259	f vs xs	1✔
260	\| IS.null vs = return (xs, id)	1✔
261	\| Just (v,vdef,ys) <- findEq vs xs = do	1✔
262	let mt1 m = IM.insert v (evalRat m vdef) m	1✔
263	(zs, mt2) <- f (IS.delete v vs) [subst1Constr v (fromRat vdef) c \| c <- ys]	1✔
264	return (zs, mt1 . mt2)	1✔
265	\| otherwise =	×
266	case IS.minView vs of	1✔
267	Nothing -> return (xs, id) -- should not happen	×
268	Just (v,vs') ->
269	case collectBounds v xs of	1✔
270	(bnd, rest) -> do	1✔
271	cond <- boundsToConstrs bnd	1✔
272	let mt1 m =	1✔
273	case Interval.simplestRationalWithin (evalBounds m bnd) of	1✔
274	Nothing -> error "ToySolver.FourierMotzkin.project': should not happen"	×
275	Just val -> IM.insert v val m	1✔
276	(ys, mt2) <- f vs' (rest ++ cond)	1✔
277	return (ys, mt1 . mt2)	1✔
278
279	findEq :: VarSet -> [Constr] -> Maybe (Var, Rat, [Constr])
280	findEq vs = msum . map f . pickup	1✔
281	where
282	pickup :: [a] -> [(a,[a])]
283	pickup [] = []	1✔
284	pickup (x:xs) = (x,xs) : [(y,x:ys) \| (y,ys) <- pickup xs]	1✔
285
286	f :: (Constr, [Constr]) -> Maybe (Var, Rat, [Constr])
287	f (IsZero e, cs) = do	1✔
288	let vs2 = IS.intersection vs (vars e)	1✔
289	guard $ not $ IS.null vs2	1✔
290	let v = IS.findMin vs2	1✔
291	(c, e') = LA.extract v e	1✔
292	return (v, (negateV e', c), cs)	1✔
293	f _ = Nothing	1✔
294
295	-- \| @'solve' {x1,…,xn} φ@ returns @Just M@ that @M ⊧_LRA φ@ when
296	-- such @M@ exists, returns @Nothing@ otherwise.
297	--
298	-- @FV(φ)@ must be a subset of @{x1,…,xn}@.
299	--
300	solve :: VarSet -> [LA.Atom Rational] -> Maybe (Model Rational)
301	solve vs cs = msum [solve' vs cs2 \| cs2 <- unDNF (constraintsToDNF cs)]	1✔
302
303	solve' :: VarSet -> [Constr] -> Maybe (Model Rational)
304	solve' vs cs = do	1✔
305	(cs2,mt) <- projectN' vs =<< simplify cs	1✔
306	let m :: Model Rational
307	m = IM.empty	1✔
308	guard $ all (eval m) cs2	1✔
309	return $ mt m	1✔
310
311	-- ---------------------------------------------------------------------------
312
313	gcd' :: [Integer] -> Integer
314	gcd' [] = 1	×
315	gcd' xs = foldl1' gcd xs	1✔
316
317	-- ---------------------------------------------------------------------------

msakai / toysolver / 767

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous