From ccff8bf591ee0b8499cdc9dc9bb2827e4d2dae69 Mon Sep 17 00:00:00 2001
From: Frédéric Besson
Date: Tue, 28 Jan 2020 15:33:22 +0100
Subject: Fix efficiency regression #11436

- The cutting plane has been (sometimes) improved to generate stronger
  cuts.
- There is now some support for profiling the Simplex
  see documentation for Show Lia Profile.
---
 plugins/micromega/coq_micromega.ml  |  30 ++++++++
 plugins/micromega/coq_micromega.mli |   2 +-
 plugins/micromega/g_micromega.mlg   |   7 +-
 plugins/micromega/simplex.ml        | 142 ++++++++++++++++++++++++++----------
 plugins/micromega/simplex.mli       |  14 ++++
 5 files changed, 152 insertions(+), 43 deletions(-)

(limited to 'plugins')

diff --git a/plugins/micromega/coq_micromega.ml b/plugins/micromega/coq_micromega.ml
index 92a2222cfa..cdadde8621 100644
--- a/plugins/micromega/coq_micromega.ml
+++ b/plugins/micromega/coq_micromega.ml
@@ -2416,6 +2416,36 @@ let nqa =
     (fun _ x -> x)
     Mc.cnfQ qq_domain_spec dump_qexpr nlinear_prover_R
 
+let print_lia_profile () =
+  Simplex.(
+    let { number_of_successes
+        ; number_of_failures
+        ; success_pivots
+        ; failure_pivots
+        ; average_pivots
+        ; maximum_pivots } =
+      Simplex.get_profile_info ()
+    in
+    Feedback.msg_notice
+      Pp.(
+        (* successes *)
+        str "number of successes: "
+        ++ int number_of_successes ++ fnl ()
+        (* success pivots *)
+        ++ str "number of success pivots: "
+        ++ int success_pivots ++ fnl ()
+        (* failure *)
+        ++ str "number of failures: "
+        ++ int number_of_failures ++ fnl ()
+        (* failure pivots *)
+        ++ str "number of failure pivots: "
+        ++ int failure_pivots ++ fnl ()
+        (* Other *)
+        ++ str "average number of pivots: "
+        ++ int average_pivots ++ fnl ()
+        ++ str "maximum number of pivots: "
+        ++ int maximum_pivots ++ fnl ()))
+
 (* Local Variables: *)
 (* coding: utf-8 *)
 (* End: *)
diff --git a/plugins/micromega/coq_micromega.mli b/plugins/micromega/coq_micromega.mli
index 37ea560241..bcfc47357f 100644
--- a/plugins/micromega/coq_micromega.mli
+++ b/plugins/micromega/coq_micromega.mli
@@ -8,7 +8,6 @@
 (*         *     (see LICENSE file for the text of the license)         *)
 (************************************************************************)
 
-(*val is_ground_tac : EConstr.constr -> unit Proofview.tactic*)
 val psatz_Z : int -> unit Proofview.tactic -> unit Proofview.tactic
 val psatz_Q : int -> unit Proofview.tactic -> unit Proofview.tactic
 val psatz_R : int -> unit Proofview.tactic -> unit Proofview.tactic
@@ -21,6 +20,7 @@ val sos_Q : unit Proofview.tactic -> unit Proofview.tactic
 val sos_R : unit Proofview.tactic -> unit Proofview.tactic
 val lra_Q : unit Proofview.tactic -> unit Proofview.tactic
 val lra_R : unit Proofview.tactic -> unit Proofview.tactic
+val print_lia_profile : unit -> unit
 
 (** {5 Use Micromega independently from tactics. } *)
 
diff --git a/plugins/micromega/g_micromega.mlg b/plugins/micromega/g_micromega.mlg
index edf8106f30..d0f70bceac 100644
--- a/plugins/micromega/g_micromega.mlg
+++ b/plugins/micromega/g_micromega.mlg
@@ -28,10 +28,6 @@ open Tacarg
 
 DECLARE PLUGIN "micromega_plugin"
 
-TACTIC EXTEND RED
-| [ "myred"  ] -> { Tactics.red_in_concl }
-END
-
 TACTIC EXTEND PsatzZ
 | [ "psatz_Z" int_or_var(i) tactic(t) ] -> {  (Coq_micromega.psatz_Z i
                                                (Tacinterp.tactic_of_value ist t))
@@ -87,3 +83,6 @@ TACTIC EXTEND PsatzQ
 | [ "psatz_Q" tactic(t) ] -> {  (Coq_micromega.psatz_Q (-1) (Tacinterp.tactic_of_value ist t)) }
 END
 
+VERNAC COMMAND EXTEND ShowLiaProfile CLASSIFIED AS QUERY
+| [ "Show" "Lia" "Profile" ] -> { Coq_micromega.print_lia_profile () }
+END
diff --git a/plugins/micromega/simplex.ml b/plugins/micromega/simplex.ml
index ade8143f3c..54976221bc 100644
--- a/plugins/micromega/simplex.ml
+++ b/plugins/micromega/simplex.ml
@@ -18,6 +18,49 @@ type ('a, 'b) sum = Inl of 'a | Inr of 'b
 
 let debug = false
 
+(** Exploiting profiling information *)
+
+let profile_info = ref []
+let nb_pivot = ref 0
+
+type profile_info =
+  { number_of_successes : int
+  ; number_of_failures : int
+  ; success_pivots : int
+  ; failure_pivots : int
+  ; average_pivots : int
+  ; maximum_pivots : int }
+
+let init_profile =
+  { number_of_successes = 0
+  ; number_of_failures = 0
+  ; success_pivots = 0
+  ; failure_pivots = 0
+  ; average_pivots = 0
+  ; maximum_pivots = 0 }
+
+let get_profile_info () =
+  let update_profile
+      { number_of_successes
+      ; number_of_failures
+      ; success_pivots
+      ; failure_pivots
+      ; average_pivots
+      ; maximum_pivots } (b, i) =
+    { number_of_successes = (number_of_successes + if b then 1 else 0)
+    ; number_of_failures = (number_of_failures + if b then 0 else 1)
+    ; success_pivots = (success_pivots + if b then i else 0)
+    ; failure_pivots = (failure_pivots + if b then 0 else i)
+    ; average_pivots = average_pivots + 1 (* number of proofs *)
+    ; maximum_pivots = max maximum_pivots i }
+  in
+  let p = List.fold_left update_profile init_profile !profile_info in
+  profile_info := [];
+  { p with
+    average_pivots =
+      ( try (p.success_pivots + p.failure_pivots) / p.average_pivots
+        with Division_by_zero -> 0 ) }
+
 type iset = unit IMap.t
 
 type tableau = Vect.t IMap.t
@@ -60,10 +103,7 @@ let output_tableau o t =
     t
 
 let output_env o t =
-  IMap.iter
-    (fun k v ->
-      Printf.fprintf o "%a : %a\n" LinPoly.pp_var k WithProof.output v)
-    t
+  IMap.iter (fun k v -> Printf.fprintf o "%i : %a\n" k WithProof.output v) t
 
 let output_vars o m =
   IMap.iter (fun k _ -> Printf.fprintf o "%a " LinPoly.pp_var k) m
@@ -224,6 +264,7 @@ let pivot_with (m : tableau) (v : var) (p : Vect.t) =
   IMap.map (fun (r : Vect.t) -> pivot_row r v p) m
 
 let pivot (m : tableau) (r : var) (c : var) =
+  incr nb_pivot;
   let row = safe_find "pivot" r m in
   let piv = solve_column c r row in
   IMap.add c piv (pivot_with (IMap.remove r m) c piv)
@@ -477,8 +518,11 @@ let make_farkas_proof (env : WithProof.t IMap.t) vm v =
           try
             let x', b = IMap.find x vm in
             let n = if b then n else Num.minus_num n in
-            WithProof.mult (Vect.cst n) (IMap.find x' env)
-          with Not_found -> WithProof.mult (Vect.cst n) (IMap.find x env)
+            let prf = IMap.find x' env in
+            WithProof.mult (Vect.cst n) prf
+          with Not_found ->
+            let prf = IMap.find x env in
+            WithProof.mult (Vect.cst n) prf
         end)
     WithProof.zero v
 
@@ -493,21 +537,43 @@ type ('a, 'b) hitkind =
 
 let cut env rmin sol vm (rst : Restricted.t) tbl (x, v) =
   let n, r = Vect.decomp_cst v in
-  let f = frac_num n in
-  if f =/ Int 0 then Forget (* The solution is integral *)
+  let fn = frac_num n in
+  if fn =/ Int 0 then Forget (* The solution is integral *)
   else
-    (* This is potentially a cut *)
-    let t =
-      if f </ Int 1 // Int 2 then
-        let t' = Int 1 // f in
-        if Num.is_integer_num t' then t' -/ Int 1 else Num.floor_num t'
-      else Int 1
-    in
-    let cut_coeff1 v =
+    (* The cut construction is from:
+       Letchford and Lodi. Strengthening Chvatal-Gomory cuts and Gomory fractional cuts.
+
+       We implement the classic Proposition 2 from the "known results"
+     *)
+
+    (* Proposition 3 requires all the variables to be restricted and is
+       therefore not always applicable. *)
+    (* let ccoeff_prop1 v = frac_num v in
+    let ccoeff_prop3 v =
+      (* mixed integer cut *)
       let fv = frac_num v in
-      if fv <=/ Int 1 -/ f then fv // (Int 1 -/ f) else (Int 1 -/ fv) // f
+      Num.min_num fv (fn */ (Int 1 -/ fv) // (Int 1 -/ fn))
     in
-    let cut_coeff2 v = frac_num (t */ v) in
+    let ccoeff_prop3 =
+      if Restricted.is_restricted x rst then ("Prop3", ccoeff_prop3)
+      else ("Prop1", ccoeff_prop1)
+    in *)
+    let n0_5 = Int 1 // Int 2 in
+    (* If the fractional part [fn] is small, we construct the t-cut.
+       If the fractional part [fn] is big, we construct the t-cut of the negated row.
+       (This is only a cut if all the fractional variables are restricted.)
+     *)
+    let ccoeff_prop2 =
+      let tmin =
+        if fn </ n0_5 then (* t-cut *)
+          Num.ceiling_num (n0_5 // fn)
+        else
+          (* multiply by -1 & t-cut *)
+          minus_num (Num.ceiling_num (n0_5 // (Int 1 -/ fn)))
+      in
+      ("Prop2", fun v -> frac_num (v */ tmin))
+    in
+    let ccoeff = ccoeff_prop2 in
     let cut_vector ccoeff =
       Vect.fold
         (fun acc x n ->
@@ -516,35 +582,31 @@ let cut env rmin sol vm (rst : Restricted.t) tbl (x, v) =
         Vect.null r
     in
     let lcut =
-      List.map
-        (fun cv -> Vect.normalise (cut_vector cv))
-        [cut_coeff1; cut_coeff2]
+      ( fst ccoeff
+      , make_farkas_proof env vm (Vect.normalise (cut_vector (snd ccoeff))) )
     in
-    let lcut = List.map (make_farkas_proof env vm) lcut in
-    let check_cutting_plane c =
+    let check_cutting_plane (p, c) =
       match WithProof.cutting_plane c with
       | None ->
         if debug then
-          Printf.printf "This is not a cutting plane for %a\n%a:" LinPoly.pp_var
-            x WithProof.output c;
+          Printf.printf "%s: This is not a cutting plane for %a\n%a:" p
+            LinPoly.pp_var x WithProof.output c;
         None
       | Some (v, prf) ->
         if debug then (
-          Printf.printf "This is a cutting plane for %a:" LinPoly.pp_var x;
+          Printf.printf "%s: This is a cutting plane for %a:" p LinPoly.pp_var x;
           Printf.printf " %a\n" WithProof.output (v, prf) );
-        if 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 (
-            if debug then
-              Printf.printf "The cut is feasible %s >= 0 \n"
-                (Num.string_of_num vl);
-            None )
-          else Some (x, (v, prf))
+        Some (x, (v, prf))
     in
-    match find_some check_cutting_plane lcut with
+    match check_cutting_plane lcut with
     | Some r -> Hit r
-    | None -> Keep (x, v)
+    | None ->
+      let has_unrestricted =
+        Vect.fold
+          (fun acc v vl -> acc || not (Restricted.is_restricted v rst))
+          false r
+      in
+      if has_unrestricted then Keep (x, v) else Forget
 
 let merge_result_old oldr f x =
   match oldr with
@@ -681,12 +743,16 @@ let integer_solver lp =
   isolve env None vr res
 
 let integer_solver lp =
+  nb_pivot := 0;
   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
+  | None ->
+    profile_info := (false, !nb_pivot) :: !profile_info;
+    None
   | Some prf ->
+    profile_info := (true, !nb_pivot) :: !profile_info;
     if debug then
       Printf.fprintf stdout "Proof %a\n" ProofFormat.output_proof prf;
     Some prf
diff --git a/plugins/micromega/simplex.mli b/plugins/micromega/simplex.mli
index 19bcce3590..ff672edafd 100644
--- a/plugins/micromega/simplex.mli
+++ b/plugins/micromega/simplex.mli
@@ -9,6 +9,20 @@
 (************************************************************************)
 open Polynomial
 
+(** Profiling *)
+
+type profile_info =
+  { number_of_successes : int
+  ; number_of_failures : int
+  ; success_pivots : int
+  ; failure_pivots : int
+  ; average_pivots : int
+  ; maximum_pivots : int }
+
+val get_profile_info : unit -> profile_info
+
+(** Simplex interface *)
+
 val optimise : Vect.t -> cstr list -> (Num.num option * Num.num option) option
 val find_point : cstr list -> Vect.t option
 val find_unsat_certificate : cstr list -> Vect.t option
-- 
cgit v1.2.3