aboutsummaryrefslogtreecommitdiff
path: root/plugins/micromega
diff options
context:
space:
mode:
Diffstat (limited to 'plugins/micromega')
-rw-r--r--plugins/micromega/ZifyInst.v19
-rw-r--r--plugins/micromega/certificate.ml72
-rw-r--r--plugins/micromega/coq_micromega.ml37
-rw-r--r--plugins/micromega/coq_micromega.mli2
-rw-r--r--plugins/micromega/g_micromega.mlg7
-rw-r--r--plugins/micromega/mutils.ml19
-rw-r--r--plugins/micromega/mutils.mli1
-rw-r--r--plugins/micromega/polynomial.ml32
-rw-r--r--plugins/micromega/polynomial.mli3
-rw-r--r--plugins/micromega/simplex.ml142
-rw-r--r--plugins/micromega/simplex.mli14
11 files changed, 255 insertions, 93 deletions
diff --git a/plugins/micromega/ZifyInst.v b/plugins/micromega/ZifyInst.v
index 97f6fe0613..edfb5a2a94 100644
--- a/plugins/micromega/ZifyInst.v
+++ b/plugins/micromega/ZifyInst.v
@@ -523,3 +523,22 @@ Instance SatProdPos : Saturate Z.mul :=
SatOk := Z.mul_pos_pos
|}.
Add Saturate SatProdPos.
+
+Lemma pow_pos_strict :
+ forall a b,
+ 0 < a -> 0 < b -> 0 < a ^ b.
+Proof.
+ intros.
+ apply Z.pow_pos_nonneg; auto.
+ apply Z.lt_le_incl;auto.
+Qed.
+
+
+Instance SatPowPos : Saturate Z.pow :=
+ {|
+ PArg1 := fun x => 0 < x;
+ PArg2 := fun y => 0 < y;
+ PRes := fun r => 0 < r;
+ SatOk := pow_pos_strict
+ |}.
+Add Saturate SatPowPos.
diff --git a/plugins/micromega/certificate.ml b/plugins/micromega/certificate.ml
index cb15274736..61234145e1 100644
--- a/plugins/micromega/certificate.ml
+++ b/plugins/micromega/certificate.ml
@@ -395,50 +395,40 @@ let saturate_by_linear_equalities sys =
output_sys sys output_sys sys';
sys'
-(* let saturate_linear_equality_non_linear sys0 =
- let (l,_) = extract_all (is_substitution false) sys0 in
- let rec elim l acc =
- match l with
- | [] -> acc
- | (v,pc)::l' ->
- let nc = saturate (non_linear_pivot sys0 pc v) (sys0@acc) in
- elim l' (nc@acc) in
- elim l []
- *)
-
-let bounded_vars (sys : WithProof.t list) =
- let l = fst (extract_all (fun ((p, o), prf) -> LinPoly.is_variable p) sys) in
- List.fold_left (fun acc (i, wp) -> IMap.add i wp acc) IMap.empty l
-
-let rec power n p = if n = 1 then p else WithProof.product p (power (n - 1) p)
-
-let bound_monomial mp m =
- if Monomial.is_var m || Monomial.is_const m then None
- else
- try
- Some
- (Monomial.fold
- (fun v i acc ->
- let wp = IMap.find v mp in
- WithProof.product (power i wp) acc)
- m (WithProof.const (Int 1)))
- with Not_found -> None
-
let bound_monomials (sys : WithProof.t list) =
- let mp = bounded_vars sys in
- let m =
+ let l =
+ extract_all
+ (fun ((p, o), _) ->
+ match LinPoly.get_bound p with
+ | None -> None
+ | Some Vect.Bound.{cst; var; coeff} ->
+ Some (Monomial.degree (LinPoly.MonT.retrieve var)))
+ sys
+ in
+ let deg =
+ List.fold_left (fun acc ((p, o), _) -> max acc (LinPoly.degree p)) 0 sys
+ in
+ let vars =
List.fold_left
- (fun acc ((p, _), _) ->
- Vect.fold
- (fun acc v _ ->
- let m = LinPoly.MonT.retrieve v in
- match bound_monomial mp m with
- | None -> acc
- | Some r -> IMap.add v r acc)
- acc p)
- IMap.empty sys
+ (fun acc ((p, o), _) -> ISet.union (LinPoly.monomials p) acc)
+ ISet.empty sys
+ in
+ let bounds =
+ saturate_bin
+ (fun (i1, w1) (i2, w2) ->
+ if i1 + i2 > deg then None
+ else
+ match WithProof.mul_bound w1 w2 with
+ | None -> None
+ | Some b -> Some (i1 + i2, b))
+ (fst l)
+ in
+ let has_mon (_, ((p, o), _)) =
+ match LinPoly.get_bound p with
+ | None -> false
+ | Some Vect.Bound.{cst; var; coeff} -> ISet.mem var vars
in
- IMap.fold (fun _ e acc -> e :: acc) m []
+ List.map snd (List.filter has_mon bounds) @ snd l
let develop_constraints prfdepth n_spec sys =
LinPoly.MonT.clear ();
diff --git a/plugins/micromega/coq_micromega.ml b/plugins/micromega/coq_micromega.ml
index 92a2222cfa..4b656f8e61 100644
--- a/plugins/micromega/coq_micromega.ml
+++ b/plugins/micromega/coq_micromega.ml
@@ -55,7 +55,6 @@ let use_csdp_cache = ref true
let () =
let int_opt l vref =
{ optdepr = false
- ; optname = List.fold_right ( ^ ) l ""
; optkey = l
; optread = (fun () -> Some !vref)
; optwrite =
@@ -63,42 +62,36 @@ let () =
in
let lia_enum_opt =
{ optdepr = false
- ; optname = "Lia Enum"
; optkey = ["Lia"; "Enum"]
; optread = (fun () -> !lia_enum)
; optwrite = (fun x -> lia_enum := x) }
in
let solver_opt =
{ optdepr = false
- ; optname = "Use the Simplex instead of Fourier elimination"
; optkey = ["Simplex"]
; optread = (fun () -> !Certificate.use_simplex)
; optwrite = (fun x -> Certificate.use_simplex := x) }
in
let dump_file_opt =
{ optdepr = false
- ; optname = "Generate Coq goals in file from calls to 'lia' 'nia'"
; optkey = ["Dump"; "Arith"]
; optread = (fun () -> !Certificate.dump_file)
; optwrite = (fun x -> Certificate.dump_file := x) }
in
let lia_cache_opt =
{ optdepr = false
- ; optname = "cache of lia (.lia.cache)"
; optkey = ["Lia"; "Cache"]
; optread = (fun () -> !use_lia_cache)
; optwrite = (fun x -> use_lia_cache := x) }
in
let nia_cache_opt =
{ optdepr = false
- ; optname = "cache of nia (.nia.cache)"
; optkey = ["Nia"; "Cache"]
; optread = (fun () -> !use_nia_cache)
; optwrite = (fun x -> use_nia_cache := x) }
in
let nra_cache_opt =
{ optdepr = false
- ; optname = "cache of nra (.nra.cache)"
; optkey = ["Nra"; "Cache"]
; optread = (fun () -> !use_nra_cache)
; optwrite = (fun x -> use_nra_cache := x) }
@@ -2416,6 +2409,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/mutils.ml b/plugins/micromega/mutils.ml
index 03f042647c..160b492d3d 100644
--- a/plugins/micromega/mutils.ml
+++ b/plugins/micromega/mutils.ml
@@ -140,6 +140,25 @@ let saturate p f sys =
Printexc.print_backtrace stdout;
raise x
+let saturate_bin (f : 'a -> 'a -> 'a option) (l : 'a list) =
+ let rec map_with acc e l =
+ match l with
+ | [] -> acc
+ | e' :: l' -> (
+ match f e e' with
+ | None -> map_with acc e l'
+ | Some r -> map_with (r :: acc) e l' )
+ in
+ let rec map2_with acc l' =
+ match l' with [] -> acc | e' :: l' -> map2_with (map_with acc e' l) l'
+ in
+ let rec iterate acc l' =
+ match map2_with [] l' with
+ | [] -> List.rev_append l' acc
+ | res -> iterate (List.rev_append l' acc) res
+ in
+ iterate [] l
+
open Num
open Big_int
diff --git a/plugins/micromega/mutils.mli b/plugins/micromega/mutils.mli
index ef8d154b13..5dcaf3be44 100644
--- a/plugins/micromega/mutils.mli
+++ b/plugins/micromega/mutils.mli
@@ -116,6 +116,7 @@ val simplify : ('a -> 'a option) -> 'a list -> 'a list option
val saturate :
('a -> 'b option) -> ('b * 'a -> 'a -> 'a option) -> 'a list -> 'a list
+val saturate_bin : ('a -> 'a -> 'a option) -> 'a list -> 'a list
val generate : ('a -> 'b option) -> 'a list -> 'b list
val app_funs : ('a -> 'b option) list -> 'a -> 'b option
val command : string -> string array -> 'a -> 'b
diff --git a/plugins/micromega/polynomial.ml b/plugins/micromega/polynomial.ml
index a4f9b60b14..b20213979b 100644
--- a/plugins/micromega/polynomial.ml
+++ b/plugins/micromega/polynomial.ml
@@ -379,6 +379,8 @@ module LinPoly = struct
else acc)
[] l
+ let get_bound p = Vect.Bound.of_vect p
+
let min_list (l : int list) =
match l with [] -> None | e :: l -> Some (List.fold_left min e l)
@@ -892,8 +894,9 @@ module WithProof = struct
if Vect.is_null r && n >/ Int 0 then
((LinPoly.product p p1, o1), ProofFormat.mul_cst_proof n prf1)
else (
- Printf.printf "mult_error %a [*] %a\n" LinPoly.pp p output
- ((p1, o1), prf1);
+ if debug then
+ Printf.printf "mult_error %a [*] %a\n" LinPoly.pp p output
+ ((p1, o1), prf1);
raise InvalidProof )
let cutting_plane ((p, o), prf) =
@@ -1027,6 +1030,31 @@ module WithProof = struct
else None
in
saturate select gen sys0
+
+ open Vect.Bound
+
+ let mul_bound w1 w2 =
+ let (p1, o1), prf1 = w1 in
+ let (p2, o2), prf2 = w2 in
+ match (LinPoly.get_bound p1, LinPoly.get_bound p2) with
+ | None, _ | _, None -> None
+ | ( Some {cst = c1; var = v1; coeff = c1'}
+ , Some {cst = c2; var = v2; coeff = c2'} ) -> (
+ let good_coeff b o =
+ match o with
+ | Eq -> Some (minus_num b)
+ | _ -> if b <=/ Int 0 then Some (minus_num b) else None
+ in
+ match (good_coeff c1 o2, good_coeff c2 o1) with
+ | None, _ | _, None -> None
+ | Some c1, Some c2 ->
+ let ext_mult c w =
+ if c =/ Int 0 then zero else mult (LinPoly.constant c) w
+ in
+ Some
+ (addition
+ (addition (product w1 w2) (ext_mult c1 w2))
+ (ext_mult c2 w1)) )
end
(* Local Variables: *)
diff --git a/plugins/micromega/polynomial.mli b/plugins/micromega/polynomial.mli
index 7e905ac69b..4b56b037e0 100644
--- a/plugins/micromega/polynomial.mli
+++ b/plugins/micromega/polynomial.mli
@@ -224,6 +224,8 @@ module LinPoly : sig
p is linear in x i.e x does not occur in b and
a is a constant such that [pred a] *)
+ val get_bound : t -> Vect.Bound.t option
+
val product : t -> t -> t
(** [product p q]
@return the product of the polynomial [p*q] *)
@@ -372,4 +374,5 @@ module WithProof : sig
val saturate_subst : bool -> t list -> t list
val is_substitution : bool -> t -> var option
+ val mul_bound : t -> t -> t option
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
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-02 23:19:12 +0000