diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 0000000..1e43ca4 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,2 @@ +# Reformatting with OCamlFormat +b5bd3cc2645f45b55dd7addd69e47189e73316e9 diff --git a/.ocamlformat b/.ocamlformat new file mode 100644 index 0000000..0e0b64b --- /dev/null +++ b/.ocamlformat @@ -0,0 +1,2 @@ +version=0.26.0 +profile=default diff --git a/Makefile b/Makefile index 85d7d41..0942730 100644 --- a/Makefile +++ b/Makefile @@ -23,12 +23,10 @@ install: build uninstall: @dune uninstall -reindent: - @find src '(' -name '*.ml' -or -name '*.mli' ')' -print0 | xargs -0 echo "reindenting: " - @find src '(' -name '*.ml' -or -name '*.mli' ')' -print0 | xargs -0 sed -i 's/[[:space:]]*$$//' - @find src '(' -name '*.ml' -or -name '*.mli' ')' -print0 | xargs -0 ocp-indent -i +fmt: + dune build @fmt opam-deps: opam install . --deps-only -.PHONY: build doc all clean test watch install uninstall reindent opam-deps +.PHONY: build doc all clean test watch install uninstall fmt opam-deps diff --git a/dune-project b/dune-project index 86b2dcd..8530fcc 100644 --- a/dune-project +++ b/dune-project @@ -35,5 +35,3 @@ (zarith :with-test) (logs (>= 0.5.0)) )) - -; See the complete stanza docs at https://dune.readthedocs.io/en/stable/dune-files.html#dune-project diff --git a/extra/pre-commit--git-hook b/extra/pre-commit--git-hook new file mode 100644 index 0000000..399971e --- /dev/null +++ b/extra/pre-commit--git-hook @@ -0,0 +1,2 @@ +#!/bin/sh +exec dune build @fmt diff --git a/src/assertBounds.ml b/src/assertBounds.ml index 6e1d9ae..a7ccda4 100644 --- a/src/assertBounds.ml +++ b/src/assertBounds.ml @@ -5,15 +5,10 @@ (******************************************************************************) module type S = sig - module Core : CoreSig.S val var : - Core.t -> - ?min:Core.bound -> - ?max:Core.bound -> - Core.Var.t -> - Core.t * bool + Core.t -> ?min:Core.bound -> ?max:Core.bound -> Core.Var.t -> Core.t * bool val poly : Core.t -> @@ -22,260 +17,235 @@ module type S = sig ?max:Core.bound -> Core.Var.t -> Core.t * bool - end -module Make(Core : CoreSig.S) : S with module Core = Core = struct - - module Core = Core - module Result = Result.Make(Core) - open Core - - let empty_info value = { - mini = None; - maxi = None; - value; - vstatus = ValueOK; - empty_dom = false; - } - - (* If the environment works on integers and bounds are strict, - we shift the bound so that it is a large bound. Same goes - on rational bounds in an integer simplex. - Ex: - * x > 5 + Ɛ -> x >= 6 - * x > 4/3 -> x >= 2 - *) - let update_bound - (get_closer : R2.t -> R2.t) - (env : Core.t) - (bnd : bound option) = - match bnd with - | Some b when env.is_int -> - Some {b with bvalue = get_closer b.bvalue} - | other -> other - - let update_min_bound = update_bound R2.ceiling - let update_max_bound = update_bound R2.floor - - let new_status_basic stt fixme s info consistent_bnds = - let has_bad_value = info.vstatus != ValueOK in - match stt, consistent_bnds with - | UNSAT _, _ -> stt, fixme - | _, false -> - UNSAT s, SX.empty - - | UNK, true -> - stt, if has_bad_value then SX.add s fixme else SX.remove s fixme - - | SAT, _ -> +module Make (Core : CoreSig.S) : S with module Core = Core = struct + module Core = Core + module Result = Result.Make (Core) + open Core + + let empty_info value = + { mini = None; maxi = None; value; vstatus = ValueOK; empty_dom = false } + + (* If the environment works on integers and bounds are strict, + we shift the bound so that it is a large bound. Same goes + on rational bounds in an integer simplex. + Ex: + * x > 5 + Ɛ -> x >= 6 + * x > 4/3 -> x >= 2 + *) + let update_bound (get_closer : R2.t -> R2.t) (env : Core.t) + (bnd : bound option) = + match bnd with + | Some b when env.is_int -> Some { b with bvalue = get_closer b.bvalue } + | other -> other + + let update_min_bound = update_bound R2.ceiling + let update_max_bound = update_bound R2.floor + + let new_status_basic stt fixme s info consistent_bnds = + let has_bad_value = info.vstatus != ValueOK in + match (stt, consistent_bnds) with + | UNSAT _, _ -> (stt, fixme) + | _, false -> (UNSAT s, SX.empty) + | UNK, true -> + (stt, if has_bad_value then SX.add s fixme else SX.remove s fixme) + | SAT, _ -> assert (fixme == SX.empty); - if has_bad_value then UNK, SX.add s fixme else stt, fixme - - let assert_basic_var env x mini maxi = - let info, poly, changed = - try - let info, poly = MX.find x env.basic in - let info, chang1 = set_min_bound info mini in - let info, chang2 = set_max_bound info maxi in - info, poly, chang1 || chang2 - with Not_found -> assert false - in - let status, fixme = - new_status_basic env.status env.fixme x info (consistent_bounds info) - in - {env with basic = MX.add x (info, poly) env.basic; status; fixme}, - changed - - (* *) - - let new_status_non_basic x stt fixme ({mini; maxi; value; _} as info) = - match stt with - | UNSAT _ -> stt, fixme - | SAT | UNK when consistent_bounds info -> + if has_bad_value then (UNK, SX.add s fixme) else (stt, fixme) + + let assert_basic_var env x mini maxi = + let info, poly, changed = + try + let info, poly = MX.find x env.basic in + let info, chang1 = set_min_bound info mini in + let info, chang2 = set_max_bound info maxi in + (info, poly, chang1 || chang2) + with Not_found -> assert false + in + let status, fixme = + new_status_basic env.status env.fixme x info (consistent_bounds info) + in + ( { env with basic = MX.add x (info, poly) env.basic; status; fixme }, + changed ) + + (* *) + + let new_status_non_basic x stt fixme ({ mini; maxi; value; _ } as info) = + match stt with + | UNSAT _ -> (stt, fixme) + | (SAT | UNK) when consistent_bounds info -> assert (not (violates_min_bound value mini)); assert (not (violates_max_bound value maxi)); assert (stt != SAT || fixme == SX.empty); - stt, fixme - | SAT | UNK -> - UNSAT x, SX.empty - - let adapt_values_of_basic_vars env _old _new x use = - let {basic; _} = env in - let diff = R2.sub _new _old in - SX.fold - (fun s env -> - let info, p = try MX.find s basic with Not_found -> assert false in - let c_x = try P.find x p with Not_found -> assert false in - let info = - {info with value = R2.add info.value (R2.mult_by_const c_x diff)} - in - let info = ajust_status_of_basic info in - let status, fixme = - new_status_basic env.status env.fixme s info true - in - {env with status; fixme; basic = MX.add s (info, p) env.basic} - )use env - - let assert_non_basic_var env x mini maxi = - let info, use = - try MX.find x env.non_basic - with Not_found -> empty_info R2.zero, SX.empty - in - let info, chang1 = set_min_bound info mini in - let info, chang2 = set_max_bound info maxi in - let old_val = info.value in - let info, changed = ajust_value_of_non_basic info in - let status, fixme = new_status_non_basic x env.status env.fixme info in - let env = - {env with - non_basic = MX.add x (info, use) env.non_basic; status; fixme} - in - let env = - if not changed then env - else adapt_values_of_basic_vars env old_val info.value x use - in - env, chang1 || chang2 - - (* exported function: check_invariants called before and after *) - let var env ?min ?max x = - debug "[entry of assert_var]" env (Result.get None); - check_invariants env (Result.get None); - let mini = update_min_bound env min in - let maxi = update_max_bound env max in - let env, changed = - if MX.mem x env.basic then - assert_basic_var env x mini maxi - else - assert_non_basic_var env x mini maxi - in - debug "[exit of assert_var]" env (Result.get None); - check_invariants env (Result.get None); - env, changed - - let register_slake slk p env = - if MX.mem slk env.slake then env, false - else {env with slake = MX.add slk p env.slake}, true - - let update_use is_fresh_slk x_status slk use = - match x_status with - | P.Exists -> + (stt, fixme) + | SAT | UNK -> (UNSAT x, SX.empty) + + let adapt_values_of_basic_vars env _old _new x use = + let { basic; _ } = env in + let diff = R2.sub _new _old in + SX.fold + (fun s env -> + let info, p = try MX.find s basic with Not_found -> assert false in + let c_x = try P.find x p with Not_found -> assert false in + let info = + { info with value = R2.add info.value (R2.mult_by_const c_x diff) } + in + let info = ajust_status_of_basic info in + let status, fixme = new_status_basic env.status env.fixme s info true in + { env with status; fixme; basic = MX.add s (info, p) env.basic }) + use env + + let assert_non_basic_var env x mini maxi = + let info, use = + try MX.find x env.non_basic + with Not_found -> (empty_info R2.zero, SX.empty) + in + let info, chang1 = set_min_bound info mini in + let info, chang2 = set_max_bound info maxi in + let old_val = info.value in + let info, changed = ajust_value_of_non_basic info in + let status, fixme = new_status_non_basic x env.status env.fixme info in + let env = + { env with non_basic = MX.add x (info, use) env.non_basic; status; fixme } + in + let env = + if not changed then env + else adapt_values_of_basic_vars env old_val info.value x use + in + (env, chang1 || chang2) + + (* exported function: check_invariants called before and after *) + let var env ?min ?max x = + debug "[entry of assert_var]" env (Result.get None); + check_invariants env (Result.get None); + let mini = update_min_bound env min in + let maxi = update_max_bound env max in + let env, changed = + if MX.mem x env.basic then assert_basic_var env x mini maxi + else assert_non_basic_var env x mini maxi + in + debug "[exit of assert_var]" env (Result.get None); + check_invariants env (Result.get None); + (env, changed) + + let register_slake slk p env = + if MX.mem slk env.slake then (env, false) + else ({ env with slake = MX.add slk p env.slake }, true) + + let update_use is_fresh_slk x_status slk use = + match x_status with + | P.Exists -> assert (SX.mem slk use); use - - | P.Removed -> + | P.Removed -> assert (SX.mem slk use); SX.remove slk use - - | P.New -> - assert (not is_fresh_slk || not (SX.mem slk use)); + | P.New -> + assert ((not is_fresh_slk) || not (SX.mem slk use)); SX.add slk use - let update_use_list is_fresh modified_stt slk non_basic = - List.fold_left - (fun non_basic (x, x_status) -> - try - let i, u = MX.find x non_basic in - MX.add x (i, update_use is_fresh x_status slk u) non_basic - with Not_found -> assert false - )non_basic modified_stt - - let normalize_polynomial is_fresh slk p env = - P.fold - (fun x c (q, env) -> + let update_use_list is_fresh modified_stt slk non_basic = + List.fold_left + (fun non_basic (x, x_status) -> + try + let i, u = MX.find x non_basic in + MX.add x (i, update_use is_fresh x_status slk u) non_basic + with Not_found -> assert false) + non_basic modified_stt + + let normalize_polynomial is_fresh slk p env = + P.fold + (fun x c (q, env) -> + try + let info, use = MX.find x env.non_basic in + let new_q, x_status = P.accumulate x c q in + let use = update_use is_fresh x_status slk use in + (new_q, { env with non_basic = MX.add x (info, use) env.non_basic }) + with Not_found -> ( try - let info, use = MX.find x env.non_basic in - let new_q, x_status = P.accumulate x c q in - let use = update_use is_fresh x_status slk use in - new_q, {env with non_basic = MX.add x (info, use) env.non_basic} - + let _, p_of_x = MX.find x env.basic in + let new_q, modified_stt = P.append q c p_of_x in + ( new_q, + { + env with + non_basic = + update_use_list is_fresh modified_stt slk env.non_basic; + } ) with Not_found -> - try - let _ , p_of_x = MX.find x env.basic in - let new_q, modified_stt = P.append q c p_of_x in - new_q, - {env with - non_basic = - update_use_list is_fresh modified_stt slk env.non_basic} - - with Not_found -> - (* var not initied -> new non_basic *) - let env, chang = assert_non_basic_var env x None None in - assert (not chang); - let new_q, x_status = P.replace x c q in - assert (x_status == P.New); - let info, use = - try MX.find x env.non_basic - with Not_found -> assert false - in - let use = update_use is_fresh x_status slk use in - new_q, - { env with - non_basic = MX.add x (info, use) env.non_basic} - )p (P.empty, env) - - - (* exported function: check_invariants called before and after *) - let poly env p ?min ?max slk = - debug "[entry of assert_poly]" env (Result.get None); - check_invariants env (Result.get None); - assert (P.is_polynomial p); - let mini = update_min_bound env min in - let maxi = update_max_bound env max in - let env, is_fresh = register_slake slk p env in - let info, is_basic, env, change = - try (* non basic existing var ? *) - let info, use = MX.find slk env.non_basic in + (* var not initied -> new non_basic *) + let env, chang = assert_non_basic_var env x None None in + assert (not chang); + let new_q, x_status = P.replace x c q in + assert (x_status == P.New); + let info, use = + try MX.find x env.non_basic with Not_found -> assert false + in + let use = update_use is_fresh x_status slk use in + (new_q, { env with non_basic = MX.add x (info, use) env.non_basic }))) + p (P.empty, env) + + (* exported function: check_invariants called before and after *) + let poly env p ?min ?max slk = + debug "[entry of assert_poly]" env (Result.get None); + check_invariants env (Result.get None); + assert (P.is_polynomial p); + let mini = update_min_bound env min in + let maxi = update_max_bound env max in + let env, is_fresh = register_slake slk p env in + let info, is_basic, env, change = + try + (* non basic existing var ? *) + let info, use = MX.find slk env.non_basic in + assert ( + let np, _ = normalize_polynomial is_fresh slk p env in + let zp, _ = P.accumulate slk R.m_one np in + P.is_empty zp); + let info, chang1 = set_min_bound info mini in + let info, chang2 = set_max_bound info maxi in + let old_val = info.value in + let info, changed = ajust_value_of_non_basic info in + let env = + { env with non_basic = MX.add slk (info, use) env.non_basic } + in + let env = + if not changed then env + else adapt_values_of_basic_vars env old_val info.value slk use + in + (info, false, env, chang1 || chang2) + with Not_found -> ( + try + (* basic existing var ? *) + let info, poly = MX.find slk env.basic in assert ( let np, _ = normalize_polynomial is_fresh slk p env in - let zp, _ = P.accumulate slk R.m_one np in - P.is_empty zp - ); + P.equal np poly); let info, chang1 = set_min_bound info mini in let info, chang2 = set_max_bound info maxi in - let old_val = info.value in - let info, changed = ajust_value_of_non_basic info in - let env = - {env with non_basic = MX.add slk (info, use) env.non_basic} - in - let env = - if not changed then env - else adapt_values_of_basic_vars env old_val info.value slk use - in - info, false, env, chang1 || chang2 + ( info, + true, + { env with basic = MX.add slk (info, poly) env.basic }, + chang1 || chang2 ) with Not_found -> - try (* basic existing var ? *) - let info, poly = MX.find slk env.basic in - assert ( - let np, _ = normalize_polynomial is_fresh slk p env in - P.equal np poly - ); - let info, chang1 = set_min_bound info mini in - let info, chang2 = set_max_bound info maxi in - info, true, {env with basic = MX.add slk (info, poly) env.basic}, - chang1 || chang2 - - with Not_found -> (* fresh basic var *) - assert (is_fresh); - let np, env = normalize_polynomial is_fresh slk p env in - let info = empty_info (evaluate_poly env np) in - let info, chang1 = set_min_bound info mini in - let info, chang2 = set_max_bound info maxi in - info, true, {env with basic = MX.add slk (info, np) env.basic}, - chang1 || chang2 - in - let status, fixme = - if is_basic then - new_status_basic env.status env.fixme slk info - (consistent_bounds info) - else - new_status_non_basic slk env.status env.fixme info - in - let env = {env with status; fixme } in - debug "[exit of assert_poly]" env (Result.get None); - check_invariants env (Result.get None); - env, change - - - end + (* fresh basic var *) + assert is_fresh; + let np, env = normalize_polynomial is_fresh slk p env in + let info = empty_info (evaluate_poly env np) in + let info, chang1 = set_min_bound info mini in + let info, chang2 = set_max_bound info maxi in + ( info, + true, + { env with basic = MX.add slk (info, np) env.basic }, + chang1 || chang2 )) + in + let status, fixme = + if is_basic then + new_status_basic env.status env.fixme slk info (consistent_bounds info) + else new_status_non_basic slk env.status env.fixme info + in + let env = { env with status; fixme } in + debug "[exit of assert_poly]" env (Result.get None); + check_invariants env (Result.get None); + (env, change) +end (* end of functor Make *) diff --git a/src/assertBounds.mli b/src/assertBounds.mli index 5e7d428..1f1d84a 100644 --- a/src/assertBounds.mli +++ b/src/assertBounds.mli @@ -5,39 +5,35 @@ (******************************************************************************) module type S = sig - module Core : CoreSig.S - (* The returned bool is [true] if the asserted bounds are not trivial - (i.e. not implied by known bounds) *) + (* The returned bool is [true] if the asserted bounds are not trivial + (i.e. not implied by known bounds) *) + + val var : + Core.t -> ?min:Core.bound -> ?max:Core.bound -> Core.Var.t -> Core.t * bool (** [var env min max x] returns a new environment obtained by changing the bounds of [x] in [env] to [min] and [max]. If the bounds were implied by other known bounds (in other words, if the environment did not change) the associated boolean will be [false]. *) - val var : + + (* The returned bool is [true] if the asserted bounds are not trivial + (i.e. not implied by known bounds) *) + + val poly : Core.t -> + Core.P.t -> ?min:Core.bound -> ?max:Core.bound -> Core.Var.t -> Core.t * bool - - (* The returned bool is [true] if the asserted bounds are not trivial - (i.e. not implied by known bounds) *) (** [poly env poly min max x] returns a new environment obtained by changing the bounds of [poly] in [env] to [min] and [max]. The polynomial is represented by the slack variable [x]. If the bounds were implied by other known bounds (in other words, if the environment did not change) the associated boolean will be [false]. *) - val poly : - Core.t -> - Core.P.t -> - ?min:Core.bound -> - ?max:Core.bound -> - Core.Var.t -> - Core.t * bool - end -module Make(Core : CoreSig.S) : S with module Core = Core +module Make (Core : CoreSig.S) : S with module Core = Core diff --git a/src/basic.ml b/src/basic.ml index 3e2535d..aa1acaa 100644 --- a/src/basic.ml +++ b/src/basic.ml @@ -6,26 +6,20 @@ open ExtSigs -module Make - (Var : Variables) - (R : Rationals) - (Ex : Explanations) : sig - - module Core : CoreSig.S - with module Var=Var - and type R.t = R.t - and type V.t = R.t - and module Ex=Ex +module Make (Var : Variables) (R : Rationals) (Ex : Explanations) : sig + module Core : + CoreSig.S + with module Var = Var + and type R.t = R.t + and type V.t = R.t + and module Ex = Ex module Assert : AssertBounds.S with module Core := Core - module Solve : SolveBounds.S with module Core := Core - module Result : Result.S with module Core := Core - + module Solve : SolveBounds.S with module Core := Core + module Result : Result.S with module Core := Core end = struct - - module Core = Core.Make(Var)(R)(Ex) - module Assert = AssertBounds.Make(Core) - module Solve = SolveBounds.Make(Core) - module Result = Result.Make(Core) - + module Core = Core.Make (Var) (R) (Ex) + module Assert = AssertBounds.Make (Core) + module Solve = SolveBounds.Make (Core) + module Result = Result.Make (Core) end diff --git a/src/basic.mli b/src/basic.mli index a6d21f3..e3929b7 100644 --- a/src/basic.mli +++ b/src/basic.mli @@ -9,20 +9,17 @@ open ExtSigs (** The main entry point of the library. It provides a functor building each key module of OcplibSimplex. *) -module Make - (Var : Variables)(R : Rationals)(Ex : Explanations) : sig - +module Make (Var : Variables) (R : Rationals) (Ex : Explanations) : sig (** The core module defines the different data types used by the project and some functions to handle them. *) - module Core : CoreSig.S - with module Var=Var - and type R.t=R.t - and type V.t = R.t - and module Ex=Ex + module Core : + CoreSig.S + with module Var = Var + and type R.t = R.t + and type V.t = R.t + and module Ex = Ex module Assert : AssertBounds.S with module Core := Core - module Solve : SolveBounds.S with module Core := Core - module Result : Result.S with module Core := Core - + module Solve : SolveBounds.S with module Core := Core + module Result : Result.S with module Core := Core end - diff --git a/src/core.ml b/src/core.ml index 1fdda9d..4170b40 100644 --- a/src/core.ml +++ b/src/core.ml @@ -7,48 +7,45 @@ open Format open ExtSigs -let src = Logs.Src.create "OcplibSimplex.Core" ~doc:"The Core of the simplex solver" +let src = + Logs.Src.create "OcplibSimplex.Core" ~doc:"The Core of the simplex solver" module MakeExpert (Var : Variables) - (R : Coefs) - (V : Value with type r = R.t) - (Ex : Explanations) - (R2 : Rat2.SIG with module R = R and module V = V) - (P : Polys.SIG with module Var = Var and module R = R) + (R : Coefs) + (V : Value with type r = R.t) + (Ex : Explanations) + (R2 : Rat2.SIG with module R = R and module V = V) + (P : Polys.SIG with module Var = Var and module R = R) (MX : MapSig with type key = Var.t) - (SX : SetSig with type elt = Var.t) - : CoreSig.S with module Var=Var and module R=R and module V = V - and module Ex=Ex and - module P = P and module MX = MX and module SX = SX = struct - + (SX : SetSig with type elt = Var.t) : + CoreSig.S + with module Var = Var + and module R = R + and module V = V + and module Ex = Ex + and module P = P + and module MX = MX + and module SX = SX = struct module Var = Var - module R = R + module R = R module V = V - module Ex = Ex - - module R2 = R2 - + module Ex = Ex + module R2 = R2 module P = P - module MX = MX module SX = SX - type bound = { - bvalue : R2.t; - explanation : Ex.t - } - + type bound = { bvalue : R2.t; explanation : Ex.t } type value_status = ValueOK | LowerKO | UpperKO - type var_info = - { - mini : bound option; - maxi : bound option; - value : R2.t; - vstatus : value_status; - empty_dom : bool; - } + type var_info = { + mini : bound option; + maxi : bound option; + value : R2.t; + vstatus : value_status; + empty_dom : bool; + } type solution = { main_vars : (Var.t * V.t) list; @@ -57,10 +54,10 @@ module MakeExpert epsilon : V.t; } - type maximum = - { max_v : bound; - is_le : bool; (* bool = true <-> large bound *) - } + type maximum = { + max_v : bound; + is_le : bool; (* bool = true <-> large bound *) + } type result = | Unknown @@ -71,35 +68,33 @@ module MakeExpert type simplex_status = UNK | UNSAT of Var.t | SAT - type t = - { - basic : (var_info * P.t) MX.t; - non_basic : (var_info * SX.t) MX.t; - slake : P.t MX.t; - fixme : SX.t; - is_int : bool; - status : simplex_status; - check_invs: bool; - nb_pivots : int ref; - } + type t = { + basic : (var_info * P.t) MX.t; + non_basic : (var_info * SX.t) MX.t; + slake : P.t MX.t; + fixme : SX.t; + is_int : bool; + status : simplex_status; + check_invs : bool; + nb_pivots : int ref; + } let empty ~is_int ~check_invs = { - basic = MX.empty; + basic = MX.empty; non_basic = MX.empty; - slake = MX.empty; - fixme = SX.empty; - status = UNK; + slake = MX.empty; + fixme = SX.empty; + status = UNK; is_int; check_invs; - nb_pivots = ref 0 + nb_pivots = ref 0; } let on_integers env = env.is_int - let equals_optimum (b : R2.t) (opt : bound option) = match opt with - | None -> false - | Some opt -> R2.compare b opt.bvalue = 0 + let equals_optimum (b : R2.t) (opt : bound option) = + match opt with None -> false | Some opt -> R2.compare b opt.bvalue = 0 let violates_min_bound (b : R2.t) (mn : bound option) = match mn with @@ -114,71 +109,74 @@ module MakeExpert let consistent_bound_value min max = R2.compare min max <= 0 let consistent_bounds_aux (mini : bound option) (maxi : bound option) = - match mini, maxi with + match (mini, maxi) with | None, None | Some _, None | None, Some _ -> true - | Some {bvalue = min; _}, Some {bvalue = max; _} -> consistent_bound_value min max + | Some { bvalue = min; _ }, Some { bvalue = max; _ } -> + consistent_bound_value min max - let consistent_bounds info = - consistent_bounds_aux info.mini info.maxi + let consistent_bounds info = consistent_bounds_aux info.mini info.maxi let set_min_bound info (bnd : bound option) = match bnd with - | None -> info, false + | None -> (info, false) | Some _new -> - let mini = info.mini in - if violates_min_bound _new.bvalue mini || equals_optimum _new.bvalue mini then - info, false - else - let empty_dom = not (consistent_bounds_aux bnd info.maxi) in - let i' = - if violates_min_bound info.value bnd then - {info with mini = bnd; vstatus = LowerKO; empty_dom} - else - {info with mini = bnd; empty_dom} - in i', true + let mini = info.mini in + if + violates_min_bound _new.bvalue mini || equals_optimum _new.bvalue mini + then (info, false) + else + let empty_dom = not (consistent_bounds_aux bnd info.maxi) in + let i' = + if violates_min_bound info.value bnd then + { info with mini = bnd; vstatus = LowerKO; empty_dom } + else { info with mini = bnd; empty_dom } + in + (i', true) let set_max_bound info (bnd : bound option) = match bnd with - | None -> info, false + | None -> (info, false) | Some _new -> - let maxi = info.maxi in - if violates_max_bound _new.bvalue maxi || equals_optimum _new.bvalue maxi then - info, false - else - let empty_dom = not (consistent_bounds_aux info.mini bnd) in - let i' = - if violates_max_bound info.value bnd then - {info with maxi = bnd; vstatus = UpperKO; empty_dom} - else - {info with maxi = bnd; empty_dom} - in - i', true + let maxi = info.maxi in + if + violates_max_bound _new.bvalue maxi || equals_optimum _new.bvalue maxi + then (info, false) + else + let empty_dom = not (consistent_bounds_aux info.mini bnd) in + let i' = + if violates_max_bound info.value bnd then + { info with maxi = bnd; vstatus = UpperKO; empty_dom } + else { info with maxi = bnd; empty_dom } + in + (i', true) let ajust_value_of_non_basic info = - if info.empty_dom then begin + if info.empty_dom then ( assert (info.vstatus != ValueOK); - info, false (* not changed if not sat_bnds *) - end + (info, false (* not changed if not sat_bnds *))) else match info.vstatus with - | ValueOK -> - info, false + | ValueOK -> (info, false) | UpperKO -> - {info with - vstatus = ValueOK; - value = match info.maxi with - | None -> assert false - | Some bnd -> bnd.bvalue}, - true - + ( { + info with + vstatus = ValueOK; + value = + (match info.maxi with + | None -> assert false + | Some bnd -> bnd.bvalue); + }, + true ) | LowerKO -> - {info with - vstatus = ValueOK; - value = match info.mini with - | None -> assert false - | Some bnd -> bnd.bvalue}, - true - + ( { + info with + vstatus = ValueOK; + value = + (match info.mini with + | None -> assert false + | Some bnd -> bnd.bvalue); + }, + true ) let ajust_status_of_basic info = let _new = @@ -186,18 +184,16 @@ module MakeExpert else if violates_max_bound info.value info.maxi then UpperKO else ValueOK in - if info.vstatus == _new then info else {info with vstatus = _new} - + if info.vstatus == _new then info else { info with vstatus = _new } - let evaluate_poly {non_basic; _} p = + let evaluate_poly { non_basic; _ } p = P.fold (fun x c acc -> - let {value = v; _}, _ = - try MX.find x non_basic with Not_found -> assert false - in - R2.add acc (R2.mult_by_const c v) - )p R2.zero - + let { value = v; _ }, _ = + try MX.find x non_basic with Not_found -> assert false + in + R2.add acc (R2.mult_by_const c v)) + p R2.zero let poly_of_slake env slk = try Some (MX.find slk env.slake) with Not_found -> None @@ -205,17 +201,16 @@ module MakeExpert let expl_of_min_bound vinfo = match vinfo.mini with | None -> Ex.empty - | Some {explanation; _} -> explanation + | Some { explanation; _ } -> explanation let expl_of_max_bound vinfo = match vinfo.maxi with | None -> Ex.empty - | Some {explanation; _} -> explanation + | Some { explanation; _ } -> explanation (* debug functions *) module Debug = struct - let string_of_status = function | ValueOK -> "OK " | UpperKO -> "KO(Upper)" @@ -225,11 +220,11 @@ module MakeExpert match i.mini with | None -> fprintf fmt "-∞ <" | Some min -> fprintf fmt "%a <=" R2.print min.bvalue + let print_max_bound fmt i = match i.maxi with | None -> fprintf fmt "< +∞" - | Some max -> - fprintf fmt "<= %a" R2.print max.bvalue + | Some max -> fprintf fmt "<= %a" R2.print max.bvalue let re_computed_status_of_info v = if violates_min_bound v.value v.mini then LowerKO @@ -239,71 +234,56 @@ module MakeExpert let print_bounds_and_values fmt mx = MX.iter (fun x (info, _) -> - let v = info.value in - let comp_status = re_computed_status_of_info info in - Format.fprintf fmt - "%a [ %a == %a ] %a (computed %s) (flag %s)@." - print_min_bound info - Var.print x - R2.print v - print_max_bound info - (string_of_status comp_status) - (string_of_status info.vstatus); - assert (info.empty_dom || comp_status == info.vstatus); - )mx + let v = info.value in + let comp_status = re_computed_status_of_info info in + Format.fprintf fmt + "%a [ %a == %a ] %a (computed %s) (flag %s)@." print_min_bound + info Var.print x R2.print v print_max_bound info + (string_of_status comp_status) + (string_of_status info.vstatus); + assert (info.empty_dom || comp_status == info.vstatus)) + mx let print_uses fmt non_basic = MX.iter (fun x (_, use_x) -> - Format.fprintf fmt - "variables that use %a are:%a@." - Var.print x - (fun fmt s -> SX.iter(fprintf fmt " %a," Var.print) s) use_x; - )non_basic + Format.fprintf fmt "variables that use %a are:%a@." Var.print x + (fun fmt s -> SX.iter (fprintf fmt " %a," Var.print) s) + use_x) + non_basic let print_solution = let aux fmt l = List.iter - (fun (x, q) -> - fprintf fmt " %a --> %a@." Var.print x V.print q; - )l + (fun (x, q) -> fprintf fmt " %a --> %a@." Var.print x V.print q) + l in fun is_int fmt s -> - if is_int - then fprintf fmt " (int solution ? %b)@." - s.int_sol; + if is_int then fprintf fmt " (int solution ? %b)@." s.int_sol; aux fmt s.main_vars; aux fmt s.slake_vars let print_result is_int fmt status = match status with - | Unknown -> fprintf fmt "Unknown" - | Sat s -> - fprintf fmt "Sat:@.%a@." (print_solution is_int) (Lazy.force s) - - | Unsat ex -> - fprintf fmt "Unsat:%a@." Ex.print (Lazy.force ex) - + | Unknown -> fprintf fmt "Unknown" + | Sat s -> fprintf fmt "Sat:@.%a@." (print_solution is_int) (Lazy.force s) + | Unsat ex -> fprintf fmt "Unsat:%a@." Ex.print (Lazy.force ex) | Unbounded s -> - fprintf fmt "Unbounded:@.%a@." (print_solution is_int) (Lazy.force s) - - | Max(mx, s) -> - let mx = Lazy.force mx in - fprintf fmt "Max: (v=%a, is_le=%b, ex=%a)@.%a@." - R2.print mx.max_v.bvalue mx.is_le Ex.print mx.max_v.explanation - (print_solution is_int) (Lazy.force s) + fprintf fmt "Unbounded:@.%a@." (print_solution is_int) (Lazy.force s) + | Max (mx, s) -> + let mx = Lazy.force mx in + fprintf fmt "Max: (v=%a, is_le=%b, ex=%a)@.%a@." R2.print + mx.max_v.bvalue mx.is_le Ex.print mx.max_v.explanation + (print_solution is_int) (Lazy.force s) let print_fixme fmt sx = match SX.elements sx with - | [] -> fprintf fmt " (fixme is empty)@."; + | [] -> fprintf fmt " (fixme is empty)@." | l -> List.iter (fprintf fmt " >> %a@." Var.print) l let print_matrix fmt env = MX.iter - (fun x (_, p) -> - fprintf fmt "%a = %a@." - Var.print x - P.print p) + (fun x (_, p) -> fprintf fmt "%a = %a@." Var.print x P.print p) env.basic let print result fmt env = @@ -328,42 +308,41 @@ module MakeExpert fprintf fmt "--- simplex status --------------------------------@."; fprintf fmt "%a@." (print_result env.is_int) result; fprintf fmt - "== end ==========================================@."; + "== end ==========================================@." end (* end of module Debug *) let print = Debug.print - let debug msg env get_result = + let debug msg env get_result = if Logs.Src.level src <> None then let result = get_result env in Logs.app ~src (fun p -> p "@.%s@.%a@." msg (print result) env) - (* + (* check invariants of the simplex: these invariants are listed in extra/simplexe_invariants.txt *) - module SP : Set.S with type elt = P.t = Set.Make(P) + module SP : Set.S with type elt = P.t = Set.Make (P) let get_all_polys env = - let sp = MX.fold (fun _ (_,p) sp -> SP.add p sp) env.basic SP.empty in + let sp = MX.fold (fun _ (_, p) sp -> SP.add p sp) env.basic SP.empty in MX.fold (fun _ p sp -> SP.add p sp) env.slake sp let get_all_vars env all_polys = let sx = env.fixme in let sx = MX.fold (fun x _ sx -> SX.add x sx) env.basic sx in let sx = - MX.fold (fun x (_, use) sx -> - SX.union use (SX.add x sx)) env.non_basic sx + MX.fold (fun x (_, use) sx -> SX.union use (SX.add x sx)) env.non_basic sx in let sx = MX.fold (fun x _ sx -> SX.add x sx) env.slake sx in SP.fold (P.fold (fun x _ sx -> SX.add x sx)) all_polys sx let info_of x env = try fst (MX.find x env.non_basic) - with Not_found -> - try fst (MX.find x env.basic) with Not_found -> assert false + with Not_found -> ( + try fst (MX.find x env.basic) with Not_found -> assert false) let _01__check_type is_int all_vars = SX.iter (fun x -> assert (is_int == Var.is_int x)) all_vars @@ -374,78 +353,73 @@ module MakeExpert let _03__check_vars_of_polys env polys = SP.iter - (P.iter - (fun x c -> - assert (R.sign c <> 0); - assert (MX.mem x env.basic || MX.mem x env.non_basic) - ))polys + (P.iter (fun x c -> + assert (R.sign c <> 0); + assert (MX.mem x env.basic || MX.mem x env.non_basic))) + polys let _04_05_06__check_use env = MX.iter (fun x (_, use) -> - SX.iter (fun s -> - assert (not (MX.mem s env.non_basic)); (* 04 *) - try assert (P.mem x (snd (MX.find s env.basic))) (*05*) - with Not_found -> assert false - ) use - )env.non_basic; + SX.iter + (fun s -> + assert (not (MX.mem s env.non_basic)); + (* 04 *) + try assert (P.mem x (snd (MX.find s env.basic))) (*05*) + with Not_found -> assert false) + use) + env.non_basic; MX.iter (fun s (_, p) -> - P.iter (fun x _ -> - try assert (SX.mem s (snd (MX.find x env.non_basic))); (*06*) - with Not_found -> assert false - )p; - )env.basic + P.iter + (fun x _ -> + try assert (SX.mem s (snd (MX.find x env.non_basic))) + with (*06*) + | Not_found -> assert false) + p) + env.basic let _07__values_ok_for_non_basic_vars env = MX.iter (fun _ (info, _) -> - if consistent_bounds info then begin - assert (not (violates_min_bound info.value info.mini)); - assert (not (violates_max_bound info.value info.maxi)); - end - else - match env.status with - | UNSAT _ -> () - | SAT | UNK -> assert false - )env.non_basic + if consistent_bounds info then ( + assert (not (violates_min_bound info.value info.mini)); + assert (not (violates_max_bound info.value info.maxi))) + else match env.status with UNSAT _ -> () | SAT | UNK -> assert false) + env.non_basic let _08_09__values_ok_when_sat env result = let check mx int_sol = MX.iter (fun _ (info, _) -> - let v = info.value in - assert (not (violates_min_bound v info.mini)); - assert (not (violates_max_bound v info.maxi)); - assert (not int_sol || V.is_int v.R2.v && R.is_zero v.R2.offset) - ) mx + let v = info.value in + assert (not (violates_min_bound v info.mini)); + assert (not (violates_max_bound v info.maxi)); + assert ((not int_sol) || (V.is_int v.R2.v && R.is_zero v.R2.offset))) + mx in match result with | Unsat _ | Unknown -> () - | Sat s | Unbounded s | Max(_,s) -> - let s = Lazy.force s in - check env.basic s.int_sol; check env.non_basic s.int_sol - + | Sat s | Unbounded s | Max (_, s) -> + let s = Lazy.force s in + check env.basic s.int_sol; + check env.non_basic s.int_sol let _10_11__check_handling_strict_ineqs env = let is_int = env.is_int in let aux _ (info, _) = - begin - match info.mini with - | None -> () - | Some m -> + (match info.mini with + | None -> () + | Some m -> let i = m.bvalue.R2.offset in - assert (not is_int || R.is_zero i); - assert (is_int || R.is_zero i || R.is_one i); - end; - begin - match info.maxi with - | None -> () - | Some m -> + assert ((not is_int) || R.is_zero i); + assert (is_int || R.is_zero i || R.is_one i)); + match info.maxi with + | None -> () + | Some m -> let i = m.bvalue.R2.offset in - assert (not is_int || R.is_zero i); - assert (is_int || R.is_zero i || R.is_m_one i); - end + assert ((not is_int) || R.is_zero i); + assert (is_int || R.is_zero i || R.is_m_one i) in MX.iter aux env.basic; MX.iter aux env.non_basic @@ -454,54 +428,55 @@ module MakeExpert let aux l env = List.iter (fun (x, v) -> - let info = info_of x env in - let v2 = R2.of_r v in - assert (not (violates_min_bound v2 info.mini)); - assert (not (violates_max_bound v2 info.maxi)); - )l - in fun env result -> + let info = info_of x env in + let v2 = R2.of_r v in + assert (not (violates_min_bound v2 info.mini)); + assert (not (violates_max_bound v2 info.maxi))) + l + in + fun env result -> match result with | Unsat _ | Unknown -> () - | Sat s | Unbounded s | Max(_,s) -> - let s = Lazy.force s in - let v = List.length s.main_vars + List.length s.slake_vars in - let w = MX.cardinal env.non_basic + MX.cardinal env.basic in - assert ( - if v <> w then - eprintf "model length = %d, but basic + non_basic = %d@." v w; - v = w); - aux s.main_vars env; - aux s.slake_vars env + | Sat s | Unbounded s | Max (_, s) -> + let s = Lazy.force s in + let v = List.length s.main_vars + List.length s.slake_vars in + let w = MX.cardinal env.non_basic + MX.cardinal env.basic in + assert ( + if v <> w then + eprintf "model length = %d, but basic + non_basic = %d@." v w; + v = w); + aux s.main_vars env; + aux s.slake_vars env let _13__check_reason_when_unsat _env = - Logs.app ~src (fun p -> - p ~header:"[check-invariants]" "@._13__check_reason_when_unsat: TODO@.@." - ) + Logs.app ~src (fun p -> + p ~header:"[check-invariants]" + "@._13__check_reason_when_unsat: TODO@.@.") let _14_15__fixme_is_subset_of_basic env = SX.iter (fun x -> - try - let info, _ = MX.find x env.basic in - assert - ((violates_min_bound info.value info.mini) || - (violates_max_bound info.value info.maxi)); (*15*) - with Not_found -> assert false (*14*) - - ) env.fixme - - let _16__fixme_containts_basic_with_bad_values_if_not_unsat - env all_vars result = + try + let info, _ = MX.find x env.basic in + assert ( + violates_min_bound info.value info.mini + || violates_max_bound info.value info.maxi) + (*15*) + with Not_found -> assert false (*14*)) + env.fixme + + let _16__fixme_containts_basic_with_bad_values_if_not_unsat env all_vars + result = match result with | Unsat _ | Unbounded _ | Max _ -> () | Unknown | Sat _ -> - SX.iter - (fun x -> - if not (SX.mem x env.fixme) then - let info = info_of x env in - assert (not (violates_min_bound info.value info.mini)); - assert (not (violates_max_bound info.value info.maxi)) - )all_vars + SX.iter + (fun x -> + if not (SX.mem x env.fixme) then ( + let info = info_of x env in + assert (not (violates_min_bound info.value info.mini)); + assert (not (violates_max_bound info.value info.maxi)))) + all_vars let _17__fixme_is_empty_if_not_unknown env result = match result with @@ -510,22 +485,27 @@ module MakeExpert let _18__vals_of_basic_vars_computation env = MX.iter - (fun _ ({value = s; _}, p) -> - let vp = evaluate_poly env p in - assert (R2.equal vp s); - )env.basic + (fun _ ({ value = s; _ }, p) -> + let vp = evaluate_poly env p in + assert (R2.equal vp s)) + env.basic let _19__check_that_vstatus_are_well_set env = let aux _ (info, _) = - if info.empty_dom then - assert (info.vstatus != ValueOK) + if info.empty_dom then assert (info.vstatus != ValueOK) else let vmin = violates_min_bound info.value info.mini in let vmax = violates_max_bound info.value info.maxi in match info.vstatus with - | ValueOK -> assert (not vmin); assert(not vmax); - | UpperKO -> assert (not vmin); assert(vmax); - | LowerKO -> assert (vmin); assert(not vmax); + | ValueOK -> + assert (not vmin); + assert (not vmax) + | UpperKO -> + assert (not vmin); + assert vmax + | LowerKO -> + assert vmin; + assert (not vmax) in MX.iter aux env.basic; MX.iter aux env.non_basic @@ -534,26 +514,27 @@ module MakeExpert match result with | Unsat _ -> () | Unknown | Sat _ | Unbounded _ | Max _ -> - let aux _ (info, _) = assert (consistent_bounds info) in - MX.iter aux env.basic; - MX.iter aux env.non_basic + let aux _ (info, _) = assert (consistent_bounds info) in + MX.iter aux env.basic; + MX.iter aux env.non_basic let _21__check_coherence_of_empty_dom = let aux mx = MX.iter (fun _ (info, _) -> - assert (consistent_bounds info == not info.empty_dom); - if info.empty_dom then - assert (violates_min_bound info.value info.mini || - violates_max_bound info.value info.maxi); - )mx + assert (consistent_bounds info == not info.empty_dom); + if info.empty_dom then + assert ( + violates_min_bound info.value info.mini + || violates_max_bound info.value info.maxi)) + mx in fun env -> aux env.non_basic; aux env.basic let check_invariants env get_result = - if env.check_invs then + if env.check_invs then ( let polys = get_all_polys env in let all_vars = get_all_vars env polys in let result = get_result env in @@ -567,28 +548,32 @@ module MakeExpert _12__check_solution_when_sat env result; _13__check_reason_when_unsat env; _14_15__fixme_is_subset_of_basic env; - _16__fixme_containts_basic_with_bad_values_if_not_unsat - env all_vars result; + _16__fixme_containts_basic_with_bad_values_if_not_unsat env all_vars + result; _17__fixme_is_empty_if_not_unknown env result; _18__vals_of_basic_vars_computation env; _19__check_that_vstatus_are_well_set env; - _20__bounds_are_consistent_if_not_unsat env result; - (*_21__check_coherence_of_empty_dom env;*) - - + _20__bounds_are_consistent_if_not_unsat env result) + (*_21__check_coherence_of_empty_dom env;*) end -module Make - (Var : Variables) - (R : Rationals) - (Ex : Explanations) - : CoreSig.S with module Var=Var and type R.t = R.t and type V.t = R.t and module Ex=Ex = struct +module Make (Var : Variables) (R : Rationals) (Ex : Explanations) : + CoreSig.S + with module Var = Var + and type R.t = R.t + and type V.t = R.t + and module Ex = Ex = struct module V' = struct include R + type r = t + let mult_by_coef = mult let div_by_coef = div end - include MakeExpert(Var)(R)(V')(Ex)(Rat2.Make(R)(V'))(Polys.Make(Var)(R)) - (Map.Make(Var))(Set.Make(Var)) + + include + MakeExpert (Var) (R) (V') (Ex) (Rat2.Make (R) (V')) (Polys.Make (Var) (R)) + (Map.Make (Var)) + (Set.Make (Var)) end diff --git a/src/core.mli b/src/core.mli index feb337a..6409e4a 100644 --- a/src/core.mli +++ b/src/core.mli @@ -8,26 +8,29 @@ open ExtSigs val src : Logs.src -module Make - (Var : Variables) - (R : Rationals) - (Ex : Explanations) - : CoreSig.S +module Make (Var : Variables) (R : Rationals) (Ex : Explanations) : + CoreSig.S with module Var = Var and type R.t = R.t and type V.t = R.t and module Ex = Ex +(** Same than Make but allows to choose the implementation of polynomials, maps + and sets *) module MakeExpert (Var : Variables) - (R : Coefs) - (V : Value with type r = R.t) - (Ex : Explanations) - (R2 : Rat2.SIG with module R = R and module V = V) - (P : Polys.SIG with module Var = Var and module R = R) + (R : Coefs) + (V : Value with type r = R.t) + (Ex : Explanations) + (R2 : Rat2.SIG with module R = R and module V = V) + (P : Polys.SIG with module Var = Var and module R = R) (MX : MapSig with type key = Var.t) - (SX : SetSig with type elt = Var.t) - : CoreSig.S with module Var=Var and module R=R and module V=V and module Ex=Ex and - module P = P and module MX = MX and module SX = SX -(** Same than Make but allows to choose the implementation of polynomials, maps - and sets *) + (SX : SetSig with type elt = Var.t) : + CoreSig.S + with module Var = Var + and module R = R + and module V = V + and module Ex = Ex + and module P = P + and module MX = MX + and module SX = SX diff --git a/src/coreSig.mli b/src/coreSig.mli index b935f99..c8afa0b 100644 --- a/src/coreSig.mli +++ b/src/coreSig.mli @@ -8,58 +8,55 @@ open ExtSigs (** Interface of the main types and auxiliary of the simplex. *) module type S = sig - (** {1 Modules} *) - (** The type of variables maipulated by the simplex algorithm. *) module Var : Variables + (** The type of variables maipulated by the simplex algorithm. *) + module Ex : Explanations (** An interface for explanations; in practice, they are labels attached to bounds used for backtracking information on how bounds were discovered. The simplex algorithm does not create explanations: it will only attach empty explanations to bounds, build the union of explanations and print them. It is the user's job to provide the initial explanations when initializing the simplex core. *) - module Ex : Explanations + module R : Coefs (** The interface for rationals provided by users for the coefficient. *) - module R : Coefs + module V : Value with type r = R.t (** The interface for values of variable and bounds provided by users. *) - module V : Value with type r = R.t (** Pairs of rationals R representing bounds with an offset [x + kƐ]. *) - module R2 : Rat2.SIG with module R = R and module V = V + module R2 : Rat2.SIG with module R = R and module V = V (** Linear relations of variables. *) module P : Polys.SIG with module Var = Var and module R = R - (** Collections of variables. *) module MX : MapSig with type key = Var.t + (** Collections of variables. *) + module SX : SetSig with type elt = Var.t (*module SLAKE : Map.S with type key = P.t*) (** {1 Types} *) - + + type bound = { bvalue : R2.t; explanation : Ex.t } (** A bound is a value of the form [x + kƐ] and an explanation. *) - type bound = { - bvalue : R2.t; - explanation : Ex.t - } type value_status = - | ValueOK (** The value is inbetween bounds. *) - | LowerKO (** The value is smaller than the lower bound. *) - | UpperKO (** The value is greater than the upper bound. *) + | ValueOK (** The value is inbetween bounds. *) + | LowerKO (** The value is smaller than the lower bound. *) + | UpperKO (** The value is greater than the upper bound. *) type var_info = { - mini : bound option; (* None -> -inf *) - maxi : bound option; (* None -> +inf *) - value : R2.t; - vstatus : value_status; - empty_dom : bool; - } + mini : bound option; (* None -> -inf *) + maxi : bound option; (* None -> +inf *) + value : R2.t; + vstatus : value_status; + empty_dom : bool; + } type solution = { main_vars : (Var.t * V.t) list; @@ -83,88 +80,90 @@ module type S = sig type simplex_status = UNK | UNSAT of Var.t | SAT type t = { - basic : (var_info * P.t) MX.t; - non_basic : (var_info * SX.t) MX.t; - slake : P.t MX.t; - fixme : SX.t; - is_int : bool; - status : simplex_status; - check_invs : bool; - nb_pivots : int ref; - } + basic : (var_info * P.t) MX.t; + non_basic : (var_info * SX.t) MX.t; + slake : P.t MX.t; + fixme : SX.t; + is_int : bool; + status : simplex_status; + check_invs : bool; + nb_pivots : int ref; + } + val empty : is_int:bool -> check_invs:bool -> t (** Returns a simplex environment with three parameters: - [is_int]: will the simplex work on an integer optimization problem or a rational problem? - [check_invs]: processes checks after the calculation (deprecated). *) - val empty : is_int : bool -> check_invs : bool -> t - (** Returns [true] if the simplex environment is on integers. *) val on_integers : t -> bool + (** Returns [true] if the simplex environment is on integers. *) - (** Equality check between bounds. *) val equals_optimum : R2.t -> bound option -> bool + (** Equality check between bounds. *) + val consistent_bounds : var_info -> bool (** Checks if the lower bound of a variable is smaller than its upper bound. *) - val consistent_bounds : var_info -> bool - (** [violates_min_bound b mb] returns [true] if [b] is smaller than [mb]. *) val violates_min_bound : R2.t -> bound option -> bool + (** [violates_min_bound b mb] returns [true] if [b] is smaller than [mb]. *) - (** [violates_max_bound b mb] returns [true] if [b] is greater than [mb]. *) val violates_max_bound : R2.t -> bound option -> bool + (** [violates_max_bound b mb] returns [true] if [b] is greater than [mb]. *) (* The returned bool is [true] if the asserted bounds are not trivial (i.e. not implied by known bounds). *) + + val set_min_bound : var_info -> bound option -> var_info * bool (** [set_min_bound vinfo b] returns a couple [(vinfo', changed)] where: - [vinfo'] is the new variable info where the new min bound [b] has been set. - [changed] is [true] if the new bound has changed the variable info *) - val set_min_bound : var_info -> bound option -> var_info * bool - (** Same as {!val:set_min_bound}, but for max bounds. *) val set_max_bound : var_info -> bound option -> var_info * bool + (** Same as {!val:set_min_bound}, but for max bounds. *) + val ajust_value_of_non_basic : var_info -> var_info * bool (** [ajust_value_of_non_basic vinfo] updates the info's value with the upper bound (resp. the lower bound), if [vinfo]'s status is {!const:UpperKO} (resp. {!const:LowerKO}). Otherwise, do nothing. The boolean returned is [true] if the new variable [var_info] has changed. *) - val ajust_value_of_non_basic: var_info -> var_info * bool (* vstatus is supposed to be well set *) + val ajust_status_of_basic : var_info -> var_info (** [ajust_status_of_basic vinfo] checks a variable info's bound matches its status. If its value violates its lower bound, its status is set to {!const:LowerKO}. In the other case, it is set to {!const:UpperKO}. If the value is between the two bounds, it is set to {!const:ValueOK}. *) - val ajust_status_of_basic : var_info -> var_info (* valuation is supposed to be well computed *) - (** Evaluates a polynomial of non basic variables. *) val evaluate_poly : t -> P.t -> R2.t + (** Evaluates a polynomial of non basic variables. *) + val poly_of_slake : t -> Var.t -> P.t option (** [poly_of_slake env slake] returns the polynomial associated to the variable [slake] in [env]. *) - val poly_of_slake : t -> Var.t -> P.t option - (** Returns the explanation associated to a variable lower bound. *) val expl_of_min_bound : var_info -> Ex.t + (** Returns the explanation associated to a variable lower bound. *) - (** Same as `expl_of_min_bound`, but for upper bounds. *) val expl_of_max_bound : var_info -> Ex.t + (** Same as `expl_of_min_bound`, but for upper bounds. *) (** {1 Debug functions} *) + (** Only use for internal debugging *) - (** Checks several invariants in the project *) val check_invariants : t -> (t -> result) -> unit + (** Checks several invariants in the project *) - (** Pretty prints the environment. *) val print : result -> Format.formatter -> t -> unit + (** Pretty prints the environment. *) - (** @deprecated *) val debug : string -> t -> (t -> result) -> unit + (** @deprecated *) end diff --git a/src/dune b/src/dune index 4b8a817..c467006 100644 --- a/src/dune +++ b/src/dune @@ -2,5 +2,4 @@ (name OcplibSimplex) (public_name ocplib-simplex) (libraries logs) - (modules_without_implementation coreSig extSigs) -) + (modules_without_implementation coreSig extSigs)) diff --git a/src/extSigs.mli b/src/extSigs.mli index 0bbba21..81819c6 100644 --- a/src/extSigs.mli +++ b/src/extSigs.mli @@ -4,34 +4,31 @@ (* Copyright (C) --- OCamlPro --- See README.md for information and licensing *) (******************************************************************************) - (*----------------------------------------------------------------------------*) (** Interface required for variables *) module type Variables = sig - - (** type of variables used in the simplex *) type t + (** type of variables used in the simplex *) - (** compare function on vars *) val compare : t -> t -> int + (** compare function on vars *) + val is_int : t -> bool (** [is_int v] returns true if the variable has integer type, and false otherwise *) - val is_int : t -> bool - (** [print fmt v] prints the given var *) val print : Format.formatter -> t -> unit - + (** [print fmt v] prints the given var *) end (*----------------------------------------------------------------------------*) (** Interface required for rationnals *) module type Rationals = sig - - (** type of rationnal numbers *) type t + (** type of rationnal numbers *) + val zero : t val one : t val m_one : t @@ -57,9 +54,9 @@ end (** Interface required for coefs *) module type Coefs = sig - - (** type of rationnal numbers *) type t + (** type of rationnal numbers *) + val zero : t val one : t val m_one : t @@ -83,9 +80,9 @@ end (** Interface required for bounds and solutions *) module type Value = sig - - (** type of rationnal numbers *) type t + (** type of rationnal numbers *) + val zero : t val one : t val m_one : t @@ -104,8 +101,9 @@ module type Value = sig val ceiling : t -> t type r - val mult_by_coef: t -> r -> t - val div_by_coef: t -> r -> t + + val mult_by_coef : t -> r -> t + val div_by_coef : t -> r -> t end (*----------------------------------------------------------------------------*) @@ -113,6 +111,7 @@ end (** Interface of explanations *) module type Explanations = sig type t + val empty : t val union : t -> t -> t val print : Format.formatter -> t -> unit @@ -126,7 +125,7 @@ module type MapSig = sig val find : key -> 'a t -> 'a val add : key -> 'a -> 'a t -> 'a t val remove : key -> 'a t -> 'a t - val mem: key -> 'a t -> bool + val mem : key -> 'a t -> bool val fold : (key -> 'a -> 'b -> 'b) -> 'a t -> 'b -> 'b val iter : (key -> 'a -> unit) -> 'a t -> unit val cardinal : 'a t -> int diff --git a/src/polys.ml b/src/polys.ml index 4803f83..5a47027 100644 --- a/src/polys.ml +++ b/src/polys.ml @@ -16,17 +16,15 @@ module type SIG = sig val empty : t val is_polynomial : t -> bool val is_empty : t -> bool - - val replace : Var.t -> R.t -> t -> t * var_status + val replace : Var.t -> R.t -> t -> t * var_status val accumulate : Var.t -> R.t -> t -> t * var_status - val append : t -> R.t -> t -> t * (Var.t * var_status) list + val append : t -> R.t -> t -> t * (Var.t * var_status) list val subst : Var.t -> t -> t -> t * (Var.t * var_status) list val from_list : (Var.t * R.t) list -> t val print : Format.formatter -> t -> unit - - val fold: (Var.t -> R.t -> 'a -> 'a) -> t -> 'a -> 'a - val iter: (Var.t -> R.t -> unit) -> t -> unit - val partition: (Var.t -> R.t -> bool) -> t -> t * t + val fold : (Var.t -> R.t -> 'a -> 'a) -> t -> 'a -> 'a + val iter : (Var.t -> R.t -> unit) -> t -> unit + val partition : (Var.t -> R.t -> bool) -> t -> t * t val compare : t -> t -> int val mem : Var.t -> t -> bool val equal : t -> t -> bool @@ -35,77 +33,78 @@ module type SIG = sig val remove : Var.t -> t -> t end -module Make(Var: Variables)(R : Rationals) : SIG - with module Var = Var and module R = R = struct - - module Var = Var - module R = R - - module MV = Map.Make(Var) - - type t = R.t MV.t - - type var_status = New | Exists | Removed - - let empty = MV.empty - let fold = MV.fold - let iter = MV.iter - let compare = MV.compare R.compare - let partition = MV.partition - - let remove = MV.remove - let find = MV.find - let bindings = MV.bindings - let equal = MV.equal R.equal - let mem = MV.mem - let is_empty = MV.is_empty - - let is_polynomial p = - try - let cpt = ref 0 in - iter (fun _ _ -> incr cpt; if !cpt > 1 then raise Exit) p; - false - with Exit -> - true +module Make (Var : Variables) (R : Rationals) : + SIG with module Var = Var and module R = R = struct + module Var = Var + module R = R + module MV = Map.Make (Var) - let replace v q t = - if R.is_zero q then MV.remove v t, Removed - else MV.add v q t, (if MV.mem v t then Exists else New) - - let accumulate v q t = - let new_q = try R.add q (find v t) with Not_found -> q in - replace v new_q t + type t = R.t MV.t + type var_status = New | Exists | Removed - (* TODO: We can maybe replace mp with a list, since keys are unique ... *) - let append_aux p coef q = - fold (fun x c (p, mp) -> + let empty = MV.empty + let fold = MV.fold + let iter = MV.iter + let compare = MV.compare R.compare + let partition = MV.partition + let remove = MV.remove + let find = MV.find + let bindings = MV.bindings + let equal = MV.equal R.equal + let mem = MV.mem + let is_empty = MV.is_empty + + let is_polynomial p = + try + let cpt = ref 0 in + iter + (fun _ _ -> + incr cpt; + if !cpt > 1 then raise Exit) + p; + false + with Exit -> true + + let replace v q t = + if R.is_zero q then (MV.remove v t, Removed) + else (MV.add v q t, if MV.mem v t then Exists else New) + + let accumulate v q t = + let new_q = try R.add q (find v t) with Not_found -> q in + replace v new_q t + + (* TODO: We can maybe replace mp with a list, since keys are unique ... *) + let append_aux p coef q = + fold + (fun x c (p, mp) -> let p, x_status = accumulate x (R.mult coef c) p in - p, MV.add x x_status mp - ) q (p, MV.empty) - - let append p coef q = - let p, mp = append_aux p coef q in p, MV.bindings mp - - let subst v p q = - try - let new_q, modified = append_aux (remove v q) (find v q) p in - new_q, MV.bindings (MV.add v Removed modified) - with Not_found -> - (* This will oblige us to enforce strong invariants !! - We should know exactly where we have to substitute !! *) - assert false - - let from_list l = - List.fold_left (fun p (x, c) -> fst (accumulate x c p)) empty l - - let print fmt p = - let l = MV.bindings p in - match l with - | [] -> Format.fprintf fmt "(empty-poly)" - | (x, q)::l -> + (p, MV.add x x_status mp)) + q (p, MV.empty) + + let append p coef q = + let p, mp = append_aux p coef q in + (p, MV.bindings mp) + + let subst v p q = + try + let new_q, modified = append_aux (remove v q) (find v q) p in + (new_q, MV.bindings (MV.add v Removed modified)) + with Not_found -> + (* This will oblige us to enforce strong invariants !! + We should know exactly where we have to substitute !! *) + assert false + + let from_list l = + List.fold_left (fun p (x, c) -> fst (accumulate x c p)) empty l + + let print fmt p = + let l = MV.bindings p in + match l with + | [] -> Format.fprintf fmt "(empty-poly)" + | (x, q) :: l -> Format.fprintf fmt "(%a) * %a" R.print q Var.print x; List.iter - (fun (x,q) -> - Format.fprintf fmt " + (%a) * %a" R.print q Var.print x) l - - end + (fun (x, q) -> + Format.fprintf fmt " + (%a) * %a" R.print q Var.print x) + l +end diff --git a/src/polys.mli b/src/polys.mli index 17df2b0..077b5dc 100644 --- a/src/polys.mli +++ b/src/polys.mli @@ -16,17 +16,15 @@ module type SIG = sig val empty : t val is_polynomial : t -> bool val is_empty : t -> bool - - val replace : Var.t -> R.t -> t -> t * var_status + val replace : Var.t -> R.t -> t -> t * var_status val accumulate : Var.t -> R.t -> t -> t * var_status - val append : t -> R.t -> t -> t * (Var.t * var_status) list + val append : t -> R.t -> t -> t * (Var.t * var_status) list val subst : Var.t -> t -> t -> t * (Var.t * var_status) list val from_list : (Var.t * R.t) list -> t val print : Format.formatter -> t -> unit - - val fold: (Var.t -> R.t -> 'a -> 'a) -> t -> 'a -> 'a - val iter: (Var.t -> R.t -> unit) -> t -> unit - val partition: (Var.t -> R.t -> bool) -> t -> t * t + val fold : (Var.t -> R.t -> 'a -> 'a) -> t -> 'a -> 'a + val iter : (Var.t -> R.t -> unit) -> t -> unit + val partition : (Var.t -> R.t -> bool) -> t -> t * t val compare : t -> t -> int val mem : Var.t -> t -> bool val equal : t -> t -> bool @@ -35,6 +33,5 @@ module type SIG = sig val remove : Var.t -> t -> t end -module Make(Var: Variables)(R : Rationals) : SIG - with module Var = Var and module R = R - +module Make (Var : Variables) (R : Rationals) : + SIG with module Var = Var and module R = R diff --git a/src/rat2.ml b/src/rat2.ml index 805b126..05dea59 100644 --- a/src/rat2.ml +++ b/src/rat2.ml @@ -10,58 +10,49 @@ module type SIG = sig module R : Coefs module V : Value with type r = R.t - type t = private {v: V.t; offset: R.t} + type t = private { v : V.t; offset : R.t } val zero : t val of_r : V.t -> t val upper : V.t -> t val lower : V.t -> t - val minus : t -> t val add : t -> t -> t val sub : t -> t -> t val mult_by_const : R.t -> t -> t val div_by_const : R.t -> t -> t - val compare : t -> t -> int val equal : t -> t -> bool val is_zero : t -> bool val is_pure_rational : t -> bool val is_int : t -> bool - val floor : t -> t val ceiling : t -> t - val print : Format.formatter -> t -> unit end -module Make(R : Rationals)(V : Value with type r = R.t) - : SIG with module R = R and module V = V = struct - +module Make (R : Rationals) (V : Value with type r = R.t) : + SIG with module R = R and module V = V = struct module R = R module V = V type t = { - v: V.t; - - offset: R.t; - (* When working on strict bounds, an epsilon is added to the bounds. - The offset represents the amount of epsilon are added. *) + v : V.t; + offset : R.t; + (* When working on strict bounds, an epsilon is added to the bounds. + The offset represents the amount of epsilon are added. *) } - let of_r v = {v; offset = R.zero} - let upper v = {v; offset = R.m_one} - let lower v = {v; offset = R.one} + let of_r v = { v; offset = R.zero } + let upper v = { v; offset = R.m_one } + let lower v = { v; offset = R.one } let zero = of_r V.zero - let is_pure_rational r = R.equal r.offset R.zero let is_int r = is_pure_rational r && V.is_int r.v - - let map f g a = {v = f a.v; offset = g a.offset} - let map2 f g a b = {v = f a.v b.v; offset = g a.offset b.offset} - - let add = map2 V.add R.add - let sub = map2 V.sub R.sub + let map f g a = { v = f a.v; offset = g a.offset } + let map2 f g a b = { v = f a.v b.v; offset = g a.offset b.offset } + let add = map2 V.add R.add + let sub = map2 V.sub R.sub let mult_by_const c e = if R.is_one c then e @@ -76,35 +67,24 @@ module Make(R : Rationals)(V : Value with type r = R.t) if c <> 0 then c else R.compare a.offset b.offset let equal a b = compare a b = 0 - let is_zero a = V.is_zero a.v && R.is_zero a.offset - let minus = map V.minus R.minus let floor r = - if V.is_int r.v - then - if is_pure_rational r - then r - else of_r (V.sub r.v V.one) + if V.is_int r.v then + if is_pure_rational r then r else of_r (V.sub r.v V.one) else of_r (V.floor r.v) let ceiling r = - if V.is_int r.v - then - if is_pure_rational r - then r - else of_r (V.add r.v V.one) + if V.is_int r.v then + if is_pure_rational r then r else of_r (V.add r.v V.one) else of_r (V.ceiling r.v) let print_offset fmt off = let c = R.compare off R.zero in if c = 0 then () - else if c > 0 - then Format.fprintf fmt "+%aƐ" R.print off + else if c > 0 then Format.fprintf fmt "+%aƐ" R.print off else Format.fprintf fmt "%aƐ" R.print off - let print fmt r = Format.fprintf fmt "%a%a" - V.print r.v - print_offset r.offset + let print fmt r = Format.fprintf fmt "%a%a" V.print r.v print_offset r.offset end diff --git a/src/rat2.mli b/src/rat2.mli index 1e4a3bd..dc2c6f2 100644 --- a/src/rat2.mli +++ b/src/rat2.mli @@ -22,72 +22,70 @@ module type SIG = sig (** {1 Type} *) type t = private { - v: V.t; - (** The raw value of the bound. *) - - offset: R.t - (** The number of epsilons to add to the bound. *) + v : V.t; (** The raw value of the bound. *) + offset : R.t; (** The number of epsilons to add to the bound. *) } (** {1 Constructors} *) - (** The zero bound with no offset. *) val zero : t + (** The zero bound with no offset. *) - (** From a rational [r], returns the Rat2 representation with no offset. *) val of_r : V.t -> t + (** From a rational [r], returns the Rat2 representation with no offset. *) - (** From a rational [r], returns [r - Ɛ]. *) val upper : V.t -> t + (** From a rational [r], returns [r - Ɛ]. *) - (** From a rational [r], returns [r + Ɛ]. *) val lower : V.t -> t + (** From a rational [r], returns [r + Ɛ]. *) (** {1 Algebraic operations} *) - (** From a bound [r + kƐ], returns [-r -kƐ]. *) val minus : t -> t + (** From a bound [r + kƐ], returns [-r -kƐ]. *) - (** Adds two bounds. *) val add : t -> t -> t + (** Adds two bounds. *) - (** Substracts two bounds. *) val sub : t -> t -> t + (** Substracts two bounds. *) + val mult_by_const : R.t -> t -> t (** Multiplies a bound by a rational constant (both v and offset are multiplied). *) - val mult_by_const : R.t -> t -> t - (** Divides a bound by a constant. Fails if the constant is zero. *) val div_by_const : R.t -> t -> t + (** Divides a bound by a constant. Fails if the constant is zero. *) (** {1 Comparison functions} *) - (** Compares two bounds; returns 0 iff the two bounds are strictly equal. *) val compare : t -> t -> int + (** Compares two bounds; returns 0 iff the two bounds are strictly equal. *) - (** Returns [true] iff the bounds are strictly equal. *) val equal : t -> t -> bool + (** Returns [true] iff the bounds are strictly equal. *) - (** Returns [true] iff the bound in argument is zero. *) val is_zero : t -> bool + (** Returns [true] iff the bound in argument is zero. *) - (** Returns [true] iff the offset is 0. *) val is_pure_rational : t -> bool + (** Returns [true] iff the offset is 0. *) - (** Returns [true] iff the offset is 0 and the field [v] is an integer. *) val is_int : t -> bool + (** Returns [true] iff the offset is 0 and the field [v] is an integer. *) (** {1 Misc} *) - (** Returns the greatest (pure) integer smaller or equal to the argument. *) val floor : t -> t + (** Returns the greatest (pure) integer smaller or equal to the argument. *) - (** Returns the smallest (pure) integer greater or equal to the argument. *) val ceiling : t -> t + (** Returns the smallest (pure) integer greater or equal to the argument. *) - (** Prints a bound. *) val print : Format.formatter -> t -> unit + (** Prints a bound. *) end -module Make(R : Rationals)(V : Value with type r = R.t) : SIG with module R = R and module V = V +module Make (R : Rationals) (V : Value with type r = R.t) : + SIG with module R = R and module V = V diff --git a/src/result.ml b/src/result.ml index 7c0f11c..7f4f47c 100644 --- a/src/result.ml +++ b/src/result.ml @@ -6,11 +6,11 @@ module type S = sig module Core : CoreSig.S + val get : (Core.P.t * bool) option -> Core.t -> Core.result end -module Make(Core : CoreSig.S) : S with module Core = Core = struct - +module Make (Core : CoreSig.S) : S with module Core = Core = struct module Core = Core open Core @@ -18,165 +18,150 @@ module Make(Core : CoreSig.S) : S with module Core = Core = struct let invert_sign = if is_lower then 1 else -1 in P.fold (fun x coef ex -> - let c = invert_sign * R.sign coef in - assert (c <> 0); - let xi, _ = try MX.find x non_basic with Not_found -> assert false in - let ex, optimum = - if c > 0 then Ex.union ex (Core.expl_of_max_bound xi), xi.maxi - else Ex.union ex (Core.expl_of_min_bound xi), xi.mini - in - assert (equals_optimum xi.value optimum); - ex - )p ex + let c = invert_sign * R.sign coef in + assert (c <> 0); + let xi, _ = try MX.find x non_basic with Not_found -> assert false in + let ex, optimum = + if c > 0 then (Ex.union ex (Core.expl_of_max_bound xi), xi.maxi) + else (Ex.union ex (Core.expl_of_min_bound xi), xi.mini) + in + assert (equals_optimum xi.value optimum); + ex) + p ex let get_unsat_core s env = try let si, p = MX.find s env.basic in - if not (consistent_bounds si) - then Ex.union - (Core.expl_of_min_bound si) - (Core.expl_of_max_bound si) - else match si.vstatus with + if not (consistent_bounds si) then + Ex.union (Core.expl_of_min_bound si) (Core.expl_of_max_bound si) + else + match si.vstatus with | ValueOK -> assert false - | LowerKO -> explain_poly env.non_basic p (Core.expl_of_min_bound si) true - | UpperKO -> explain_poly env.non_basic p (Core.expl_of_max_bound si) false + | LowerKO -> + explain_poly env.non_basic p (Core.expl_of_min_bound si) true + | UpperKO -> + explain_poly env.non_basic p (Core.expl_of_max_bound si) false with Not_found -> let si, _ = MX.find s env.non_basic in assert (not (consistent_bounds si)); - Ex.union - (Core.expl_of_min_bound si) - (Core.expl_of_max_bound si) + Ex.union (Core.expl_of_min_bound si) (Core.expl_of_max_bound si) let get_int_solution env = let is_int_sol = ref true in let sol = - MX.fold (fun x (xi,_) sol -> + MX.fold + (fun x (xi, _) sol -> let v = xi.value in assert (R2.is_pure_rational v); is_int_sol := !is_int_sol && V.is_int v.R2.v; - (x, v.R2.v) :: sol - ) env.non_basic [] + (x, v.R2.v) :: sol) + env.non_basic [] in let sol = - MX.fold (fun x (xi, _) sol -> + MX.fold + (fun x (xi, _) sol -> let v = xi.value in assert (R2.is_pure_rational v); is_int_sol := !is_int_sol && V.is_int v.R2.v; - (x, v.R2.v) :: sol - )env.basic sol + (x, v.R2.v) :: sol) + env.basic sol in let slake = env.slake in let sol_slk, sol = List.partition (fun (x, _) -> MX.mem x slake) sol in - { main_vars = sol; slake_vars = sol_slk; int_sol = !is_int_sol; epsilon = V.zero} - - - - let eval_eps - (eps : V.t) - (inf : R2.t) - (sup : R2.t) = - let {R2.v = inf_r; offset = inf_d} : R2.t = inf in - let {R2.v = sup_r; offset = sup_d} : R2.t = sup in + { + main_vars = sol; + slake_vars = sol_slk; + int_sol = !is_int_sol; + epsilon = V.zero; + } + + let eval_eps (eps : V.t) (inf : R2.t) (sup : R2.t) = + let ({ R2.v = inf_r; offset = inf_d } : R2.t) = inf in + let ({ R2.v = sup_r; offset = sup_d } : R2.t) = sup in let c = V.compare inf_r sup_r in assert (c <= 0); if c = 0 || R.compare inf_d sup_d <= 0 then eps else V.min eps (V.div_by_coef (V.sub sup_r inf_r) (R.sub inf_d sup_d)) - let eps_for_mini - (i : Core.var_info) - (min : R2.t) - (eps : V.t) : V.t = + let eps_for_mini (i : Core.var_info) (min : R2.t) (eps : V.t) : V.t = assert (R2.is_pure_rational min || R.equal min.R2.offset R.one); let eps = eval_eps eps min i.value in eps - let eps_for_maxi - (i : Core.var_info) - (max : R2.t) - (eps : V.t) : V.t = + let eps_for_maxi (i : Core.var_info) (max : R2.t) (eps : V.t) : V.t = assert (R2.is_pure_rational max || R.equal max.R2.offset R.m_one); let eps = eval_eps eps i.value max in eps let get_rat_solution = let compute_epsilon mp eps = - MX.fold (fun _ (i, _) eps -> + MX.fold + (fun _ (i, _) eps -> let eps' = - match i.mini , i.maxi with + match (i.mini, i.maxi) with | None, None -> eps - - | Some min, None -> - eps_for_mini i min.bvalue eps - - | None, Some max -> - eps_for_maxi i max.bvalue eps - + | Some min, None -> eps_for_mini i min.bvalue eps + | None, Some max -> eps_for_maxi i max.bvalue eps | Some min, Some max -> - eps - |> eps_for_mini i min.bvalue - |> eps_for_maxi i max.bvalue + eps |> eps_for_mini i min.bvalue |> eps_for_maxi i max.bvalue in assert (V.compare eps' V.zero > 0); - eps' - ) mp eps + eps') + mp eps in let compute_solution slake mp eps acc = MX.fold (fun x (info, _) (m, s) -> - let {R2.v = q1; offset = q2} = info.value in - let q = V.add q1 (V.mult_by_coef eps q2) in - let q2 = R2.of_r q in - assert (not (violates_min_bound q2 info.mini)); - assert (not (violates_max_bound q2 info.maxi)); - if MX.mem x slake then m, (x, q) :: s else (x,q) :: m, s - )mp acc + let { R2.v = q1; offset = q2 } = info.value in + let q = V.add q1 (V.mult_by_coef eps q2) in + let q2 = R2.of_r q in + assert (not (violates_min_bound q2 info.mini)); + assert (not (violates_max_bound q2 info.maxi)); + if MX.mem x slake then (m, (x, q) :: s) else ((x, q) :: m, s)) + mp acc in fun env -> let eps = compute_epsilon env.basic V.one in let eps = compute_epsilon env.non_basic eps in let acc = compute_solution env.slake env.basic eps ([], []) in - let m,s = compute_solution env.slake env.non_basic eps acc in - { main_vars = m ; slake_vars = s; int_sol = false; epsilon = eps } - + let m, s = compute_solution env.slake env.non_basic eps acc in + { main_vars = m; slake_vars = s; int_sol = false; epsilon = eps } let get_solution env = if env.is_int then get_int_solution env else get_rat_solution env - let get_max_info {non_basic; _} p = + let get_max_info { non_basic; _ } p = let max_v : Core.bound = Core.P.fold (fun x c (max_v : Core.bound) -> - let xi, _ = try MX.find x non_basic with Not_found -> assert false in - let bnd = - if R.sign c > 0 - then xi.maxi - else xi.mini - in - let ex = match bnd with - | None -> Core.Ex.empty - | Some {explanation; _} -> explanation - in - let value = xi.value in - assert (equals_optimum value bnd); - { - bvalue = R2.add max_v.bvalue (R2.mult_by_const c value); - explanation = Ex.union max_v.explanation ex - } - ) + let xi, _ = + try MX.find x non_basic with Not_found -> assert false + in + let bnd = if R.sign c > 0 then xi.maxi else xi.mini in + let ex = + match bnd with + | None -> Core.Ex.empty + | Some { explanation; _ } -> explanation + in + let value = xi.value in + assert (equals_optimum value bnd); + { + bvalue = R2.add max_v.bvalue (R2.mult_by_const c value); + explanation = Ex.union max_v.explanation ex; + }) p - {bvalue = R2.zero; explanation = Ex.empty} + { bvalue = R2.zero; explanation = Ex.empty } in - {max_v; is_le = R.is_zero max_v.bvalue.R2.offset} - + { max_v; is_le = R.is_zero max_v.bvalue.R2.offset } let get opt env = match env.status with - | UNK -> Unknown + | UNK -> Unknown | UNSAT s -> Unsat (lazy (get_unsat_core s env)) - | SAT -> - match opt with - | None -> Sat (lazy (get_solution env)) - | Some(_, false) -> Unbounded (lazy (get_solution env)) - | Some(p, true) -> Max (lazy(get_max_info env p), lazy(get_solution env)) - + | SAT -> ( + match opt with + | None -> Sat (lazy (get_solution env)) + | Some (_, false) -> Unbounded (lazy (get_solution env)) + | Some (p, true) -> + Max (lazy (get_max_info env p), lazy (get_solution env))) end diff --git a/src/result.mli b/src/result.mli index 0e9122e..dedb1bd 100644 --- a/src/result.mli +++ b/src/result.mli @@ -6,6 +6,7 @@ module type S = sig module Core : CoreSig.S + val get : (Core.P.t * bool) option -> Core.t -> Core.result (** [get (objective, is_max_bounded) env] retrieves the result from a simplex [env]. @@ -17,8 +18,6 @@ module type S = sig @return solution that satisfies the constraints if any *) - - end -module Make(Core : CoreSig.S) : S with module Core = Core +module Make (Core : CoreSig.S) : S with module Core = Core diff --git a/src/solveBounds.ml b/src/solveBounds.ml index 9392eaf..49c3782 100644 --- a/src/solveBounds.ml +++ b/src/solveBounds.ml @@ -4,23 +4,22 @@ (* Copyright (C) --- OCamlPro --- See README.md for information and licensing *) (******************************************************************************) -let src = - Logs.Src.create "OcplibSimplex.SolveBounds" +let src = + Logs.Src.create "OcplibSimplex.SolveBounds" ~doc:"A module providing arithmetical utilities for solving the simplex" module type S = sig module Core : CoreSig.S + val solve : Core.t -> Core.t val maximize : Core.t -> Core.P.t -> Core.t * (Core.P.t * bool) option end -module Make(Core : CoreSig.S) : S with module Core = Core = struct - +module Make (Core : CoreSig.S) : S with module Core = Core = struct module Core = Core - module Result = Result.Make(Core) + module Result = Result.Make (Core) open Core - let gauss_pivot s p x c = let p, _ = P.replace s R.m_one (P.remove x p) in let c = R.div R.m_one c in @@ -32,30 +31,29 @@ module Make(Core : CoreSig.S) : S with module Core = Core = struct let look_for_next_pivot si pi non_basic = let status = si.vstatus in let is_lower = - match status with - | ValueOK -> assert false | LowerKO -> 1 | UpperKO -> -1 + match status with ValueOK -> assert false | LowerKO -> 1 | UpperKO -> -1 in try P.iter (fun x coef -> - let xi,use = try MX.find x non_basic with Not_found -> assert false in - let c = is_lower * R.sign coef in - assert (c <> 0); - if c > 0 && not (equals_optimum xi.value xi.maxi) then - raise (Out (x, coef, xi, use)); - if c < 0 && not (equals_optimum xi.value xi.mini) then - raise (Out (x, coef, xi, use)); - )pi; + let xi, use = + try MX.find x non_basic with Not_found -> assert false + in + let c = is_lower * R.sign coef in + assert (c <> 0); + if c > 0 && not (equals_optimum xi.value xi.maxi) then + raise (Out (x, coef, xi, use)); + if c < 0 && not (equals_optimum xi.value xi.mini) then + raise (Out (x, coef, xi, use))) + pi; None with Out (x, c, xi, use) -> Some (x, c, xi, use) - let adapt_valuation_of_newly_basic old_si new_si old_xi c_x = let diff = R2.div_by_const c_x (R2.sub new_si.value old_si.value) in { old_xi with value = R2.add diff old_xi.value } - -(* + (* let string_of_var_status stt = match stt with | P.Removed -> "Removed" @@ -68,59 +66,59 @@ module Make(Core : CoreSig.S) : S with module Core = Core = struct let rec solve_rec env round = Core.debug (Format.sprintf "[solve] round %d" round) env (Result.get None); Core.check_invariants env (Result.get None); - if SX.is_empty env.fixme then {env with status = SAT} + if SX.is_empty env.fixme then { env with status = SAT } else let s = SX.choose env.fixme in let fixme = SX.remove s env.fixme in let si, p = try MX.find s env.basic with Not_found -> assert false in match look_for_next_pivot si p env.non_basic with - | None -> {env with fixme = SX.empty; status = UNSAT s} - - | Some(x, c, xi, use_x) -> - Logs.debug ~src (fun p -> - p ~header:"[solve_rec]" - "pivot basic %a and non-basic %a@." - Var.print s Var.print x - ); - let basic = MX.remove s env.basic in - let non_basic = MX.remove x env.non_basic in - let q = gauss_pivot s p x c in - assert (SX.mem s use_x); - let use_x = SX.remove s use_x in - let old_si = si in - let si, changed = Core.ajust_value_of_non_basic si in - assert (changed); - - let old_xi = xi in - let xi = adapt_valuation_of_newly_basic old_si si xi c in - let xi = ajust_status_of_basic xi in - let diff_xi_val = R2.sub xi.value old_xi.value in - let fixme = (* do this earlier to detect bad pivots *) - if xi.vstatus == ValueOK then fixme - else SX.add x fixme - in - let non_basic = - P.fold - (fun y _ non_basic -> - let yi, use_y = - try MX.find y non_basic with Not_found -> assert false in - MX.add y (yi, SX.add x (SX.remove s use_y)) non_basic - )(P.remove s q) non_basic - in - - let non_basic = MX.add s (si, SX.add x use_x) non_basic in - - let basic, non_basic, fixme = - SX.fold - (fun t (basic, non_basic, fixme) -> - let ti0, r = try MX.find t basic with Not_found -> assert false in - let cx = try P.find x r with Not_found -> assert false in - (*should update_ti*) - let diff_cx = R2.mult_by_const cx diff_xi_val in - let ti = {ti0 with value = R2.add ti0.value diff_cx} in - let ti = ajust_status_of_basic ti in - let r', changed = P.subst x q r in - (* + | None -> { env with fixme = SX.empty; status = UNSAT s } + | Some (x, c, xi, use_x) -> + Logs.debug ~src (fun p -> + p ~header:"[solve_rec]" "pivot basic %a and non-basic %a@." + Var.print s Var.print x); + let basic = MX.remove s env.basic in + let non_basic = MX.remove x env.non_basic in + let q = gauss_pivot s p x c in + assert (SX.mem s use_x); + let use_x = SX.remove s use_x in + let old_si = si in + let si, changed = Core.ajust_value_of_non_basic si in + assert changed; + + let old_xi = xi in + let xi = adapt_valuation_of_newly_basic old_si si xi c in + let xi = ajust_status_of_basic xi in + let diff_xi_val = R2.sub xi.value old_xi.value in + let fixme = + (* do this earlier to detect bad pivots *) + if xi.vstatus == ValueOK then fixme else SX.add x fixme + in + let non_basic = + P.fold + (fun y _ non_basic -> + let yi, use_y = + try MX.find y non_basic with Not_found -> assert false + in + MX.add y (yi, SX.add x (SX.remove s use_y)) non_basic) + (P.remove s q) non_basic + in + + let non_basic = MX.add s (si, SX.add x use_x) non_basic in + + let basic, non_basic, fixme = + SX.fold + (fun t (basic, non_basic, fixme) -> + let ti0, r = + try MX.find t basic with Not_found -> assert false + in + let cx = try P.find x r with Not_found -> assert false in + (*should update_ti*) + let diff_cx = R2.mult_by_const cx diff_xi_val in + let ti = { ti0 with value = R2.add ti0.value diff_cx } in + let ti = ajust_status_of_basic ti in + let r', changed = P.subst x q r in + (* Format.eprintf "update poly of basic %a@." Var.print t; List.iter (fun (v, vstt) -> @@ -128,49 +126,46 @@ module Make(Core : CoreSig.S) : S with module Core = Core = struct Var.print v (string_of_var_status vstt); )changed; *) - let non_basic = - List.fold_left - (fun non_basic (z, vstt) -> + let non_basic = + List.fold_left + (fun non_basic (z, vstt) -> match vstt with | P.Exists -> non_basic | P.New -> - let zi, use_z = - try MX.find z non_basic with Not_found -> assert false - in - MX.add z (zi, SX.add t use_z) non_basic - - | P.Removed -> - if Var.compare z x = 0 then non_basic - else let zi, use_z = - try MX.find z non_basic with Not_found -> assert false + try MX.find z non_basic + with Not_found -> assert false in - MX.add z (zi, SX.remove t use_z) non_basic - )non_basic changed - in - (*val subst : Var.t -> t -> t -> t * (Var.t * var_status) list*) - - let basic = MX.add t (ti, r') basic in - let fixme = - if ti.vstatus == ValueOK then - if ti0.vstatus == ValueOK then fixme - else SX.remove t fixme - else SX.add t fixme - in - basic, non_basic, fixme - )use_x (basic, non_basic, fixme) - in - - (* ... *) - - let basic = MX.add x (xi, q) basic in - - (* ... *) - - let env = {env with fixme; basic; non_basic} in - env.nb_pivots := !(env.nb_pivots) + 1; - solve_rec env (round + 1) - + MX.add z (zi, SX.add t use_z) non_basic + | P.Removed -> + if Var.compare z x = 0 then non_basic + else + let zi, use_z = + try MX.find z non_basic + with Not_found -> assert false + in + MX.add z (zi, SX.remove t use_z) non_basic) + non_basic changed + in + + (*val subst : Var.t -> t -> t -> t * (Var.t * var_status) list*) + let basic = MX.add t (ti, r') basic in + let fixme = + if ti.vstatus == ValueOK then + if ti0.vstatus == ValueOK then fixme else SX.remove t fixme + else SX.add t fixme + in + (basic, non_basic, fixme)) + use_x (basic, non_basic, fixme) + in + + (* ... *) + let basic = MX.add x (xi, q) basic in + + (* ... *) + let env = { env with fixme; basic; non_basic } in + env.nb_pivots := !(env.nb_pivots) + 1; + solve_rec env (round + 1) let solve env = Core.debug "[entry of solve]" env (Result.get None); @@ -184,311 +179,295 @@ module Make(Core : CoreSig.S) : S with module Core = Core = struct Core.check_invariants env (Result.get None); env - - - - - - let non_basic_to_maximize {non_basic=n_b; _} opt = + let non_basic_to_maximize { non_basic = n_b; _ } opt = let acc = ref None in try P.iter (fun x c -> - let xi, use = try MX.find x n_b with Not_found -> assert false in - let sg = R.sign c in - if sg > 0 && not (equals_optimum xi.value xi.maxi) || - sg < 0 && not (equals_optimum xi.value xi.mini) then begin - acc := Some (x, c, xi, use, sg > 0); - raise Exit - end - )opt; + let xi, use = try MX.find x n_b with Not_found -> assert false in + let sg = R.sign c in + if + (sg > 0 && not (equals_optimum xi.value xi.maxi)) + || (sg < 0 && not (equals_optimum xi.value xi.mini)) + then ( + acc := Some (x, c, xi, use, sg > 0); + raise Exit)) + opt; !acc with Exit -> !acc - - type 'a maximiza_basic = - | Free - | Stuck - | Progress of 'a - + type 'a maximiza_basic = Free | Stuck | Progress of 'a let basic_var_to_pivot_for_maximization = let choose_best_pivot acc s si p c_px bnd_opt is_min = match bnd_opt with - | None -> - if !acc = Stuck then acc := Free (* !!! to check *) - - | Some {bvalue = bnd; _} -> - let tmp = if is_min then R2.sub si.value bnd else R2.sub bnd si.value in - let ratio = R2.div_by_const (R.abs c_px) tmp in - begin - match !acc with - | Free | Stuck -> - acc := Progress (ratio, s, si, p, c_px, bnd, is_min) - - | Progress (r_old,_,_,_,_,_,_) -> - if R2.compare r_old ratio > 0 then - acc := Progress (ratio, s, si, p, c_px, bnd, is_min) - end; - if R2.is_zero ratio then raise Exit (* in the case, the pivot is found*) + | None -> if !acc = Stuck then acc := Free (* !!! to check *) + | Some { bvalue = bnd; _ } -> + let tmp = + if is_min then R2.sub si.value bnd else R2.sub bnd si.value + in + let ratio = R2.div_by_const (R.abs c_px) tmp in + (match !acc with + | Free | Stuck -> acc := Progress (ratio, s, si, p, c_px, bnd, is_min) + | Progress (r_old, _, _, _, _, _, _) -> + if R2.compare r_old ratio > 0 then + acc := Progress (ratio, s, si, p, c_px, bnd, is_min)); + if R2.is_zero ratio then + raise Exit (* in the case, the pivot is found*) in - fun {basic; _} x use_x should_incr_x -> + fun { basic; _ } x use_x should_incr_x -> (* Initially, we assume that we are stuck, unless, use_x is empty *) let acc = ref (if SX.is_empty use_x then Free else Stuck) in try SX.iter (fun s -> - let si, p = try MX.find s basic with Not_found -> assert false in - let c_px = try P.find x p with Not_found -> assert false in - let sg = R.sign c_px in - assert (sg <> 0); - match should_incr_x, sg > 0, si.mini, si.maxi with - | true , true , _, mx_opt -> - (* by increasing x, s will increase and max(s) <> +infty *) - choose_best_pivot acc s si p c_px mx_opt false - - | true , false, mn_opt, _ -> - (* by increasing x, s will decrease and min(s) <> -infty *) - choose_best_pivot acc s si p c_px mn_opt true - - | false, true , mn_opt, _ -> - (* by decreasing x, s will decreease and min(s) <> -infty *) - choose_best_pivot acc s si p c_px mn_opt true - - | false, false, _, mx_opt -> - (* by decreasning x, s will increase and max(s) <> +infty *) - choose_best_pivot acc s si p c_px mx_opt false - - (*| true, true, _, None - | true, false, None, _ - | false, true, None, _ - | false, false, _, None -> - (* for the cases where max or max = infty, we keep acc unchanged. - if acc = None at the end, the problem is unbounded *) - () - *) - )use_x; + let si, p = try MX.find s basic with Not_found -> assert false in + let c_px = try P.find x p with Not_found -> assert false in + let sg = R.sign c_px in + assert (sg <> 0); + match (should_incr_x, sg > 0, si.mini, si.maxi) with + | true, true, _, mx_opt -> + (* by increasing x, s will increase and max(s) <> +infty *) + choose_best_pivot acc s si p c_px mx_opt false + | true, false, mn_opt, _ -> + (* by increasing x, s will decrease and min(s) <> -infty *) + choose_best_pivot acc s si p c_px mn_opt true + | false, true, mn_opt, _ -> + (* by decreasing x, s will decreease and min(s) <> -infty *) + choose_best_pivot acc s si p c_px mn_opt true + | false, false, _, mx_opt -> + (* by decreasning x, s will increase and max(s) <> +infty *) + choose_best_pivot acc s si p c_px mx_opt false + (*| true, true, _, None + | true, false, None, _ + | false, true, None, _ + | false, false, _, None -> + (* for the cases where max or max = infty, we keep acc unchanged. + if acc = None at the end, the problem is unbounded *) + () + *)) + use_x; !acc with Exit -> !acc - let can_fix_valuation_without_pivot should_incr xi ratio_opt = if should_incr then - match xi.maxi, ratio_opt with + match (xi.maxi, ratio_opt) with | None, _ -> None - | Some {bvalue = bnd; _}, Some ratio -> - let diff = R2.sub bnd xi.value in - if R2.compare diff ratio < 0 then Some ({xi with value = bnd}, diff) - else None - - | Some {bvalue = bnd; _}, None -> - let diff = R2.sub bnd xi.value in - Some ({xi with value = bnd}, diff) - + | Some { bvalue = bnd; _ }, Some ratio -> + let diff = R2.sub bnd xi.value in + if R2.compare diff ratio < 0 then Some ({ xi with value = bnd }, diff) + else None + | Some { bvalue = bnd; _ }, None -> + let diff = R2.sub bnd xi.value in + Some ({ xi with value = bnd }, diff) else - match xi.mini, ratio_opt with + match (xi.mini, ratio_opt) with | None, _ -> None - | Some {bvalue = bnd; _}, Some ratio -> - let diff = R2.sub xi.value bnd in - if R2.compare diff ratio < 0 then Some ({xi with value = bnd}, diff) - else None - - | Some {bvalue = bnd; _}, None -> - let diff = R2.sub xi.value bnd in - Some ({xi with value = bnd}, diff) - - - let update_valuation_without_pivot - ({basic; non_basic; _ } as env) x use_x new_xi diff _should_incr = + | Some { bvalue = bnd; _ }, Some ratio -> + let diff = R2.sub xi.value bnd in + if R2.compare diff ratio < 0 then Some ({ xi with value = bnd }, diff) + else None + | Some { bvalue = bnd; _ }, None -> + let diff = R2.sub xi.value bnd in + Some ({ xi with value = bnd }, diff) + + let update_valuation_without_pivot ({ basic; non_basic; _ } as env) x use_x + new_xi diff _should_incr = let non_basic = MX.add x (new_xi, use_x) non_basic in let diff = if _should_incr then diff else R2.minus diff in let basic = SX.fold (fun s basic -> - let si, p = try MX.find s basic with Not_found -> assert false in - let cx = try P.find x p with Not_found -> assert false in - assert (not (R.is_zero cx)); - let delta = R2.mult_by_const cx diff in - let si = {si with value = R2.add si.value delta} in - MX.add s (si, p) basic - )use_x basic + let si, p = try MX.find s basic with Not_found -> assert false in + let cx = try P.find x p with Not_found -> assert false in + assert (not (R.is_zero cx)); + let delta = R2.mult_by_const cx diff in + let si = { si with value = R2.add si.value delta } in + MX.add s (si, p) basic) + use_x basic in - {env with basic; non_basic} + { env with basic; non_basic } let rec maximize_rec env opt rnd = - Logs.debug ~src (fun p -> p "[maximize_rec] round %d // OPT = %a@." rnd P.print opt); + Logs.debug ~src (fun p -> + p "[maximize_rec] round %d // OPT = %a@." rnd P.print opt); Core.debug - (Format.sprintf "[maximize_rec] round %d" rnd) env (Result.get None); + (Format.sprintf "[maximize_rec] round %d" rnd) + env (Result.get None); Core.check_invariants env (Result.get None); match non_basic_to_maximize env opt with - | None -> Logs.debug (fun p -> p "max reached@."); - rnd, env, Some (opt, true) (* max reached *) - - | Some (_x, _c, _xi, _use_x, _should_incr) -> - Logs.debug (fun p -> p "pivot non basic var %a ?@." Var.print _x); - match basic_var_to_pivot_for_maximization env _x _use_x _should_incr with - | Free -> - Logs.debug (fun p -> p "non basic %a not constrained by basic vars: Set it to max@." - Var.print _x); - begin - match can_fix_valuation_without_pivot _should_incr _xi None with - | Some (new_xi, diff) -> - Logs.debug (fun p -> p - "No --> I can set value of %a to min/max WO pivot@." - Var.print _x); + | None -> + Logs.debug (fun p -> p "max reached@."); + (rnd, env, Some (opt, true)) (* max reached *) + | Some (_x, _c, _xi, _use_x, _should_incr) -> ( + Logs.debug (fun p -> p "pivot non basic var %a ?@." Var.print _x); + match + basic_var_to_pivot_for_maximization env _x _use_x _should_incr + with + | Free -> ( + Logs.debug (fun p -> + p "non basic %a not constrained by basic vars: Set it to max@." + Var.print _x); + match can_fix_valuation_without_pivot _should_incr _xi None with + | Some (new_xi, diff) -> + Logs.debug (fun p -> + p "No --> I can set value of %a to min/max WO pivot@." + Var.print _x); + let env, opt = + ( update_valuation_without_pivot env _x _use_x new_xi diff + _should_incr, + opt ) + in + (* no pivot *) + maximize_rec env opt (rnd + 1) + | None -> + Logs.debug (fun p -> + p "no pivot finally(no upper bnd), pb unbounded@."); + (rnd, env, Some (opt, false)) + (* unbounded *)) + | Stuck -> + Logs.debug (fun p -> p "no pivot finally, pb unbounded@."); + (rnd, env, Some (opt, false)) (* unbounded *) + | Progress (ratio, s, si, p, c_px, bnd, _is_min) -> + Logs.debug (fun p -> p "pivot with basic var %a ?@." Var.print s); let env, opt = - update_valuation_without_pivot - env _x _use_x new_xi diff _should_incr, opt - in - (* no pivot *) - maximize_rec env opt (rnd + 1) - | None -> - Logs.debug (fun p -> p - "no pivot finally(no upper bnd), pb unbounded@."); - rnd, env, Some (opt, false) (* unbounded *) - end - | Stuck -> - Logs.debug (fun p -> p "no pivot finally, pb unbounded@."); - rnd, env, Some (opt, false) (* unbounded *) - - | Progress (ratio, s, si, p, c_px, bnd, _is_min) -> - Logs.debug (fun p -> p "pivot with basic var %a ?@." Var.print s); - let env, opt = - match - can_fix_valuation_without_pivot _should_incr _xi (Some ratio) with - | Some (new_xi, diff) -> - Logs.debug (fun p -> p - "No --> I can set value of %a to min/max WO pivot@." - Var.print _x); - update_valuation_without_pivot - env _x _use_x new_xi diff _should_incr, opt - - | None -> - let x = _x in - let c = c_px in - let use_x = _use_x in - let xi = _xi in - Logs.debug (fun p -> p - "[maximize_rec] pivot basic %a and non-basic %a@." - Var.print s Var.print x); - let basic = MX.remove s env.basic in - let non_basic = MX.remove x env.non_basic in - let q = gauss_pivot s p x c in - assert (SX.mem s use_x); - let use_x = SX.remove s use_x in - let old_si = si in - - let si = {si with value = bnd} in (* difference wrt solve *) - (* + match + can_fix_valuation_without_pivot _should_incr _xi (Some ratio) + with + | Some (new_xi, diff) -> + Logs.debug (fun p -> + p "No --> I can set value of %a to min/max WO pivot@." + Var.print _x); + ( update_valuation_without_pivot env _x _use_x new_xi diff + _should_incr, + opt ) + | None -> + let x = _x in + let c = c_px in + let use_x = _use_x in + let xi = _xi in + Logs.debug (fun p -> + p "[maximize_rec] pivot basic %a and non-basic %a@." + Var.print s Var.print x); + let basic = MX.remove s env.basic in + let non_basic = MX.remove x env.non_basic in + let q = gauss_pivot s p x c in + assert (SX.mem s use_x); + let use_x = SX.remove s use_x in + let old_si = si in + + let si = { si with value = bnd } in + + (* difference wrt solve *) + (* because the code of solve below, assumes that value in si violotas a bound let si, changed = Core.ajust_value_of_non_basic si in assert (changed); *) + let old_xi = xi in + let xi = adapt_valuation_of_newly_basic old_si si xi c in + let xi = ajust_status_of_basic xi in + let diff_xi_val = R2.sub xi.value old_xi.value in + assert (xi.vstatus == ValueOK); + let non_basic = + P.fold + (fun y _ non_basic -> + let yi, use_y = + try MX.find y non_basic + with Not_found -> assert false + in + MX.add y (yi, SX.add x (SX.remove s use_y)) non_basic) + (P.remove s q) non_basic + in - let old_xi = xi in - let xi = adapt_valuation_of_newly_basic old_si si xi c in - let xi = ajust_status_of_basic xi in - let diff_xi_val = R2.sub xi.value old_xi.value in - assert(xi.vstatus == ValueOK); - let non_basic = - P.fold - (fun y _ non_basic -> - let yi, use_y = - try MX.find y non_basic with Not_found -> assert false in - MX.add y (yi, SX.add x (SX.remove s use_y)) non_basic - )(P.remove s q) non_basic - in - - let non_basic = MX.add s (si, SX.add x use_x) non_basic in - - let basic, non_basic = - SX.fold - (fun t (basic, non_basic) -> - let ti0, r = - try MX.find t basic with Not_found -> assert false in - let cx = try P.find x r with Not_found -> assert false in - (*should update_ti*) - let diff_cx = R2.mult_by_const cx diff_xi_val in - let ti = {ti0 with value = R2.add ti0.value diff_cx} in - let ti = ajust_status_of_basic ti in - let r', changed = P.subst x q r in - - let non_basic = - List.fold_left - (fun non_basic (z, vstt) -> - match vstt with - | P.Exists -> non_basic - | P.New -> - let zi, use_z = - try MX.find z non_basic - with Not_found -> assert false - in - MX.add z (zi, SX.add t use_z) non_basic - - | P.Removed -> - if Var.compare z x = 0 then non_basic - else - let zi, use_z = - try MX.find z non_basic - with Not_found -> assert false - in - MX.add z (zi, SX.remove t use_z) non_basic - )non_basic changed - in - - let basic = MX.add t (ti, r') basic in - assert(ti.vstatus == ValueOK); - basic, non_basic - )use_x (basic, non_basic) - in - - (* ... *) - let basic = MX.add x (xi, q) basic in + let non_basic = MX.add s (si, SX.add x use_x) non_basic in - (* ... *) - {env with basic; non_basic}, (fst (P.subst x q opt)) + let basic, non_basic = + SX.fold + (fun t (basic, non_basic) -> + let ti0, r = + try MX.find t basic with Not_found -> assert false + in + let cx = + try P.find x r with Not_found -> assert false + in + (*should update_ti*) + let diff_cx = R2.mult_by_const cx diff_xi_val in + let ti = + { ti0 with value = R2.add ti0.value diff_cx } + in + let ti = ajust_status_of_basic ti in + let r', changed = P.subst x q r in + + let non_basic = + List.fold_left + (fun non_basic (z, vstt) -> + match vstt with + | P.Exists -> non_basic + | P.New -> + let zi, use_z = + try MX.find z non_basic + with Not_found -> assert false + in + MX.add z (zi, SX.add t use_z) non_basic + | P.Removed -> + if Var.compare z x = 0 then non_basic + else + let zi, use_z = + try MX.find z non_basic + with Not_found -> assert false + in + MX.add z (zi, SX.remove t use_z) non_basic) + non_basic changed + in + let basic = MX.add t (ti, r') basic in + assert (ti.vstatus == ValueOK); + (basic, non_basic)) + use_x (basic, non_basic) + in + (* ... *) + let basic = MX.add x (xi, q) basic in - in - env.nb_pivots := !(env.nb_pivots) + 1; - maximize_rec env opt (rnd + 1) + (* ... *) + ({ env with basic; non_basic }, fst (P.subst x q opt)) + in + env.nb_pivots := !(env.nb_pivots) + 1; + maximize_rec env opt (rnd + 1)) let maximize env opt0 = let env = solve env in match env.status with | UNK -> assert false - | UNSAT _ -> env, None + | UNSAT _ -> (env, None) | SAT -> - Logs.debug (fun p -> p "[maximize] pb SAT! try to maximize %a@." P.print opt0); - let {basic; non_basic; _} = env in - let unbnd = ref false in - let opt = - P.fold - (fun x c acc -> - if MX.mem x non_basic then fst (P.accumulate x c acc) - else - try fst (P.append acc c (snd (MX.find x basic))) - with Not_found -> - unbnd := true; - fst (P.accumulate x c acc) - )opt0 P.empty - in - if !unbnd then env, Some (opt, false) (* unbounded *) - else - begin + Logs.debug (fun p -> + p "[maximize] pb SAT! try to maximize %a@." P.print opt0); + let { basic; non_basic; _ } = env in + let unbnd = ref false in + let opt = + P.fold + (fun x c acc -> + if MX.mem x non_basic then fst (P.accumulate x c acc) + else + try fst (P.append acc c (snd (MX.find x basic))) + with Not_found -> + unbnd := true; + fst (P.accumulate x c acc)) + opt0 P.empty + in + if !unbnd then (env, Some (opt, false)) (* unbounded *) + else ( Logs.debug (fun p -> p "start maximization@."); let rnd, env, is_max = maximize_rec env opt 1 in Core.check_invariants env (Result.get is_max); - Logs.debug (fun p -> p "[maximize] pb SAT! Max found ? %b for %a == %a@." - (is_max != None) P.print opt0 P.print opt); + Logs.debug (fun p -> + p "[maximize] pb SAT! Max found ? %b for %a == %a@." + (is_max != None) P.print opt0 P.print opt); Logs.debug (fun p -> p "maximization done after %d steps@." rnd); - env, is_max - end - - - - + (env, is_max)) end diff --git a/src/solveBounds.mli b/src/solveBounds.mli index f7842df..ca0cf96 100644 --- a/src/solveBounds.mli +++ b/src/solveBounds.mli @@ -24,4 +24,4 @@ module type S = sig *) end -module Make(Core : CoreSig.S) : S with module Core = Core +module Make (Core : CoreSig.S) : S with module Core = Core diff --git a/src/version.ml b/src/version.ml index 00c77dc..b6ddd08 100644 --- a/src/version.ml +++ b/src/version.ml @@ -4,4 +4,4 @@ (* Copyright (C) --- OCamlPro --- See README.md for information and licensing *) (******************************************************************************) -let version="0.5" +let version = "0.5" diff --git a/tests/dune b/tests/dune index 4fcdd22..42e1428 100644 --- a/tests/dune +++ b/tests/dune @@ -1,51 +1,28 @@ (tests - (names - solveEmpty - standalone_minimal_maximization - standalone_minimal - standalone_test_strict_bounds) - (modules - simplex - solveEmpty - standalone_minimal_maximization - standalone_minimal + (names solveEmpty standalone_minimal_maximization standalone_minimal standalone_test_strict_bounds) + (modules simplex solveEmpty standalone_minimal_maximization + standalone_minimal standalone_test_strict_bounds) (libraries ocplib-simplex zarith)) (rule (alias runtest) (action - (diff solveEmpty.expected solveEmpty.output) - ) -) + (diff solveEmpty.expected solveEmpty.output))) (rule (alias runtest) (action - (diff - standalone_minimal_maximization.expected - standalone_minimal_maximization.output - ) - ) -) + (diff standalone_minimal_maximization.expected + standalone_minimal_maximization.output))) (rule (alias runtest) (action - (diff - standalone_minimal.expected - standalone_minimal.output - ) - ) -) + (diff standalone_minimal.expected standalone_minimal.output))) (rule (alias runtest) (action - (diff - standalone_test_strict_bounds.expected - standalone_test_strict_bounds.output - ) - ) -) - + (diff standalone_test_strict_bounds.expected + standalone_test_strict_bounds.output))) diff --git a/tests/simplex.ml b/tests/simplex.ml index 65f295d..36286e6 100644 --- a/tests/simplex.ml +++ b/tests/simplex.ml @@ -8,14 +8,13 @@ module Var = struct type t = string let print fmt s = Format.fprintf fmt "%s" s - let compare = String.compare - let is_int _ = true end module Rat = struct type t = Q.t + let add = Q.add let minus = Q.neg let mult = Q.mul @@ -27,7 +26,6 @@ module Rat = struct let m_one = Q.minus_one let is_zero n = Q.equal n Q.zero let to_string = Q.to_string - let print = Q.pp_print let is_int v = Z.equal (Q.den v) Z.one let div = Q.div @@ -36,31 +34,28 @@ module Rat = struct let is_m_one v = Q.equal v Q.minus_one let sign = Q.sign let min = Q.min - let floor v = Z.fdiv (Q.num v) (Q.den v) |> Q.of_bigint let ceiling v = Z.cdiv (Q.num v) (Q.den v) |> Q.of_bigint - end module Ex = struct - include Set.Make(String) + include Set.Make (String) - let print fmt s = match elements s with + let print fmt s = + match elements s with | [] -> Format.fprintf fmt "()" - | e::l -> - Format.fprintf fmt "%s" e; - List.iter (Format.fprintf fmt ", %s") l + | e :: l -> + Format.fprintf fmt "%s" e; + List.iter (Format.fprintf fmt ", %s") l end -module Ty = OcplibSimplex.Core.Make(Var)(Rat)(Ex) -module AB = OcplibSimplex.AssertBounds.Make(Ty) - -module Sim = OcplibSimplex.Basic.Make(Var)(Rat)(Ex) +module Ty = OcplibSimplex.Core.Make (Var) (Rat) (Ex) +module AB = OcplibSimplex.AssertBounds.Make (Ty) +module Sim = OcplibSimplex.Basic.Make (Var) (Rat) (Ex) let aux header (sim, opt) = Format.printf "%a" header (); - Format.printf "%a" - (Sim.Core.print (Sim.Result.get opt sim)) sim + Format.printf "%a" (Sim.Core.print (Sim.Result.get opt sim)) sim let () = Logs.Src.set_level OcplibSimplex.Core.src (Some Debug); diff --git a/tests/solveEmpty.ml b/tests/solveEmpty.ml index 9666092..848be05 100644 --- a/tests/solveEmpty.ml +++ b/tests/solveEmpty.ml @@ -1,10 +1,6 @@ - open Simplex let () = let sim = Sim.Core.empty ~is_int:false ~check_invs:true in let sim = Sim.Solve.solve sim in - aux ( - fun fmt () -> - Format.fprintf fmt "\n### Test Solve Empty@." - ) (sim, None) + aux (fun fmt () -> Format.fprintf fmt "\n### Test Solve Empty@.") (sim, None) diff --git a/tests/standalone_minimal.ml b/tests/standalone_minimal.ml index 7a32efd..0350244 100644 --- a/tests/standalone_minimal.ml +++ b/tests/standalone_minimal.ml @@ -1,32 +1,28 @@ - open Simplex let () = let sim = Sim.Core.empty ~is_int:true ~check_invs:true in let zero = Sim.Core.R2.zero in - let m_one = (Sim.Core.R2.of_r Rat.m_one) in + let m_one = Sim.Core.R2.of_r Rat.m_one in (* x >= 0 *) let sim, _ = Sim.Assert.var sim "x" - ~min:{Sim.Core.bvalue = zero; explanation = Ex.singleton "x>=0"} + ~min:{ Sim.Core.bvalue = zero; explanation = Ex.singleton "x>=0" } in (* y >= 0 *) let sim, _ = Sim.Assert.var sim "y" - ~min:{Sim.Core.bvalue = zero; explanation = Ex.singleton "y>=0"} + ~min:{ Sim.Core.bvalue = zero; explanation = Ex.singleton "y>=0" } in - let x_y = Sim.Core.P.from_list ["x", Rat.one; "y", Rat.one] in + let x_y = Sim.Core.P.from_list [ ("x", Rat.one); ("y", Rat.one) ] in (* z == x + y <= -1 *) let sim, _ = Sim.Assert.poly sim x_y "z" - ~max:{Sim.Core.bvalue = m_one; explanation = Ex.singleton "x+y<=-1"} + ~max:{ Sim.Core.bvalue = m_one; explanation = Ex.singleton "x+y<=-1" } in let sim = Sim.Solve.solve sim in - aux ( - fun fmt () -> - Format.fprintf fmt "\n### Test Solve Unsat@." - ) (sim, None) + aux (fun fmt () -> Format.fprintf fmt "\n### Test Solve Unsat@.") (sim, None) diff --git a/tests/standalone_minimal_maximization.ml b/tests/standalone_minimal_maximization.ml index b19327b..3a59893 100644 --- a/tests/standalone_minimal_maximization.ml +++ b/tests/standalone_minimal_maximization.ml @@ -1,33 +1,30 @@ - open Simplex let large i = Sim.Core.R2.of_r (Q.of_int i) let upper i = Sim.Core.R2.upper (Q.of_int i) let lower i = Sim.Core.R2.lower (Q.of_int i) - -let bnd r e = {Sim.Core.bvalue = r; explanation = e} +let bnd r e = { Sim.Core.bvalue = r; explanation = e } let () = let sim = Sim.Core.empty ~is_int:true ~check_invs:true in - let x_y = Sim.Core.P.from_list ["x", Rat.one; "y", Rat.one] in - let y1 = Sim.Core.P.from_list ["y", Rat.one] in - let ym1 = Sim.Core.P.from_list ["y", Rat.m_one] in + let x_y = Sim.Core.P.from_list [ ("x", Rat.one); ("y", Rat.one) ] in + let y1 = Sim.Core.P.from_list [ ("y", Rat.one) ] in + let ym1 = Sim.Core.P.from_list [ ("y", Rat.m_one) ] in (* s == x + y >= 10 - let sim = Sim.Assert.poly sim x_y "s" (large 10) Ex.empty None Ex.empty in + let sim = Sim.Assert.poly sim x_y "s" (large 10) Ex.empty None Ex.empty in *) (* x <= 5 *) let sim, _ = - Sim.Assert.var sim "x" - ~min:(bnd (large 3) (Ex.singleton "x>=3")) + Sim.Assert.var sim "x" ~min:(bnd (large 3) (Ex.singleton "x>=3")) in (* s == x + y <= 10 *) let sim, _ = - Sim.Assert.poly sim x_y "s" - ~max:(bnd (large 10) (Ex.singleton "x+y<=10")) in + Sim.Assert.poly sim x_y "s" ~max:(bnd (large 10) (Ex.singleton "x+y<=10")) + in let max_hdr pb fmt () = Format.fprintf fmt "### Problem 'max %a'@." Sim.Core.P.print pb diff --git a/tests/standalone_test_strict_bounds.ml b/tests/standalone_test_strict_bounds.ml index 6a339fd..b9bf6e5 100644 --- a/tests/standalone_test_strict_bounds.ml +++ b/tests/standalone_test_strict_bounds.ml @@ -5,8 +5,7 @@ let sep () = let pp_epsilon fmt (max_v, eps) = let pp = - if Sim.Core.R2.is_pure_rational max_v.Sim.Core.bvalue - then Format.ifprintf + if Sim.Core.R2.is_pure_rational max_v.Sim.Core.bvalue then Format.ifprintf else Format.fprintf in pp fmt "(epsilon: %a)" Rat.print eps @@ -15,43 +14,33 @@ let aux sim opt_p = let sim, opt = Sim.Solve.maximize sim opt_p in sep (); Format.printf "The problem 'max %a' ...@." Sim.Core.P.print opt_p; - begin - match Sim.Result.get opt sim with - | Sim.Core.Unknown -> assert false - | Sim.Core.Sat _ -> assert false - - | Sim.Core.Unsat ex -> - Format.printf " is unsat (reason = %a)@." Ex.print (Lazy.force ex); - - | Sim.Core.Unbounded _ -> Format.printf " is unbounded@." - - | Sim.Core.Max (mx,sol) -> - let {Sim.Core.max_v; is_le} = Lazy.force mx in + (match Sim.Result.get opt sim with + | Sim.Core.Unknown -> assert false + | Sim.Core.Sat _ -> assert false + | Sim.Core.Unsat ex -> + Format.printf " is unsat (reason = %a)@." Ex.print (Lazy.force ex) + | Sim.Core.Unbounded _ -> Format.printf " is unbounded@." + | Sim.Core.Max (mx, sol) -> + let { Sim.Core.max_v; is_le } = Lazy.force mx in let sol = Lazy.force sol in - Format.printf - " has an upper bound: %a (is_le = %b)(reason: %a)%a@." - Sim.Core.R2.print max_v.Sim.Core.bvalue - is_le - Ex.print max_v.Sim.Core.explanation - pp_epsilon (max_v, sol.Sim.Core.epsilon) - ; - end; + Format.printf " has an upper bound: %a (is_le = %b)(reason: %a)%a@." + Sim.Core.R2.print max_v.Sim.Core.bvalue is_le Ex.print + max_v.Sim.Core.explanation pp_epsilon + (max_v, sol.Sim.Core.epsilon)); sep (); Format.printf "@." let large i = Sim.Core.R2.of_r (Q.of_int i) let upper i = Sim.Core.R2.upper (Q.of_int i) let lower i = Sim.Core.R2.lower (Q.of_int i) - -let bnd r e = {Sim.Core.bvalue = r; explanation = e} - +let bnd r e = { Sim.Core.bvalue = r; explanation = e } let r_two = Rat.add Rat.one Rat.one let () = let sim = Sim.Core.empty ~is_int:true ~check_invs:false in - let x_m_y = Sim.Core.P.from_list ["x", Rat.one; "y", Rat.m_one] in - let tx_ty = Sim.Core.P.from_list ["x", r_two; "y", r_two] in + let x_m_y = Sim.Core.P.from_list [ ("x", Rat.one); ("y", Rat.m_one) ] in + let tx_ty = Sim.Core.P.from_list [ ("x", r_two); ("y", r_two) ] in (* 3 < y < 5*) let sim, _ = @@ -69,13 +58,13 @@ let () = (* 0 <= x - y *) let sim, _ = - Sim.Assert.poly sim x_m_y "s'" - ~min:(bnd (large 1) (Ex.singleton "x-y>=1")) + Sim.Assert.poly sim x_m_y "s'" ~min:(bnd (large 1) (Ex.singleton "x-y>=1")) in (* s == 2x + 2y <= 20 *) let sim, _ = Sim.Assert.poly sim tx_ty "s" - ~max:(bnd (large 20) (Ex.singleton "2x+2y<=20")) in + ~max:(bnd (large 20) (Ex.singleton "2x+2y<=20")) + in aux sim tx_ty