diff options
Diffstat (limited to 'plugins/micromega/simplex.ml')
| -rw-r--r-- | plugins/micromega/simplex.ml | 221 |
1 files changed, 140 insertions, 81 deletions
diff --git a/plugins/micromega/simplex.ml b/plugins/micromega/simplex.ml index ade8143f3c..eaa26ded62 100644 --- a/plugins/micromega/simplex.ml +++ b/plugins/micromega/simplex.ml @@ -1,28 +1,70 @@ (************************************************************************) (* * The Coq Proof Assistant / The Coq Development Team *) -(* v * INRIA, CNRS and contributors - Copyright 1999-2019 *) -(* <O___,, * (see CREDITS file for the list of authors) *) +(* v * Copyright INRIA, CNRS and contributors *) +(* <O___,, * (see version control and CREDITS file for authors & dates) *) (* \VV/ **************************************************************) (* // * This file is distributed under the terms of the *) (* * GNU Lesser General Public License Version 2.1 *) (* * (see LICENSE file for the text of the license) *) (************************************************************************) +open NumCompat +open Q.Notations open Polynomial -open Num - -(*open Util*) open Mutils 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 (** Mapping basic variables to their equation. All variables >= than a threshold rst are restricted.*) +type tableau = Vect.t IMap.t module Restricted = struct type t = @@ -60,10 +102,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 @@ -78,7 +117,7 @@ let output_vars o m = let unfeasible (rst : Restricted.t) tbl = Restricted.fold rst - (fun k v m -> if Vect.get_cst v >=/ Int 0 then m else IMap.add k () m) + (fun k v m -> if Vect.get_cst v >=/ Q.zero then m else IMap.add k () m) tbl IMap.empty let is_feasible rst tb = IMap.is_empty (unfeasible rst tb) @@ -98,7 +137,7 @@ let is_feasible rst tb = IMap.is_empty (unfeasible rst tb) let is_maximised_vect rst v = Vect.for_all (fun xi ai -> - if ai >/ Int 0 then false else Restricted.is_restricted xi rst) + if ai >/ Q.zero then false else Restricted.is_restricted xi rst) v (** [is_maximised rst v] @@ -121,11 +160,11 @@ let is_maximised rst v = *) type result = - | Max of num (** Maximum is reached *) + | Max of Q.t (** Maximum is reached *) | Ubnd of var (** Problem is unbounded *) | Feas (** Problem is feasible *) -type pivot = Done of result | Pivot of int * int * num +type pivot = Done of result | Pivot of int * int * Q.t type simplex = Opt of tableau * result (** For a row, x = ao.xo+...+ai.xi @@ -140,7 +179,7 @@ let rec find_pivot_column (rst : Restricted.t) (r : Vect.t) = match Vect.choose r with | None -> failwith "find_pivot_column" | Some (xi, ai, r') -> - if ai </ Int 0 then + if ai </ Q.zero then if Restricted.is_restricted xi rst then find_pivot_column rst r' (* ai.xi cannot be improved *) else (xi, -1) (* r is not restricted, sign of ai does not matter *) @@ -167,9 +206,9 @@ let find_pivot_row rst tbl j sgn = Restricted.fold rst (fun i' v res -> let aij = Vect.get j v in - if Int sgn */ aij </ Int 0 then + if Q.of_int sgn */ aij </ Q.zero then (* This would improve *) - let score' = Num.abs_num (Vect.get_cst v // aij) in + let score' = Q.abs (Vect.get_cst v // aij) in min_score res (i', score') else res) tbl None @@ -206,10 +245,10 @@ let find_pivot vr (rst : Restricted.t) tbl = let solve_column (c : var) (r : var) (e : Vect.t) : Vect.t = let a = Vect.get c e in - if a =/ Int 0 then failwith "Cannot solve column" + if a =/ Q.zero then failwith "Cannot solve column" else - let a' = Int (-1) // a in - Vect.mul a' (Vect.set r (Int (-1)) (Vect.set c (Int 0) e)) + let a' = Q.neg_one // a in + Vect.mul a' (Vect.set r Q.neg_one (Vect.set c Q.zero e)) (** [pivot_row r c e] @param c is such that c = e @@ -218,18 +257,19 @@ let solve_column (c : var) (r : var) (e : Vect.t) : Vect.t = let pivot_row (row : Vect.t) (c : var) (e : Vect.t) : Vect.t = let g = Vect.get c row in - if g =/ Int 0 then row else Vect.mul_add g e (Int 1) (Vect.set c (Int 0) row) + if g =/ Q.zero then row else Vect.mul_add g e Q.one (Vect.set c Q.zero row) 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) let adapt_unbounded vr x rst tbl = - if Vect.get_cst (IMap.find vr tbl) >=/ Int 0 then tbl else pivot tbl vr x + if Vect.get_cst (IMap.find vr tbl) >=/ Q.zero then tbl else pivot tbl vr x module BaseSet = Set.Make (struct type t = iset @@ -254,7 +294,7 @@ let simplex opt vr rst tbl = output_tableau stdout tbl; Printf.fprintf stdout "Error for variables %a\n" output_vars m end; - if (not opt) && Vect.get_cst (IMap.find vr tbl) >=/ Int 0 then + if (not opt) && Vect.get_cst (IMap.find vr tbl) >=/ Q.zero then Opt (tbl, Feas) else match find_pivot vr rst tbl with @@ -267,7 +307,7 @@ let simplex opt vr rst tbl = | Feas -> raise (Invalid_argument "find_pivot") ) | Pivot (i, j, s) -> if debug then begin - Printf.fprintf stdout "Find pivot for x%i(%s)\n" vr (string_of_num s); + Printf.fprintf stdout "Find pivot for x%i(%s)\n" vr (Q.to_string s); Printf.fprintf stdout "Leaving variable x%i\n" i; Printf.fprintf stdout "Entering variable x%i\n" j end; @@ -318,18 +358,17 @@ let push_real (opt : bool) (nw : var) (v : Vect.t) (rst : Restricted.t) | Feas -> Sat (t', None) | Max n -> if debug then begin - Printf.printf "The objective is maximised %s\n" (string_of_num n); + Printf.printf "The objective is maximised %s\n" (Q.to_string n); Printf.printf "%a = %a\n" LinPoly.pp_var nw pp_row (IMap.find nw t') end; - if n >=/ Int 0 then Sat (t', None) + if n >=/ Q.zero then Sat (t', None) else let v' = safe_find "push_real" nw t' in - Unsat - (Vect.set nw (Int 1) (Vect.set 0 (Int 0) (Vect.mul (Int (-1)) v'))) ) + Unsat (Vect.set nw Q.one (Vect.set 0 Q.zero (Vect.mul Q.neg_one v'))) ) -open Mutils (** One complication is that equalities needs some pre-processing. *) +open Mutils open Polynomial @@ -340,7 +379,7 @@ let make_certificate vm l = (Vect.fold (fun acc x n -> let x', b = IMap.find x vm in - Vect.set x' (if b then n else Num.minus_num n) acc) + Vect.set x' (if b then n else Q.neg n) acc) Vect.null l) (** [eliminate_equalities vr0 l] @@ -356,11 +395,11 @@ let eliminate_equalities (vr0 : var) (l : Polynomial.cstr list) = | c :: l -> ( match c.op with | Ge -> - let v = Vect.set 0 (minus_num c.cst) c.coeffs in + let v = Vect.set 0 (Q.neg c.cst) c.coeffs in elim (idx + 1) (vr + 1) (IMap.add vr (idx, true) vm) l ((vr, v) :: acc) | Eq -> - let v1 = Vect.set 0 (minus_num c.cst) c.coeffs in - let v2 = Vect.mul (Int (-1)) v1 in + let v1 = Vect.set 0 (Q.neg c.cst) c.coeffs in + let v2 = Vect.mul Q.neg_one v1 in let vm = IMap.add vr (idx, true) (IMap.add (vr + 1) (idx, false) vm) in elim (idx + 1) (vr + 2) vm l ((vr, v1) :: (vr + 1, v2) :: acc) | Gt -> raise Strict ) @@ -378,7 +417,7 @@ let find_full_solution rst tbl = IMap.fold (fun vr v res -> Vect.set vr (Vect.get_cst v) res) tbl Vect.null let choose_conflict (sol : Vect.t) (l : (var * Vect.t) list) = - let esol = Vect.set 0 (Int 1) sol in + let esol = Vect.set 0 Q.one sol in let rec most_violating l e (x, v) rst = match l with | [] -> Some ((x, v), rst) @@ -435,7 +474,7 @@ let optimise obj l = let _, vm, l' = eliminate_equalities (vr0 + 1) l in let bound pos res = match res with - | Opt (_, Max n) -> Some (if pos then n else minus_num n) + | Opt (_, Max n) -> Some (if pos then n else Q.neg n) | Opt (_, Ubnd _) -> None | Opt (_, Feas) -> None in @@ -460,9 +499,7 @@ let make_farkas_certificate (env : WithProof.t IMap.t) vm v = begin try let x', b = IMap.find x vm in - mul_cst_proof - (if b then n else Num.minus_num n) - (snd (IMap.find x' env)) + mul_cst_proof (if b then n else Q.neg n) (snd (IMap.find x' env)) with Not_found -> (* This is an introduced hypothesis *) mul_cst_proof n (snd (IMap.find x env)) @@ -476,13 +513,16 @@ let make_farkas_proof (env : WithProof.t IMap.t) vm v = begin 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 n = if b then n else Q.neg n in + 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 -let frac_num n = n -/ Num.floor_num n +let frac_num n = n -/ Q.floor n type ('a, 'b) hitkind = | Forget @@ -493,21 +533,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 =/ Q.zero 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 = - let fv = frac_num v in - if fv <=/ Int 1 -/ f then fv // (Int 1 -/ f) else (Int 1 -/ fv) // f + (* 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 + Num.min_num fv (fn */ (Q.one -/ fv) // (Q.one -/ fn)) + in + let ccoeff_prop3 = + if Restricted.is_restricted x rst then ("Prop3", ccoeff_prop3) + else ("Prop1", ccoeff_prop1) + in *) + let n0_5 = Q.one // Q.two 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 *) + Q.ceiling (n0_5 // fn) + else + (* multiply by -1 & t-cut *) + Q.neg (Q.ceiling (n0_5 // (Q.one -/ fn))) + in + ("Prop2", fun v -> frac_num (v */ tmin)) in - let cut_coeff2 v = frac_num (t */ v) in + let ccoeff = ccoeff_prop2 in let cut_vector ccoeff = Vect.fold (fun acc x n -> @@ -516,35 +578,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 @@ -589,7 +647,7 @@ let eliminate_variable (bounded, vr, env, tbl) x = let tv = var_of_vect t in (* x = z - t *) let xdef = Vect.add z (Vect.uminus t) in - let xp = ((Vect.set x (Int 1) (Vect.uminus xdef), Eq), Def vr) in + let xp = ((Vect.set x Q.one (Vect.uminus xdef), Eq), Def vr) in let zp = ((z, Ge), Def zv) in let tp = ((t, Ge), Def tv) in (* Pivot the current tableau using xdef *) @@ -600,11 +658,8 @@ let eliminate_variable (bounded, vr, env, tbl) x = (fun lp -> let (v, o), p = lp in let ai = Vect.get x v in - if ai =/ Int 0 then lp - else - WithProof.addition - (WithProof.mult (Vect.cst (Num.minus_num ai)) xp) - lp) + if ai =/ Q.zero then lp + else WithProof.addition (WithProof.mult (Vect.cst (Q.neg ai)) xp) lp) env in (* Add the variables to the environment *) @@ -681,12 +736,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 |