# File ‹~~/src/Tools/Argo/argo_simplex.ML›

```(*  Title:      Tools/Argo/argo_simplex.ML
Author:     Sascha Boehme

Linear arithmetic reasoning based on the simplex algorithm. It features:

* simplification and normalization of arithmetic expressions
* decision procedure for reals

* propagating implied inequality literals while assuming external knowledge
* propagating equalities for fixed variables to all other theory solvers
* pruning the tableau after new atoms have been added: eliminate unnecessary
variables

The implementation is inspired by:

Bruno Dutertre and Leonardo de Moura. A fast linear-arithmetic solver
for DPLL(T). In Computer Aided Verification, pages 81-94. Springer, 2006.
*)

signature ARGO_SIMPLEX =
sig
(* context *)
type context
val context: context

(* enriching the context *)
val add_atom: Argo_Term.term -> context -> Argo_Lit.literal option * context

(* main operations *)
val prepare: context -> context
val assume: Argo_Common.literal -> context -> Argo_Lit.literal Argo_Common.implied * context
val check: context -> Argo_Lit.literal Argo_Common.implied * context
val explain: Argo_Lit.literal -> context -> (Argo_Cls.clause * context) option
val backtrack: context -> context
end

structure Argo_Simplex: ARGO_SIMPLEX =
struct

(* extended rationals *)

(*
Extended rationals (c, k) are reals (c + k * e) where e is some small positive real number.
Extended rationals are used to represent a strict inequality by a non-strict inequality:
c < x  ~~  c + k * e <= e
x < c  ~~  x <= c - k * e
*)

type erat = Rat.rat * Rat.rat

val erat_zero = (@0, @0)

fun add (c1, k1) (c2, k2) = (c1 + c2, k1 + k2)
fun sub (c1, k1) (c2, k2) = (c1 - c2, k1 - k2)
fun mul n (c, k) = (n * c, n * k)

val erat_ord = prod_ord Rat.ord Rat.ord

fun less_eq n1 n2 = is_less_equal (erat_ord (n1, n2))
fun less n1 n2 = is_less (erat_ord (n1, n2))

(* term functions *)

fun dest_monom (Argo_Term.T (_, Argo_Expr.Mul, [Argo_Term.T (_, Argo_Expr.Num n, _), t])) = (t, n)
| dest_monom t = (t, @1)

datatype node = Var of Argo_Term.term | Num of Rat.rat
datatype ineq = Lower of Argo_Term.term * erat | Upper of Argo_Term.term * erat

fun dest_node (Argo_Term.T (_, Argo_Expr.Num n, _)) = Num n
| dest_node t = Var t

fun dest_atom true (k as Argo_Expr.Le) t1 t2 = SOME (k, dest_node t1, dest_node t2)
| dest_atom true (k as Argo_Expr.Lt) t1 t2 = SOME (k, dest_node t1, dest_node t2)
| dest_atom false Argo_Expr.Le t1 t2 = SOME (Argo_Expr.Lt, dest_node t2, dest_node t1)
| dest_atom false Argo_Expr.Lt t1 t2 = SOME (Argo_Expr.Le, dest_node t2, dest_node t1)
| dest_atom _ _ _ _ = NONE

fun ineq_of pol (Argo_Term.T (_, k, [t1, t2])) =
(case dest_atom pol k t1 t2 of
SOME (Argo_Expr.Le, Var x, Num n) => SOME (Upper (x, (n, @0)))
| SOME (Argo_Expr.Le, Num n, Var x) => SOME (Lower (x, (n, @0)))
| SOME (Argo_Expr.Lt, Var x, Num n) => SOME (Upper (x, (n, @~1)))
| SOME (Argo_Expr.Lt, Num n, Var x) => SOME (Lower (x, (n, @1)))
| _ => NONE)
| ineq_of _ _ = NONE

(* proofs *)

(*
comment missing
*)

fun mk_ineq is_lt = if is_lt then Argo_Expr.mk_lt else Argo_Expr.mk_le
fun ineq_rule_of is_lt = if is_lt then Argo_Proof.Lt else Argo_Proof.Le

fun rewrite_top f = Argo_Rewr.rewrite_top (f Argo_Rewr.context)

fun unnegate_conv (e as Argo_Expr.E (Argo_Expr.Not, [Argo_Expr.E (Argo_Expr.Le, [e1, e2])])) =
Argo_Rewr.rewr (Argo_Proof.Rewr_Not_Ineq Argo_Proof.Le) (Argo_Expr.mk_lt e2 e1) e
| unnegate_conv (e as Argo_Expr.E (Argo_Expr.Not, [Argo_Expr.E (Argo_Expr.Lt, [e1, e2])])) =
Argo_Rewr.rewr (Argo_Proof.Rewr_Not_Ineq Argo_Proof.Lt) (Argo_Expr.mk_le e2 e1) e
| unnegate_conv e = Argo_Rewr.keep e

fun scale_conv r mk n e1 e2 =
let
fun scale e = Argo_Expr.mk_mul (Argo_Expr.mk_num n) e
val (e1, e2) = if n > @0 then (scale e1, scale e2) else (scale e2, scale e1)
val conv = Argo_Rewr.rewr (Argo_Proof.Rewr_Ineq_Mul (r, n)) (mk e1 e2)
in Argo_Rewr.seq [conv, Argo_Rewr.args (rewrite_top Argo_Rewr.norm_mul)] end

fun dest_ineq (Argo_Expr.E (Argo_Expr.Le, [e1, e2])) = SOME (false, e1, e2)
| dest_ineq (Argo_Expr.E (Argo_Expr.Lt, [e1, e2])) = SOME (true, e1, e2)
| dest_ineq _ = NONE

fun scale_ineq_conv n e =
if n = @1 then Argo_Rewr.keep e
else
(case dest_ineq e of
NONE => raise Fail "bad inequality"
| SOME (is_lt, e1, e2) => scale_conv (ineq_rule_of is_lt) (mk_ineq is_lt) n e1 e2 e)

fun simp_lit (n, (lit, p)) =
let val conv = Argo_Rewr.seq [unnegate_conv, scale_ineq_conv n]
in Argo_Rewr.with_proof conv (Argo_Lit.signed_expr_of lit, p) end

fun reduce_conv r = Argo_Rewr.rewr (Argo_Proof.Rewr_Ineq_Nums (r, false)) Argo_Expr.false_expr

fun simp_combine es p prf =
let
fun dest e (is_lt, (es1, es2)) =
let val (is_lt', e1, e2) = the (dest_ineq e)
in (is_lt orelse is_lt', (e1 :: es1, e2 :: es2)) end
val (is_lt, (es1, es2)) = fold_rev dest es (false, ([], []))
val e = uncurry (mk_ineq is_lt) (apply2 Argo_Expr.mk_add (es1, es2))
val conv = Argo_Rewr.seq [Argo_Rewr.args combine_conv, reduce_conv (ineq_rule_of is_lt)]
in prf |> Argo_Rewr.with_proof conv (e, p) |>> snd end

fun linear_combination nlps prf =
let val ((es, ps), prf) = fold_map simp_lit nlps prf |>> split_list
in prf |> Argo_Proof.mk_linear_comb ps |-> simp_combine es |-> Argo_Proof.mk_lemma [] end

fun proof_of (lit, SOME p) (ls, prf) = ((lit, p), (ls, prf))
| proof_of (lit, NONE) (ls, prf) =
let val (p, prf) = Argo_Proof.mk_hyp lit prf
in ((lit, p), (Argo_Lit.negate lit :: ls, prf)) end

(* tableau *)

(*
The tableau consists of equations x_i = a_i1 * x_1 + ... a_ik * x_k where
the variable on the left-hand side is called a basic variable and
the variables on the right-hand side are called non-basic variables.

For each basic variable, the polynom on the right-hand side is stored as a map
from variables to coefficients. Only variables with non-zero coefficients are stored.
The map is sorted by the term order of the variables for a deterministic order when
analyzing a polynom.

Additionally, for each basic variable a boolean flag is kept that, when false,
indicates that the current value of the basic variable might be outside its bounds.
The value of a non-basic variable is always within its bounds.

The tableau is stored as a table indexed by variables. For each variable,
both basic and non-basic, its current value is stored as extended rational
along with either the equations or the occurrences.
*)

type basic = bool * (Argo_Term.term * Rat.rat) Ord_List.T
type entry = erat * basic option
type tableau = entry Argo_Termtab.table

fun dirty ms = SOME (false, ms)
fun checked ms = SOME (true, ms)

fun basic_entry ms = (erat_zero, dirty ms)
val non_basic_entry: entry = (erat_zero, NONE)

fun value_of tableau x =
(case Argo_Termtab.lookup tableau x of
NONE => erat_zero
| SOME (v, _) => v)

fun first_unchecked_basic tableau =
Argo_Termtab.get_first (fn (y, (v, SOME (false, ms))) => SOME (y, v, ms) | _ => NONE) tableau

local

fun coeff_of ms x = the (AList.lookup Argo_Term.eq_term ms x)

val eq_var = Argo_Term.eq_term
fun monom_ord sp = prod_ord Argo_Term.term_ord (K EQUAL) sp

fun add_monom m ms = Ord_List.insert monom_ord m ms
fun update_monom (m as (x, a)) = if a = @0 then AList.delete eq_var x else AList.update eq_var m

fun add_scaled_monom n (x, a) ms =
(case AList.lookup eq_var ms x of
NONE => add_monom (x, n * a) ms
| SOME b => update_monom (x, n * a + b) ms)

fun replace_polynom x n ms' ms = fold (add_scaled_monom n) ms' (AList.delete eq_var x ms)

fun map_basic f (v, SOME (_, ms)) = f v ms
| map_basic _ e = e

fun map_basic_entries x f =
let
fun apply (e as (v, SOME (_, ms))) = if AList.defined eq_var ms x then f v ms else e
| apply ve = ve
in Argo_Termtab.map (K apply) end

fun put_entry x e = Argo_Termtab.update (x, e)

let val ms = Ord_List.make monom_ord (map dest_monom ts)
in fold (fn (x, _) => put_entry x non_basic_entry) ms (put_entry y (basic_entry ms) tableau) end
| add_new_entry x tableau = put_entry x non_basic_entry tableau

fun with_non_basic update_basic x f tableau =
(case Argo_Termtab.lookup tableau x of
NONE => tableau
| SOME (v, NONE) => f v tableau
| SOME (v, SOME (_, ms)) => if update_basic then put_entry x (v, dirty ms) tableau else tableau)

in

if Argo_Termtab.defined tableau x then tableau

fun basic_within_bounds y = Argo_Termtab.map_entry y (map_basic (fn v => fn ms => (v, checked ms)))

fun eliminate _ tableau = tableau

fun update_non_basic pred x v' = with_non_basic true x (fn v =>
let fun update_basic n v ms = (add v (mul (coeff_of ms x) n), dirty ms)
in pred v ? put_entry x (v', NONE) o map_basic_entries x (update_basic (sub v' v)) end)

fun update_pivot y vy ms x c v = with_non_basic false x (fn vx =>
let
val a = Rat.inv c
val v' = mul a (sub v vy)

fun scale_or_drop (x', b) = if Argo_Term.eq_term (x', x) then NONE else SOME (x', ~ a * b)
val ms = add_monom (y, a) (map_filter scale_or_drop ms)

fun update_basic v ms' =
let val n = coeff_of ms' x
in (add v (mul n v'), dirty (replace_polynom x n ms ms')) end
in
put_entry x (add vx v', dirty ms) #>
put_entry y (v, NONE) #>
map_basic_entries x update_basic
end)

end

(* bounds *)

(*
comment missing
*)

type bound = (erat * Argo_Common.literal) option
type atoms = (erat * Argo_Term.term) list
type bounds_atoms = ((bound * bound) * (atoms * atoms))
type bounds = bounds_atoms Argo_Termtab.table

val empty_bounds_atoms: bounds_atoms = ((NONE, NONE), ([], []))

fun on_some pred (SOME (n, _)) = pred n
| on_some _ NONE = false

fun none_or_some pred (SOME (n, _)) = pred n
| none_or_some _ NONE = true

fun bound_of (SOME (n, _)) = n
| bound_of NONE = raise Fail "bad bound"

fun reason_of (SOME (_, r)) = r
| reason_of NONE = raise Fail "bad reason"

fun bounds_atoms_of bounds x = the_default empty_bounds_atoms (Argo_Termtab.lookup bounds x)
fun bounds_of bounds x = fst (bounds_atoms_of bounds x)

fun put_bounds x bs bounds = Argo_Termtab.map_default (x, empty_bounds_atoms) (apfst (K bs)) bounds

fun has_bound_atoms bounds x =
(case Argo_Termtab.lookup bounds x of
NONE => false
| SOME (_, ([], [])) => false
| _ => true)

fun add_new_atom f x n t =
let val ins = f (insert (eq_snd Argo_Term.eq_term) (n, t))
in Argo_Termtab.map_default (x, empty_bounds_atoms) (apsnd ins) end

fun del_atom x t =
let fun eq_atom (t1, (_, t2)) = Argo_Term.eq_term (t1, t2)
in Argo_Termtab.map_entry x (apsnd (apply2 (remove eq_atom t))) end

(* context *)

type context = {
tableau: tableau, (* values of variables and tableau entries for each variable *)
bounds: bounds, (* bounds and unassigned atoms for each variable *)
prf: Argo_Proof.context, (* proof context *)
back: bounds list} (* stack storing previous bounds and unassigned atoms *)

fun mk_context tableau bounds prf back: context =
{tableau=tableau, bounds=bounds, prf=prf, back=back}

val context = mk_context Argo_Termtab.empty Argo_Termtab.empty Argo_Proof.simplex_context []

(* declaring atoms *)

fun add_ineq_atom f t x n ({tableau, bounds, prf, back}: context) =
(* TODO: check whether the atom is already known to hold *)
(NONE, mk_context (add_entry x tableau) (add_new_atom f x n t bounds) prf back)

(case ineq_of true t of
SOME (Lower (x, n)) => add_ineq_atom apfst t x n cx
| SOME (Upper (x, n)) => add_ineq_atom apsnd t x n cx
| NONE => (NONE, cx))

(* preparing the solver after new atoms have been added *)

(*
Variables that do not directly occur in atoms can be eliminated from the tableau
since no bounds will ever limit their value. This can reduce the tableau size
substantially.
*)

fun prepare ({tableau, bounds, prf, back}: context) =
let fun drop (xe as (x, _)) = not (has_bound_atoms bounds x) ? eliminate xe
in mk_context (Argo_Termtab.fold drop tableau tableau) bounds prf back end

(* assuming external knowledge *)

fun bounds_conflict r1 r2 ({tableau, bounds, prf, back}: context) =
let
val ((lp2, lp1), (lits, prf)) = ([], prf) |> proof_of r2 ||>> proof_of r1
val (p, prf) = linear_combination [(@~1, lp1), (@1, lp2)] prf
in (Argo_Common.Conflict (lits, p), mk_context tableau bounds prf back) end

fun assume_bounds order x c bs ({tableau, bounds, prf, back}: context) =
let
val lits = []
val bounds = put_bounds x bs bounds
val tableau = update_non_basic (fn v => erat_ord (v, c) = order) x c tableau
in (Argo_Common.Implied lits, mk_context tableau bounds prf back) end

fun assume_lower r x c (low, upp) cx =
if on_some (fn l => less_eq c l) low then (Argo_Common.Implied [], cx)
else if on_some (fn u => less u c) upp then bounds_conflict r (reason_of upp) cx
else assume_bounds LESS x c (SOME (c, r), upp) cx

fun assume_upper r x c (low, upp) cx =
if on_some (fn u => less_eq u c) upp then (Argo_Common.Implied [], cx)
else if on_some (fn l => less c l) low then bounds_conflict (reason_of low) r cx
else assume_bounds GREATER x c (low, SOME (c, r)) cx

fun with_bounds r t f x n ({tableau, bounds, prf, back}: context) =
f r x n (bounds_of bounds x) (mk_context tableau (del_atom x t bounds) prf back)

fun choose f (SOME (Lower (x, n))) cx = f assume_lower x n cx
| choose f (SOME (Upper (x, n))) cx = f assume_upper x n cx
| choose _ NONE cx = (Argo_Common.Implied [], cx)

fun assume (r as (lit, _)) cx =
let val (t, pol) = Argo_Lit.dest lit
in choose (with_bounds r t) (ineq_of pol t) cx end

(* checking for consistency and pending implications *)

fun basic_bounds_conflict lower y ms ({tableau, bounds, prf, back}: context) =
let
val (a, low, upp) = if lower then (@1, fst, snd) else (@~1, snd, fst)
fun coeff_proof f a x = apfst (pair a) o proof_of (reason_of (f (bounds_of bounds x)))
fun monom_proof (x, a) = coeff_proof (if a < @0 then low else upp) a x
val ((alp, alps), (lits, prf)) = ([], prf) |> coeff_proof low a y ||>> fold_map monom_proof ms
val (p, prf) = linear_combination (alp :: alps) prf
in (Argo_Common.Conflict (lits, p), mk_context tableau bounds prf back) end

fun can_compensate ord tableau bounds (x, a) =
let val (low, upp) = bounds_of bounds x
in
if Rat.ord (a, @0) = ord then none_or_some (fn u => less (value_of tableau x) u) upp
else none_or_some (fn l => less l (value_of tableau x)) low
end

fun check (cx as {tableau, bounds, prf, back}: context) =
(case first_unchecked_basic tableau of
NONE => (Argo_Common.Implied [], cx)
| SOME (y, v, ms) =>
let val (low, upp) = bounds_of bounds y
in
if on_some (fn l => less v l) low then adjust GREATER true y v ms (bound_of low) cx
else if on_some (fn u => less u v) upp then adjust LESS false y v ms (bound_of upp) cx
else check (mk_context (basic_within_bounds y tableau) bounds prf back)
end)

and adjust ord lower y vy ms v (cx as {tableau, bounds, prf, back}: context) =
(case find_first (can_compensate ord tableau bounds) ms of
NONE => basic_bounds_conflict lower y ms cx
| SOME (x, a) => check (mk_context (update_pivot y vy ms x a v tableau) bounds prf back))

(* explanations *)

fun explain _ _ = NONE

(* backtracking *)

fun add_level ({tableau, bounds, prf, back}: context) =
mk_context tableau bounds prf (bounds :: back)

fun backtrack ({back=[], ...}: context) = raise Empty
| backtrack ({tableau, prf, back=bounds :: back, ...}: context) =
mk_context tableau bounds prf back

end
```