From 6f1634d2f822037a482436a64d3ef3bfb2fac2a0 Mon Sep 17 00:00:00 2001 From: Frédéric Besson Date: Fri, 8 Mar 2019 09:56:49 +0100 Subject: Several improvements and fixes of Lia - Improved reification for Micromega (support for #8764) - Fixes #9268: Do not take universes into account in lia reification Improve #9291 by threading the evar_map during reification. Universes are unified. - Remove (potentially cyclic) dependency over lra for Rle_abs - Towards a complete simplex-based lia fixes #9615 Lia is now exclusively using cutting plane proofs. For this to always work, all the variables need to be positive. Therefore, lia is pre-processing the goal for each variable x it introduces the constraints x = y - z , y>=0 , z >= 0 for some fresh variable y and z. For scalability, nia is currently NOT performing this pre-processing. - Lia is using the FSet library manual merge of commit #230899e87c51c12b2f21b6fedc414d099a1425e4 to work around a "leaked" hint breaking compatibility of eauto --- plugins/micromega/simplex.ml | 216 +++++++++++++++++++++++++++---------------- 1 file changed, 134 insertions(+), 82 deletions(-) (limited to 'plugins/micromega/simplex.ml') diff --git a/plugins/micromega/simplex.ml b/plugins/micromega/simplex.ml index 4465aa1ee1..4ddeb6c2c0 100644 --- a/plugins/micromega/simplex.ml +++ b/plugins/micromega/simplex.ml @@ -11,9 +11,11 @@ (** A naive simplex *) open Polynomial open Num -open Util +(*open Util*) open Mutils +type ('a,'b) sum = Inl of 'a | Inr of 'b + let debug = false type iset = unit IMap.t @@ -130,12 +132,6 @@ let is_maximised rst v = violating a restriction. *) -(* let is_unbounded rst tbl vr = - IMap.for_all (fun x v -> if Vect.get vr v =/ Int 0 - then None else Some(x,v) in - let (c,r) = extract is_conflict l in - match c with - | Some (c,_) -> Some (c,r) - | None -> match l with - | [] -> None - | e::l -> Some(e,l) - -(*let remove_redundant rst t = - IMap.fold (fun k v m -> if Restricted.is_restricted k rst && Vect.for_all (fun x _ -> x == 0 || Restricted.is_restricted x rst) v - then begin - if debug then - Printf.printf "%a is redundant\n" LinPoly.pp_var k; - IMap.remove k m - end - else m) t t - *) + + let rec most_violating l e (x,v) rst = + match l with + | [] -> Some((x,v),rst) + | (x',v')::l -> + let e' = Vect.dotproduct esol v' in + if e' <=/ e + then most_violating l e' (x',v') ((x,v)::rst) + else most_violating l e (x,v) ((x',v')::rst) in + + match l with + | [] -> None + | (x,v)::l -> let e = Vect.dotproduct esol v in + most_violating l e (x,v) [] + let rec solve opt l (rst:Restricted.t) (t:tableau) = @@ -515,65 +509,117 @@ let make_farkas_proof (env: WithProof.t IMap.t) vm v = WithProof.mult (Vect.cst n) (IMap.find x env) end) WithProof.zero v -(* -let incr_cut rmin x = - match rmin with - | None -> true - | Some r -> Int.compare x r = 1 - *) -let cut env rmin sol vm (rst:Restricted.t) (x,v) = -(* if not (incr_cut rmin x) - then None - else *) - let (n,r) = Vect.decomp_cst v in +let frac_num n = n -/ Num.floor_num n - let nf = Num.floor_num n in - if nf =/ n + +(* [resolv_var v rst tbl] returns (if it exists) a restricted variable vr such that v = vr *) +exception FoundVar of int + +let resolve_var v rst tbl = + let v = Vect.set v (Int 1) Vect.null in + try + IMap.iter (fun k vect -> + if Restricted.is_restricted k rst + then if Vect.equal v vect then raise (FoundVar k) + else ()) tbl ; None + with FoundVar k -> Some k + +let prepare_cut env rst tbl x v = + (* extract the unrestricted part *) + let (unrst,rstv) = Vect.partition (fun x vl -> not (Restricted.is_restricted x rst) && frac_num vl <>/ Int 0) (Vect.set 0 (Int 0) v) in + if Vect.is_null unrst + then Some rstv + else Some (Vect.fold (fun acc k i -> + match resolve_var k rst tbl with + | None -> acc (* Should not happen *) + | Some v' -> Vect.set v' i acc) + rstv unrst) + +let cut env rmin sol vm (rst:Restricted.t) tbl (x,v) = + begin + (* Printf.printf "Trying to cut %i\n" x;*) + let (n,r) = Vect.decomp_cst v in + + + let f = frac_num n in + + if f =/ Int 0 then None (* The solution is integral *) else (* This is potentially a cut *) - let cut = Vect.normalise - (Vect.fold (fun acc x n -> - if Restricted.is_restricted x rst then - Vect.set x (n -/ (Num.floor_num n)) acc - else acc - ) Vect.null r) in - if debug then Printf.fprintf stdout "Cut vector for %a : %a\n" LinPoly.pp_var x LinPoly.pp cut ; - let cut = make_farkas_proof env vm cut in - - match WithProof.cutting_plane cut with - | None -> None - | Some (v,prf) -> - if debug then begin - Printf.printf "This is a cutting plane:\n" ; - Printf.printf "%a -> %a\n" WithProof.output cut WithProof.output (v,prf); - end; - if Pervasives.(=) (snd v) Eq - then (* Unsat *) Some (x,(v,prf)) - else if eval_op Ge (Vect.dotproduct (fst v) (Vect.set 0 (Int 1) sol)) (Int 0) - then begin - (* Can this happen? *) - if debug then Printf.printf "The cut is feasible - drop it\n"; - None - end - else Some(x,(v,prf)) - -let find_cut env u sol vm rst tbl = - (* find first *) - IMap.fold (fun x v acc -> - match acc with - | None -> cut env u sol vm rst (x,v) - | Some c -> acc) tbl None - -(* -let find_cut env u sol vm rst tbl = - IMap.fold (fun x v acc -> - match acc with - | Some c -> Some c - | None -> cut env u sol vm rst (x,v) - ) tbl None - *) + let t = + if f Vect.null + | Some r -> + (*Printf.printf "Potential cut %a\n" LinPoly.pp r;*) + Vect.fold (fun acc x n -> Vect.set x (ccoeff n) acc) Vect.null r + in + + let lcut = List.map (fun cv -> Vect.normalise (cut_vector cv)) [cut_coeff1 ; cut_coeff2] in + + let lcut = List.map (make_farkas_proof env vm) lcut in + + let check_cutting_plane c = + match WithProof.cutting_plane c with + | None -> + if debug then Printf.printf "This is not cutting plane for %a\n%a:" LinPoly.pp_var x WithProof.output c; + None + | Some(v,prf) -> + if debug then begin + Printf.printf "This is a cutting plane for %a:" LinPoly.pp_var x; + Printf.printf " %a\n" WithProof.output (v,prf); + end; + if Pervasives.(=) (snd v) Eq + then (* Unsat *) Some (x,(v,prf)) + else + let vl = (Vect.dotproduct (fst v) (Vect.set 0 (Int 1) sol)) in + if eval_op Ge vl (Int 0) + then begin + (* Can this happen? *) + if debug then Printf.printf "The cut is feasible %s >= 0 ??\n" (Num.string_of_num vl); + None + end + else Some(x,(v,prf)) in + + find_some check_cutting_plane lcut + end + +let find_cut nb env u sol vm rst tbl = + if nb = 0 + then + IMap.fold (fun x v acc -> + match acc with + | None -> cut env u sol vm rst tbl (x,v) + | Some c -> Some c) tbl None + else + IMap.fold (fun x v acc -> + match cut env u sol vm rst tbl (x,v) , acc with + | None , Some r | Some r , None -> Some r + | None , None -> None + | Some (v,((lp,o),p1)) , Some (v',((lp',o'),p2)) -> + Some (if ProofFormat.pr_size p1 Inl (Restricted.restrict vr rst,t',x) | Unsat c -> Inr c in + let nb = ref 0 in + let rec isolve env cr vr res = + incr nb; match res with | Inr c -> Some (Step (vr, make_farkas_certificate env vm (Vect.normalise c),Done)) | Inl (rst,tbl,x) -> @@ -595,10 +644,11 @@ let integer_solver lp = Printf.fprintf stdout "Looking for a cut\n"; Printf.fprintf stdout "Restricted %a ...\n" Restricted.pp rst; Printf.fprintf stdout "Current Tableau\n%a\n" output_tableau tbl; + (* Printf.fprintf stdout "Bounding box\n%a\n" output_box (bounding_box (vr+1) rst tbl l')*) end; let sol = find_solution rst tbl in - match find_cut env cr (*x*) sol vm rst tbl with + match find_cut (!nb mod 2) env cr (*x*) sol vm rst tbl with | None -> None | Some(cr,((v,op),cut)) -> if Pervasives.(=) op Eq @@ -615,6 +665,8 @@ let integer_solver lp = isolve env None vr res let integer_solver lp = + if debug then Printf.printf "Input integer solver\n%a\n" WithProof.output_sys (List.map WithProof.of_cstr lp); + match integer_solver lp with | None -> None | Some prf -> if debug -- cgit v1.2.3