diff --git a/src/prolog/lib/clpb.pl b/src/prolog/lib/clpb.pl new file mode 100644 index 00000000..0ddfc986 --- /dev/null +++ b/src/prolog/lib/clpb.pl @@ -0,0 +1,1775 @@ +/* CLP(B): Constraint Logic Programming over Boolean Variables + + Copyright (C): 2019 Markus Triska + All rights reserved. + + E-mail: triska@metalevel.at + WWW: http://www.metalevel.at + +*/ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + CLP(B): Constraint Logic Programming over Boolean Variables. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Public operators. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +:- op(300, fy, ~). +:- op(500, yfx, #). + +:- module(clpb, [ + sat/1, + taut/2, + labeling/1, + sat_count/2, + weighted_maximum/3, + random_labeling/2 + ]). + +:- use_module(library(assoc)). +:- use_module(library(between)). +:- use_module(library(atts)). +:- use_module(library(lists)). +:- use_module(library(non_iso)). +:- use_module(library(dcgs)). +%:- use_module(library(types)). + +:- attribute + clpb/1, + clpb_bdd/1, + clpb_atom/1, + clpb_hash/1, + clpb_max/1, + clpb_omit_boolean/1, + clpb_visited/1. + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Compatibility predicates. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +group_pairs_by_key([], []). +group_pairs_by_key([M-N|T0], [M-[N|TN]|T]) :- + same_key(M, T0, TN, T1), + group_pairs_by_key(T1, T). + +same_key(M0, [M-N|T0], [N|TN], T) :- + M0 == M, + !, + same_key(M, T0, TN, T). +same_key(_, L, [], L). + +must_be(What, Term) :- must_be(What, unknown(Term)-1, Term). + +must_be(acyclic, Where, Term) :- !, + ( acyclic_term(Term) -> + true + ; domain_error(acyclic_term, Term, Where) + ). +must_be(list, Where, Term) :- !, + ( acyclic_term(Term), clpz_list(Term, Where) -> true + ; type_error(list, Term, Where) + ). +must_be(list(What), Where, Term) :- !, + must_be(list, Where, Term), + maplist(must_be(What, Where), Term). +must_be(ground, _, Term) :- !, + functor(Term, _, _). + +must_be(Type, Goal-Arg, Term) :- + must_be(Term, Type, Goal, Arg). + +clpz_list(Nil, _) :- Nil == []. +clpz_list(Ls, Where) :- + ( var(Ls) -> + instantiation_error(Ls, Where) + ; Ls = [_|Rest], + clpz_list(Rest, Where) + ). + + + +instantiation_error(Term) :- instantiation_error(Term, unknown(Term)-1). + +instantiation_error(_, Goal-Arg) :- + throw(error(instantiation_error, instantiation_error(Goal, Arg))). + + +domain_error(Expectation, Term) :- + domain_error(Expectation, Term, unknown(Term)-1). + +domain_error(Expectation, Term, Goal-Arg) :- + throw(error(domain_error(Expectation, Term), domain_error(Goal, Arg, Expectation, Term))). + + +type_error(Expectation, Term) :- + type_error(Expectation, Term, unknown(Term)-1). + +type_error(Expectation, Term, Goal-Arg) :- + throw(error(type_error(Expectation, Term), type_error(Goal, Arg, Expectation, Term))). + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + foldl/4 +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +foldl(Goal_3, Ls, A0, A) :- + foldl_(Ls, Goal_3, A0, A). + +foldl_([], _, A, A). +foldl_([L|Ls], G_3, A0, A) :- + call(G_3, L, A0, A1), + foldl_(Ls, G_3, A1, A). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + foldl/5 +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +foldl(Goal_4, Xs, Ys, A0, A) :- + foldl_(Xs, Ys, Goal_4, A0, A). + +foldl_([], [], _, A, A). +foldl_([X|Xs], [Y|Ys], G_4, A0, A) :- + call(G_4, X, Y, A0, A1), + foldl_(Xs, Ys, G_4, A1, A). + + +partition(Pred, Ls0, As, Bs) :- + include(Pred, Ls0, As), + exclude(Pred, Ls0, Bs). + +sum_list(Ls, S) :- sumlist(Ls, S). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Pairs. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +pairs_keys_values([], [], []). +pairs_keys_values([A-B|ABs], [A|As], [B|Bs]) :- + pairs_keys_values(ABs, As, Bs). + +pairs_keys(Ps, Ks) :- pairs_keys_values(Ps, Ks, _). + +pairs_values(Ps, Vs) :- pairs_keys_values(Ps, _, Vs). + +map_list_to_pairs(Pred, Ls, Ps) :- + map_list_to_pairs2(Ls, Pred, Ps). + +map_list_to_pairs2([], _, []). +map_list_to_pairs2([H|T0], Pred, [K-H|T]) :- + call(Pred, H, K), + map_list_to_pairs2(T0, Pred, T). + + + + + +goal_expansion(get_attr(Var, Module, Value), (var(Var),get_atts(Var, Access))) :- + Access =.. [Module,Value]. + +goal_expansion(put_attr(Var, Module, Value), put_atts(Var, Access)) :- + Access =.. [Module,Value]. + +goal_expansion(del_attr(Var, Module), (var(Var) -> put_atts(Var, -Access);true)) :- + Access =.. [Module,_]. + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Each CLP(B) variable belongs to exactly one BDD. Each CLP(B) + variable gets an attribute (in module "clpb") of the form: + + index_root(Index,Root) + + where Index is the variable's unique integer index, and Root is the + root of the BDD that the variable belongs to. + + Each CLP(B) variable also gets an attribute in module clpb_hash: an + association table node(LID,HID) -> Node, to keep the BDD reduced. + The association table of each variable must be rebuilt on occasion + to remove nodes that are no longer reachable. We rebuild the + association tables of involved variables after BDDs are merged to + build a new root. This only serves to reclaim memory: Keeping a + node in a local table even when it no longer occurs in any BDD does + not affect the solver's correctness. However, apply_shortcut/4 + relies on the invariant that every node that occurs in the relevant + BDDs is also registered in the table of its branching variable. + + A root is a logical variable with a single attribute ("clpb_bdd") + of the form: + + Sat-BDD + + where Sat is the SAT formula (in original form) that corresponds to + BDD. Sat is necessary to rebuild the BDD after variable aliasing, + and to project all remaining constraints to a list of sat/1 goals. + + Finally, a BDD is either: + + *) The integers 0 or 1, denoting false and true, respectively, or + *) A node of the form + + node(ID, Var, Low, High, Aux) + Where ID is the node's unique integer ID, Var is the + node's branching variable, and Low and High are the + node's low (Var = 0) and high (Var = 1) children. Aux + is a free variable, one for each node, that can be used + to attach attributes and store intermediate results. + + Variable aliasing is treated as a conjunction of corresponding SAT + formulae. + + You should think of CLP(B) as a potentially vast collection of BDDs + that can range from small to gigantic in size, and which can merge. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Type checking. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +is_sat(V) :- var(V), !, non_monotonic(V). +is_sat(v(V)) :- var(V), !. +is_sat(v(I)) :- integer(I), between(0, 1, I). +is_sat(I) :- integer(I), between(0, 1, I). +is_sat(A) :- atom(A). +is_sat(~A) :- is_sat(A). +is_sat(A*B) :- is_sat(A), is_sat(B). +is_sat(A+B) :- is_sat(A), is_sat(B). +is_sat(A#B) :- is_sat(A), is_sat(B). +is_sat(A=:=B) :- is_sat(A), is_sat(B). +is_sat(A=\=B) :- is_sat(A), is_sat(B). +is_sat(A==B) :- is_sat(A), is_sat(B). +is_sat(AB) :- is_sat(A), is_sat(B). +is_sat(+(Ls)) :- must_be(list, Ls), maplist(is_sat, Ls). +is_sat(*(Ls)) :- must_be(list, Ls), maplist(is_sat, Ls). +is_sat(X^F) :- var(X), is_sat(F). +is_sat(card(Is,Fs)) :- + must_be(list(ground), Is), + must_be(list, Fs), + maplist(is_sat, Fs). + +:- dynamic(monotonic/0). + +non_monotonic(X) :- + ( var_index(X, _) -> + % OK: already constrained to a CLP(B) variable + true + ; monotonic -> + instantiation_error(X) + ; true + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Rewriting to canonical expressions. + Atoms are converted to variables with a special attribute. + A global lookup table maintains the correspondence between atoms and + their variables throughout different sat/1 goals. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +% elementary +sat_rewrite(V, V) :- var(V), !. +sat_rewrite(I, I) :- integer(I), !. +sat_rewrite(A, V) :- atom(A), !, clpb_atom_var(A, V). +sat_rewrite(v(V), V). +sat_rewrite(P0*Q0, P*Q) :- sat_rewrite(P0, P), sat_rewrite(Q0, Q). +sat_rewrite(P0+Q0, P+Q) :- sat_rewrite(P0, P), sat_rewrite(Q0, Q). +sat_rewrite(P0#Q0, P#Q) :- sat_rewrite(P0, P), sat_rewrite(Q0, Q). +sat_rewrite(X^F0, X^F) :- sat_rewrite(F0, F). +sat_rewrite(card(Is,Fs0), card(Is,Fs)) :- + maplist(sat_rewrite, Fs0, Fs). +% synonyms +sat_rewrite(~P, R) :- sat_rewrite(1 # P, R). +sat_rewrite(P =:= Q, R) :- sat_rewrite(~P # Q, R). +sat_rewrite(P =\= Q, R) :- sat_rewrite(P # Q, R). +sat_rewrite(P =< Q, R) :- sat_rewrite(~P + Q, R). +sat_rewrite(P >= Q, R) :- sat_rewrite(Q =< P, R). +sat_rewrite(P < Q, R) :- sat_rewrite(~P * Q, R). +sat_rewrite(P > Q, R) :- sat_rewrite(Q < P, R). +sat_rewrite(+(Ls), R) :- foldl(or, Ls, 0, F), sat_rewrite(F, R). +sat_rewrite(*(Ls), R) :- foldl(and, Ls, 1, F), sat_rewrite(F, R). + +or(A, B, B + A). + +and(A, B, B * A). + +must_be_sat(Sat) :- + must_be(acyclic, Sat), + ( is_sat(Sat) -> true + ; no_truth_value(Sat) + ). + +no_truth_value(Term) :- domain_error(clpb_expr, Term). + +parse_sat(Sat0, Sat) :- + must_be_sat(Sat0), + sat_rewrite(Sat0, Sat), + term_variables(Sat, Vs), + maplist(enumerate_variable, Vs). + +enumerate_variable(V) :- + ( var_index_root(V, _, _) -> true + ; clpb_next_id('$clpb_next_var', Index), + put_attr(V, clpb, index_root(Index,_)), + put_empty_hash(V) + ). + +var_index(V, I) :- var_index_root(V, I, _). + +var_index_root(V, I, Root) :- get_attr(V, clpb, index_root(I,Root)). + +put_empty_hash(V) :- + empty_assoc(H0), + put_attr(V, clpb_hash, H0). + +sat_roots(Sat, Roots) :- + term_variables(Sat, Vs), + maplist(var_index_root, Vs, _, Roots0), + term_variables(Roots0, Roots). + +%% sat(+Expr) is semidet. +% +% True iff Expr is a satisfiable Boolean expression. + +sat(Sat0) :- + ( phrase(sat_ands(Sat0), Ands), Ands = [_,_|_] -> + maplist(sat, Ands) + ; parse_sat(Sat0, Sat), + sat_bdd(Sat, BDD), + sat_roots(Sat, Roots), + roots_and(Roots, Sat0-BDD, And-BDD1), + maplist(del_bdd, Roots), + maplist(=(Root), Roots), + root_put_formula_bdd(Root, And, BDD1), + is_bdd(BDD1), + satisfiable_bdd(BDD1) + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Posting many small sat/1 constraints is better than posting a huge + conjunction (or negated disjunction), because unneeded nodes are + removed from node tables after BDDs are merged. This is not + possible in sat_bdd/2 because the nodes may occur in other BDDs. A + better version of sat_bdd/2 or a proper implementation of a unique + table including garbage collection would make this obsolete and + also improve taut/2 and sat_count/2 in such cases. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +sat_ands(X) --> + ( { var(X) } -> [X] + ; { X = (A*B) } -> sat_ands(A), sat_ands(B) + ; { X = *(Ls) } -> sat_ands_(Ls) + ; { X = ~Y } -> not_ors(Y) + ; [X] + ). + +sat_ands_([]) --> []. +sat_ands_([L|Ls]) --> [L], sat_ands_(Ls). + +not_ors(X) --> + ( { var(X) } -> [~X] + ; { X = (A+B) } -> not_ors(A), not_ors(B) + ; { X = +(Ls) } -> not_ors_(Ls) + ; [~X] + ). + +not_ors_([]) --> []. +not_ors_([L|Ls]) --> [~L], not_ors_(Ls). + +del_bdd(Root) :- del_attr(Root, clpb_bdd). + +root_get_formula_bdd(Root, F, BDD) :- get_attr(Root, clpb_bdd, F-BDD). + +root_put_formula_bdd(Root, F, BDD) :- put_attr(Root, clpb_bdd, F-BDD). + +roots_and(Roots, Sat0-BDD0, Sat-BDD) :- + foldl(root_and, Roots, Sat0-BDD0, Sat-BDD), + rebuild_hashes(BDD). + +root_and(Root, Sat0-BDD0, Sat-BDD) :- + ( root_get_formula_bdd(Root, F, B) -> + Sat = F*Sat0, + bdd_and(B, BDD0, BDD) + ; Sat = Sat0, + BDD = BDD0 + ). + +bdd_and(NA, NB, And) :- + apply(*, NA, NB, And), + is_bdd(And). + +%% taut(+Expr, -T) is semidet +% +% Tautology check. Succeeds with T = 0 if the Boolean expression Expr +% cannot be satisfied, and with T = 1 if Expr is always true with +% respect to the current constraints. Fails otherwise. + +taut(Sat0, T) :- + parse_sat(Sat0, Sat), + ( T = 0, \+ sat(Sat) -> true + ; T = 1, tautology(Sat) -> true + ; false + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + The algebraic equivalence: tautology(F) <=> \+ sat(~F) does NOT + hold in CLP(B) because the quantifiers of universally quantified + variables always implicitly appear in front of the *entire* + expression. Thus we have for example: X+a is not a tautology, but + ~(X+a), meaning forall(a, ~(X+a)), is unsatisfiable: + + sat(~(X+a)) = sat(~X * ~a) = sat(~X), sat(~a) = X=0, false + + The actual negation of X+a, namely ~forall(A,X+A), in terms of + CLP(B): ~ ~exists(A, ~(X+A)), is of course satisfiable: + + ?- sat(~ ~A^ ~(X+A)). + %@ X = 0, + %@ sat(A=:=A). + + Instead, of such rewriting, we test whether the BDD of the negated + formula is 0. Critically, this avoids constraint propagation. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +tautology(Sat) :- + ( phrase(sat_ands(Sat), Ands), Ands = [_,_|_] -> + maplist(tautology, Ands) + ; catch((sat_roots(Sat, Roots), + roots_and(Roots, _-1, _-Ands), + sat_bdd(1#Sat, BDD), + bdd_and(BDD, Ands, B), + B == 0, + % reset all attributes + throw(tautology)), + tautology, + true) + ). + +satisfiable_bdd(BDD) :- + ( BDD == 0 -> false + ; BDD == 1 -> true + ; ( bdd_nodes(var_unbound, BDD, Nodes) -> + bdd_variables_classification(BDD, Nodes, Classes), + partition(var_class, Classes, Eqs, Bs, Ds), + domain_consistency(Eqs, Goal), + aliasing_consistency(Bs, Ds, Goals), + maplist(unification, [Goal|Goals]) + ; % if any variable is instantiated, we do not perform + % any propagation for now + true + ) + ). + +var_class(_=_, <). +var_class(further_branching(_,_), =). +var_class(negative_decisive(_), >). + +unification(true). +unification(A=B) :- A = B. % safe_goal/1 detects safety of this call + +var_unbound(Node) :- + node_var_low_high(Node, Var, _, _), + var(Var). + +universal_var(Var) :- get_attr(Var, clpb_atom, _). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + By aliasing consistency, we mean that all unifications X=Y, where + taut(X=:=Y, 1) holds, are posted. + + To detect this, we distinguish two kinds of variables among those + variables that are not skipped in any branch: further-branching and + negative-decisive. X is negative-decisive iff every node where X + appears as a branching variable has 0 as one of its children. X is + further-branching iff 1 is not a direct child of any node where X + appears as a branching variable. + + Any potential aliasing must involve one further-branching, and one + negative-decisive variable. X=Y must hold if, for each low branch + of nodes with X as branching variable, Y has high branch 0, and for + each high branch of nodes involving X, Y has low branch 0. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +aliasing_consistency(Bs, Ds, Goals) :- + phrase(aliasings(Bs, Ds), Goals). + +aliasings([], _) --> []. +aliasings([further_branching(B,Nodes)|Bs], Ds) --> + { var_index(B, BI) }, + aliasings_(Ds, B, BI, Nodes), + aliasings(Bs, Ds). + +aliasings_([], _, _, _) --> []. +aliasings_([negative_decisive(D)|Ds], B, BI, Nodes) --> + { var_index(D, DI) }, + ( { DI > BI, + always_false(high, DI, Nodes), + always_false(low, DI, Nodes), + var_or_atom(D, DA), var_or_atom(B, BA) } -> + [DA=BA] + ; [] + ), + aliasings_(Ds, B, BI, Nodes). + +var_or_atom(Var, VA) :- + ( get_attr(Var, clpb_atom, VA) -> true + ; VA = Var + ). + +always_false(Which, DI, Nodes) :- + phrase(nodes_always_false(Nodes, Which, DI), Opposites), + maplist(with_aux(unvisit), Opposites). + +nodes_always_false([], _, _) --> []. +nodes_always_false([Node|Nodes], Which, DI) --> + { which_node_child(Which, Node, Child), + opposite(Which, Opposite) }, + opposite_always_false(Opposite, DI, Child), + nodes_always_false(Nodes, Which, DI). + +which_node_child(low, Node, Child) :- + node_var_low_high(Node, _, Child, _). +which_node_child(high, Node, Child) :- + node_var_low_high(Node, _, _, Child). + +opposite(low, high). +opposite(high, low). + +opposite_always_false(Opposite, DI, Node) --> + ( { node_visited(Node) } -> [] + ; { node_var_low_high(Node, Var, Low, High), + with_aux(put_visited, Node), + var_index(Var, VI) }, + [Node], + ( { VI =:= DI } -> + { which_node_child(Opposite, Node, Child), + Child == 0 } + ; opposite_always_false(Opposite, DI, Low), + opposite_always_false(Opposite, DI, High) + ) + ). + +further_branching(Node) :- + node_var_low_high(Node, _, Low, High), + Low \== 1, + High \== 1. + +negative_decisive(Node) :- + node_var_low_high(Node, _, Low, High), + ( Low == 0 -> true + ; High == 0 -> true + ; false + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Instantiate all variables that only admit a single Boolean value. + + This is the case if: The variable is not skipped in any branch + leading to 1 (its being skipped means that it may be assigned + either 0 or 1 and can thus not be fixed yet), and all nodes where + it occurs as a branching variable have either lower or upper child + fixed to 0 consistently. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +domain_consistency(Eqs, Goal) :- + maplist(eq_a_b, Eqs, Vs, Values), + Goal = (Vs = Values). % propagate all assignments at once + +eq_a_b(A=B, A, B). + +consistently_false_(Which, Node) :- + which_node_child(Which, Node, Child), + Child == 0. + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + In essentially one sweep of the BDD, all variables can be classified: + Unification with 0 or 1, further branching and/or negative decisive. + + Strategy: Breadth-first traversal of the BDD, failing (and thus + clearing all attributes) if the variable is skipped in some branch, + and moving the frontier along each time. + + A formula is only satisfiable if it is a tautology after all (also + implicitly) existentially quantified variables are projected away. + However, we only need to check this explicitly if at least one + universally quantified variable appears. Otherwise, we know that + the formula is satisfiable at this point, because its BDD is not 0. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +bdd_variables_classification(BDD, Nodes, Classes) :- + nodes_variables(Nodes, Vs0), + variables_in_index_order(Vs0, Vs), + ( partition(universal_var, Vs, [_|_], Es) -> + foldl(existential, Es, BDD, 1) + ; true + ), + phrase(variables_classification(Vs, [BDD]), Classes), + maplist(with_aux(unvisit), Nodes). + +variables_classification([], _) --> []. +variables_classification([V|Vs], Nodes0) --> + { var_index(V, Index) }, + ( { phrase(nodes_with_variable(Nodes0, Index), Nodes) } -> + ( { maplist(consistently_false_(low), Nodes) } -> [V=1] + ; { maplist(consistently_false_(high), Nodes) } -> [V=0] + ; [] + ), + ( { maplist(further_branching, Nodes) } -> + [further_branching(V, Nodes)] + ; [] + ), + ( { maplist(negative_decisive, Nodes) } -> + [negative_decisive(V)] + ; [] + ), + { maplist(with_aux(unvisit), Nodes) }, + variables_classification(Vs, Nodes) + ; variables_classification(Vs, Nodes0) + ). + +nodes_with_variable([], _) --> []. +nodes_with_variable([Node|Nodes], VI) --> + { Node \== 1 }, + ( { node_visited(Node) } -> nodes_with_variable(Nodes, VI) + ; { with_aux(put_visited, Node), + node_var_low_high(Node, OVar, Low, High), + var_index(OVar, OVI) }, + { OVI =< VI }, + ( { OVI =:= VI } -> [Node] + ; nodes_with_variable([Low,High], VI) + ), + nodes_with_variable(Nodes, VI) + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Node management. Always use an existing node, if there is one. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +make_node(Var, Low, High, Node) :- + ( Low == High -> Node = Low + ; low_high_key(Low, High, Key), + ( lookup_node(Var, Key, Node) -> true + ; clpb_next_id('$clpb_next_node', ID), + Node = node(ID,Var,Low,High,_Aux), + register_node(Var, Key, Node) + ) + ). + +make_node(Var, Low, High, Node) --> + % make it conveniently usable within DCGs + { make_node(Var, Low, High, Node) }. + + +% The key of a node for hashing is determined by the IDs of its +% children. + +low_high_key(Low, High, node(LID,HID)) :- + node_id(Low, LID), + node_id(High, HID). + + +rebuild_hashes(BDD) :- + bdd_nodes(nodevar_put_empty_hash, BDD, Nodes), + maplist(re_register_node, Nodes). + +nodevar_put_empty_hash(Node) :- + node_var_low_high(Node, Var, _, _), + empty_assoc(H0), + put_attr(Var, clpb_hash, H0). + +re_register_node(Node) :- + node_var_low_high(Node, Var, Low, High), + low_high_key(Low, High, Key), + register_node(Var, Key, Node). + +register_node(Var, Key, Node) :- + get_attr(Var, clpb_hash, H0), + put_assoc(Key, H0, Node, H), + put_attr(Var, clpb_hash, H). + +lookup_node(Var, Key, Node) :- + get_attr(Var, clpb_hash, H0), + get_assoc(Key, H0, Node). + + +node_id(0, false). +node_id(1, true). +node_id(node(ID,_,_,_,_), ID). + +node_aux(Node, Aux) :- arg(5, Node, Aux). + +node_var_low_high(Node, Var, Low, High) :- + Node = node(_,Var,Low,High,_). + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + sat_bdd/2 converts a SAT formula in canonical form to an ordered + and reduced BDD. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +sat_bdd(V, Node) :- var(V), !, make_node(V, 0, 1, Node). +sat_bdd(I, I) :- integer(I), !. +sat_bdd(V^Sat, Node) :- !, sat_bdd(Sat, BDD), existential(V, BDD, Node). +sat_bdd(card(Is,Fs), Node) :- !, counter_network(Is, Fs, Node). +sat_bdd(Sat, Node) :- !, + Sat =.. [F,A,B], + sat_bdd(A, NA), + sat_bdd(B, NB), + apply(F, NA, NB, Node). + +existential(V, BDD, Node) :- + var_index(V, Index), + bdd_restriction(BDD, Index, 0, NA), + bdd_restriction(BDD, Index, 1, NB), + apply(+, NA, NB, Node). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Counter network for card(Is,Fs). +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +same_length([], []). +same_length([_|As], [_|Bs]) :- + same_length(As, Bs). + +counter_network(Cs, Fs, Node) :- + same_length([_|Fs], Indicators), + fill_indicators(Indicators, 0, Cs), + phrase(formulas_variables(Fs, Vars0), ExBDDs), + maplist(unvisit, Vars0), + % The counter network is built bottom-up, so variables with + % highest index must be processed first. + variables_in_index_order(Vars0, Vars1), + reverse(Vars1, Vars), + counter_network_(Vars, Indicators, Node0), + foldl(existential_and, ExBDDs, Node0, Node). + +% Introduce fresh variables for expressions that are not variables. +% These variables are later existentially quantified to remove them. +% Also, new variables are introduced for variables that are used more +% than once, as in card([0,1],[X,X,Y]), to keep the BDD ordered. + +formulas_variables([], []) --> []. +formulas_variables([F|Fs], [V|Vs]) --> + ( { var(F), \+ is_visited(F) } -> + { V = F, + put_visited(F) } + ; { enumerate_variable(V), + sat_rewrite(V =:= F, Sat), + sat_bdd(Sat, BDD) }, + [V-BDD] + ), + formulas_variables(Fs, Vs). + +counter_network_([], [Node], Node). +counter_network_([Var|Vars], [I|Is0], Node) :- + foldl(indicators_pairing(Var), Is0, Is, I, _), + counter_network_(Vars, Is, Node). + +indicators_pairing(Var, I, Node, Prev, I) :- make_node(Var, Prev, I, Node). + +fill_indicators([], _, _). +fill_indicators([I|Is], Index0, Cs) :- + ( memberchk(Index0, Cs) -> I = 1 + ; member(A-B, Cs), between(A, B, Index0) -> I = 1 + ; I = 0 + ), + Index1 is Index0 + 1, + fill_indicators(Is, Index1, Cs). + +existential_and(Ex-BDD, Node0, Node) :- + bdd_and(BDD, Node0, Node1), + existential(Ex, Node1, Node), + % remove attributes to avoid residual goals for variables that + % are only used as substitutes for formulas + del_attrs(Ex). + + +del_attrs(Var) :- + ( var(Var) -> + put_atts(Var, [ + -clpb(_), + -clpb_bdd(_), + -clpb_atom(_), + -clpb_hash(_), + -clpb_max(_), + -clpb_omit_boolean(_), + -clpb_visited(_) + ]) + ; true + ). + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Compute F(NA, NB). + + We use a DCG to thread through an implicit argument G0, an + association table F(IDA,IDB) -> Node, used for memoization. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +apply(F, NA, NB, Node) :- + empty_assoc(G0), + phrase(apply(F, NA, NB, Node), [G0], _). + +apply(F, NA, NB, Node) --> + ( { integer(NA), integer(NB) } -> { once(bool_op(F, NA, NB, Node)) } + ; { apply_shortcut(F, NA, NB, Node) } -> [] + ; { node_id(NA, IDA), node_id(NB, IDB), Key =.. [F,IDA,IDB] }, + ( state(G0), { get_assoc(Key, G0, Node) } -> [] + ; apply_(F, NA, NB, Node), + state(G0, G), + { put_assoc(Key, G0, Node, G) } + ) + ). + +apply_shortcut(+, NA, NB, Node) :- + ( NA == 0 -> Node = NB + ; NA == 1 -> Node = 1 + ; NB == 0 -> Node = NA + ; NB == 1 -> Node = 1 + ; false + ). + +apply_shortcut(*, NA, NB, Node) :- + ( NA == 0 -> Node = 0 + ; NA == 1 -> Node = NB + ; NB == 0 -> Node = 0 + ; NB == 1 -> Node = NA + ; false + ). + + +apply_(F, NA, NB, Node) --> + { var_less_than(NA, NB), + !, + node_var_low_high(NA, VA, LA, HA) }, + apply(F, LA, NB, Low), + apply(F, HA, NB, High), + make_node(VA, Low, High, Node). +apply_(F, NA, NB, Node) --> + { node_var_low_high(NA, VA, LA, HA), + node_var_low_high(NB, VB, LB, HB), + VA == VB }, + !, + apply(F, LA, LB, Low), + apply(F, HA, HB, High), + make_node(VA, Low, High, Node). +apply_(F, NA, NB, Node) --> % NB < NA + { node_var_low_high(NB, VB, LB, HB) }, + apply(F, NA, LB, Low), + apply(F, NA, HB, High), + make_node(VB, Low, High, Node). + + +node_varindex(Node, VI) :- + node_var_low_high(Node, V, _, _), + var_index(V, VI). + +var_less_than(NA, NB) :- + ( integer(NB) -> true + ; node_varindex(NA, VAI), + node_varindex(NB, VBI), + VAI < VBI + ). + +bool_op(+, 0, 0, 0). +bool_op(+, 0, 1, 1). +bool_op(+, 1, 0, 1). +bool_op(+, 1, 1, 1). + +bool_op(*, 0, 0, 0). +bool_op(*, 0, 1, 0). +bool_op(*, 1, 0, 0). +bool_op(*, 1, 1, 1). + +bool_op(#, 0, 0, 0). +bool_op(#, 0, 1, 1). +bool_op(#, 1, 0, 1). +bool_op(#, 1, 1, 0). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Access implicit state in DCGs. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +state(S) --> state(S, S). + +state(S0, S), [S] --> [S0]. + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Unification. X = Expr is equivalent to sat(X =:= Expr). +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +verify_attributes(Var, Other, Gs) :- + % format("~w = ~w\n", [Var,Other]), + ( get_attr(Var, clpb, index_root(I,Root)) -> + ( integer(Other) -> + ( between(0, 1, Other) -> + root_get_formula_bdd(Root, Sat, BDD0), + bdd_restriction(BDD0, I, Other, BDD), + root_put_formula_bdd(Root, Sat, BDD), + Gs = [satisfiable_bdd(BDD)] + ; no_truth_value(Other) + ) + ; atom(Other) -> + root_get_formula_bdd(Root, Sat0, _), + Gs = [root_rebuild_bdd(Root, Sat0)] + ; % due to variable aliasing, any BDDs may become unordered, + % so we need to rebuild the new BDD from the conjunction + % after the unification is in place + root_get_formula_bdd(Root, Sat0, _), + Sat = Sat0*OtherSat, + parse_sat(Other, OtherSat), + sat_roots(Sat, Roots), + phrase(formulas_(Roots), [F|Fs]), + foldl(and, Fs, F, And), + maplist(del_bdd, Roots), + maplist(=(NewRoot), Roots), + Gs = [root_rebuild_bdd(NewRoot, And)] + ) + ; Gs = [] + ). + + +formulas_([]) --> []. +formulas_([Root|Roots]) --> + ( { root_get_formula_bdd(Root, F, _) } -> + [F] + ; [] + ), + formulas_(Roots). + +root_rebuild_bdd(Root, Formula) :- + parse_sat(Formula, Sat), + sat_bdd(Sat, BDD), + is_bdd(BDD), + root_put_formula_bdd(Root, Formula, BDD), + satisfiable_bdd(BDD). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Support for project_attributes/2. + + This is called by the toplevel as + + project_attributes(+QueryVars, +AttrVars) + + in order to project all remaining constraints onto QueryVars. + + All CLP(B) variables that do not occur in QueryVars + need to be existentially quantified, so that they do not occur in + residual goals. This is very easy to do in the case of CLP(B). +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +project_attributes(QueryVars0, _) :- + include(clpb_variable, QueryVars0, QueryVars), + maplist(var_index_root, QueryVars, _, Roots0), + sort(Roots0, Roots), + maplist(remove_hidden_variables(QueryVars), Roots). + +clpb_variable(Var) :- var_index(Var, _). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + All CLP(B) variables occurring in BDDs but not in query variables + become existentially quantified. This must also be reflected in the + formula. In addition, an attribute is attached to these variables + to suppress superfluous sat(V=:=V) goals. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +remove_hidden_variables(QueryVars, Root) :- + root_get_formula_bdd(Root, Formula, BDD0), + maplist(put_visited, QueryVars), + bdd_variables(BDD0, HiddenVars0), + exclude(universal_var, HiddenVars0, HiddenVars), + maplist(unvisit, QueryVars), + foldl(existential, HiddenVars, BDD0, BDD), + foldl(quantify_existantially, HiddenVars, Formula, ExFormula), + root_put_formula_bdd(Root, ExFormula, BDD). + +quantify_existantially(E, E0, E^E0) :- put_attr(E, clpb_omit_boolean, true). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + BDD restriction. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +bdd_restriction(Node, VI, Value, Res) :- + empty_assoc(G0), + phrase(bdd_restriction_(Node, VI, Value, Res), [G0], _), + is_bdd(Res). + +bdd_restriction_(Node, VI, Value, Res) --> + ( { integer(Node) } -> { Res = Node } + ; { node_var_low_high(Node, Var, Low, High) } -> + ( { integer(Var) } -> + ( { Var =:= 0 } -> bdd_restriction_(Low, VI, Value, Res) + ; { Var =:= 1 } -> bdd_restriction_(High, VI, Value, Res) + ; { no_truth_value(Var) } + ) + ; { var_index(Var, I0), + node_id(Node, ID) }, + ( { I0 =:= VI } -> + ( { Value =:= 0 } -> { Res = Low } + ; { Value =:= 1 } -> { Res = High } + ) + ; { I0 > VI } -> { Res = Node } + ; state(G0), { get_assoc(ID, G0, Res) } -> [] + ; bdd_restriction_(Low, VI, Value, LRes), + bdd_restriction_(High, VI, Value, HRes), + make_node(Var, LRes, HRes, Res), + state(G0, G), + { put_assoc(ID, G0, Res, G) } + ) + ) + ; { domain_error(node, Node) } + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Relating a BDD to its elements (nodes and variables). + + Note that BDDs can become quite big (easily millions of nodes), and + memory space is a major bottleneck for many problems. If possible, + we therefore do not duplicate the entire BDD in memory (as in + bdd_ites/2), but only extract its features as needed. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +bdd_nodes(BDD, Ns) :- bdd_nodes(ignore_node, BDD, Ns). + +ignore_node(_). + +% VPred is a unary predicate that is called for each node that has a +% branching variable (= each inner node). + +bdd_nodes(VPred, BDD, Ns) :- + phrase(bdd_nodes_(VPred, BDD), Ns), + maplist(with_aux(unvisit), Ns). + +bdd_nodes_(VPred, Node) --> + ( { node_visited(Node) } -> [] + ; { call(VPred, Node), + with_aux(put_visited, Node), + node_var_low_high(Node, _, Low, High) }, + [Node], + bdd_nodes_(VPred, Low), + bdd_nodes_(VPred, High) + ). + +node_visited(Node) :- integer(Node). +node_visited(Node) :- with_aux(is_visited, Node). + +bdd_variables(BDD, Vs) :- + bdd_nodes(BDD, Nodes), + nodes_variables(Nodes, Vs). + +nodes_variables(Nodes, Vs) :- + phrase(nodes_variables_(Nodes), Vs), + maplist(unvisit, Vs). + +nodes_variables_([]) --> []. +nodes_variables_([Node|Nodes]) --> + { node_var_low_high(Node, Var, _, _) }, + ( { integer(Var) } -> [] + ; { is_visited(Var) } -> [] + ; { put_visited(Var) }, + [Var] + ), + nodes_variables_(Nodes). + +unvisit(V) :- del_attr(V, clpb_visited). + +is_visited(V) :- get_attr(V, clpb_visited, true). + +put_visited(V) :- put_attr(V, clpb_visited, true). + +with_aux(Pred, Node) :- + node_aux(Node, Aux), + call(Pred, Aux). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Internal consistency checks. + + To enable these checks, assert clpb:validation/0. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +:- dynamic(validation/0). + +is_bdd(BDD) :- + ( validation -> + bdd_ites(BDD, ITEs), + pairs_values(ITEs, Ls0), + sort(Ls0, Ls1), + ( same_length(Ls0, Ls1) -> true + ; domain_error(reduced_ites, (ITEs,Ls0,Ls1)) + ), + ( member(ITE, ITEs), \+ registered_node(ITE) -> + domain_error(registered_node, ITE) + ; true + ), + ( member(I, ITEs), \+ clpb_ordered(I) -> + domain_error(ordered_node, I) + ; true + ) + ; true + ). + +clpb_ordered(_-ite(Var,High,Low)) :- + ( var_index(Var, VI) -> + greater_varindex_than(High, VI), + greater_varindex_than(Low, VI) + ; true + ). + +greater_varindex_than(Node, VI) :- + ( integer(Node) -> true + ; node_var_low_high(Node, Var, _, _), + ( var_index(Var, OI) -> + OI > VI + ; true + ) + ). + +registered_node(Node-ite(Var,High,Low)) :- + ( var(Var) -> + low_high_key(Low, High, Key), + lookup_node(Var, Key, Node0), + Node == Node0 + ; true + ). + +bdd_ites(BDD, ITEs) :- + bdd_nodes(BDD, Nodes), + maplist(node_ite, Nodes, ITEs). + +node_ite(Node, Node-ite(Var,High,Low)) :- + node_var_low_high(Node, Var, Low, High). + +%% labeling(+Vs) is multi. +% +% Enumerate concrete solutions. Assigns truth values to the Boolean +% variables Vs such that all stated constraints are satisfied. + +labeling(Vs0) :- + must_be(list, Vs0), + maplist(labeling_var, Vs0), + variables_in_index_order(Vs0, Vs), + maplist(indomain, Vs). + +labeling_var(V) :- var(V), !. +labeling_var(V) :- V == 0, !. +labeling_var(V) :- V == 1, !. +labeling_var(V) :- domain_error(clpb_variable, V). + +variables_in_index_order(Vs0, Vs) :- + maplist(var_with_index, Vs0, IVs0), + keysort(IVs0, IVs), + pairs_values(IVs, Vs). + +var_with_index(V, I-V) :- + ( var_index_root(V, I, _) -> true + ; I = 0 + ). + +indomain(0). +indomain(1). + + +%% sat_count(+Expr, -Count) is det. +% +% Count the number of admissible assignments. Count is the number of +% different assignments of truth values to the variables in the +% Boolean expression Expr, such that Expr is true and all posted +% constraints are satisfiable. +% +% A common form of invocation is `sat_count(+[1|Vs], Count)`: This +% counts the number of admissible assignments to `Vs` without imposing +% any further constraints. +% +% Examples: +% +% == +% ?- sat(A =< B), Vs = [A,B], sat_count(+[1|Vs], Count). +% Vs = [A, B], +% Count = 3, +% sat(A=:=A*B). +% +% ?- length(Vs, 120), +% sat_count(+Vs, CountOr), +% sat_count(*(Vs), CountAnd). +% Vs = [...], +% CountOr = 1329227995784915872903807060280344575, +% CountAnd = 1. +% == + +sat_count(Sat0, N) :- + catch((parse_sat(Sat0, Sat), + sat_bdd(Sat, BDD), + sat_roots(Sat, Roots), + roots_and(Roots, _-BDD, _-BDD1), + % we mark variables that occur in Sat0 as visited ... + term_variables(Sat0, Vs), + maplist(put_visited, Vs), + % ... so that they do not appear in Vs1 ... + bdd_variables(BDD1, Vs1), + partition(universal_var, Vs1, Univs, Exis), + % ... and then remove remaining variables: + foldl(universal, Univs, BDD1, BDD2), + foldl(existential, Exis, BDD2, BDD3), + variables_in_index_order(Vs, IVs), + foldl(renumber_variable, IVs, 1, VNum), + bdd_count(BDD3, VNum, Count0), + var_u(BDD3, VNum, P), + % Do not unify N directly, because we are not prepared + % for propagation here in case N is a CLP(B) variable. + N0 is 2^(P - 1)*Count0, + % reset all attributes and Aux variables + throw(count(N0))), + count(N0), + N = N0). + +universal(V, BDD, Node) :- + var_index(V, Index), + bdd_restriction(BDD, Index, 0, NA), + bdd_restriction(BDD, Index, 1, NB), + apply(*, NA, NB, Node). + +renumber_variable(V, I0, I) :- + put_attr(V, clpb, index_root(I0,_)), + I is I0 + 1. + +bdd_count(Node, VNum, Count) :- + ( integer(Node) -> Count = Node + ; node_aux(Node, Count), + ( integer(Count) -> true + ; node_var_low_high(Node, V, Low, High), + bdd_count(Low, VNum, LCount), + bdd_count(High, VNum, HCount), + bdd_pow(Low, V, VNum, LPow), + bdd_pow(High, V, VNum, HPow), + Count0 is LPow*LCount + HPow*HCount, + Count = Count0 + ) + ). + + +bdd_pow(Node, V, VNum, Pow) :- + var_index(V, Index), + var_u(Node, VNum, P), + Pow is 2^(P - Index - 1). + +var_u(Node, VNum, Index) :- + ( integer(Node) -> Index = VNum + ; node_varindex(Node, Index) + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Pick a solution in such a way that each solution is equally likely. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +%% random_labeling(+Seed, +Vs) is det. +% +% Select a single random solution. An admissible assignment of truth +% values to the Boolean variables in Vs is chosen in such a way that +% each admissible assignment is equally likely. Seed is an integer, +% used as the initial seed for the random number generator. + +single_bdd(Vars0) :- + maplist(monotonic_variable, Vars0, Vars), + % capture all variables with a single BDD + sat(+[1|Vars]). + +random_labeling(Seed, Vars) :- + must_be(list, Vars), + set_random(seed(Seed)), + ( ground(Vars) -> true + ; catch((single_bdd(Vars), + once((member(Var, Vars),var(Var))), + var_index_root(Var, _, Root), + root_get_formula_bdd(Root, _, BDD), + bdd_variables(BDD, Vs), + variables_in_index_order(Vs, IVs), + foldl(renumber_variable, IVs, 1, VNum), + phrase(random_bindings(VNum, BDD), Bs), + maplist(del_attrs, Vs), + % reset all attribute modifications + throw(randsol(Vars, Bs))), + randsol(Vars, Bs), + true), + maplist(call, Bs), + % set remaining variables to 0 or 1 with equal probability + include(var, Vars, Remaining), + maplist(maybe_zero, Remaining) + ). + +maybe_zero(Var) :- + ( maybe -> Var = 0 + ; Var = 1 + ). + +random_bindings(_, Node) --> { Node == 1 }, !. +random_bindings(VNum, Node) --> + { node_var_low_high(Node, Var, Low, High), + bdd_count(Node, VNum, Total), + bdd_count(Low, VNum, LCount) }, + ( { maybe(LCount, Total) } -> + [Var=0], random_bindings(VNum, Low) + ; [Var=1], random_bindings(VNum, High) + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Find solutions with maximum weight. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +%% weighted_maximum(+Weights, +Vs, -Maximum) is multi. +% +% Enumerate weighted optima over admissible assignments. Maximize a +% linear objective function over Boolean variables Vs with integer +% coefficients Weights. This predicate assigns 0 and 1 to the +% variables in Vs such that all stated constraints are satisfied, and +% Maximum is the maximum of sum(Weight_i*V_i) over all admissible +% assignments. On backtracking, all admissible assignments that +% attain the optimum are generated. +% +% This predicate can also be used to _minimize_ a linear Boolean +% program, since negative integers can appear in Weights. +% +% Example: +% +% == +% ?- sat(A#B), weighted_maximum([1,2,1], [A,B,C], Maximum). +% A = 0, B = 1, C = 1, Maximum = 3. +% == + +weighted_maximum(Ws, Vars, Max) :- + must_be(list(integer), Ws), + must_be(list(var), Vars), + single_bdd(Vars), + Vars = [Var|_], + var_index_root(Var, _, Root), + root_get_formula_bdd(Root, _, BDD0), + bdd_variables(BDD0, Vs), + % existentially quantify variables that are not considered + maplist(put_visited, Vars), + exclude(is_visited, Vs, Unvisited), + maplist(unvisit, Vars), + foldl(existential, Unvisited, BDD0, BDD), + maplist(var_with_index, Vars, IVs), + pairs_keys_values(Pairs0, IVs, Ws), + keysort(Pairs0, Pairs1), + pairs_keys_values(Pairs1, IVs1, WeightsIndexOrder), + pairs_values(IVs1, VarsIndexOrder), + % Pairs is a list of Var-Weight terms, in index order of Vars + pairs_keys_values(Pairs, VarsIndexOrder, WeightsIndexOrder), + bdd_maximum(BDD, Pairs, Max), + max_labeling(BDD, Pairs). + +max_labeling(1, Pairs) :- max_upto(Pairs, _, _). +max_labeling(node(_,Var,Low,High,Aux), Pairs0) :- + max_upto(Pairs0, Var, Pairs), + get_attr(Aux, clpb_max, max(_,Dir)), + direction_labeling(Dir, Var, Low, High, Pairs). + +max_upto([], _, _). +max_upto([Var0-Weight|VWs0], Var, VWs) :- + ( Var == Var0 -> VWs = VWs0 + ; Weight =:= 0 -> + ( Var0 = 0 ; Var0 = 1 ), + max_upto(VWs0, Var, VWs) + ; Weight < 0 -> Var0 = 0, max_upto(VWs0, Var, VWs) + ; Var0 = 1, max_upto(VWs0, Var, VWs) + ). + +direction_labeling(low, 0, Low, _, Pairs) :- max_labeling(Low, Pairs). +direction_labeling(high, 1, _, High, Pairs) :- max_labeling(High, Pairs). + +bdd_maximum(1, Pairs, Max) :- + pairs_values(Pairs, Weights0), + include(<(0), Weights0, Weights), + sum_list(Weights, Max). +bdd_maximum(node(_,Var,Low,High,Aux), Pairs0, Max) :- + ( get_attr(Aux, clpb_max, max(Max,_)) -> true + ; ( skip_to_var(Var, Weight, Pairs0, Pairs), + ( Low == 0 -> + bdd_maximum_(High, Pairs, MaxHigh, MaxToHigh), + Max is MaxToHigh + MaxHigh + Weight, + Dir = high + ; High == 0 -> + bdd_maximum_(Low, Pairs, MaxLow, MaxToLow), + Max is MaxToLow + MaxLow, + Dir = low + ; bdd_maximum_(Low, Pairs, MaxLow, MaxToLow), + bdd_maximum_(High, Pairs, MaxHigh, MaxToHigh), + Max0 is MaxToLow + MaxLow, + Max1 is MaxToHigh + MaxHigh + Weight, + Max is max(Max0,Max1), + ( Max0 =:= Max1 -> Dir = _Any + ; Max0 < Max1 -> Dir = high + ; Dir = low + ) + ), + store_maximum(Aux, Max, Dir) + ) + ). + +bdd_maximum_(Node, Pairs, Max, MaxTo) :- + bdd_maximum(Node, Pairs, Max), + between_weights(Node, Pairs, MaxTo). + +store_maximum(Aux, Max, Dir) :- put_attr(Aux, clpb_max, max(Max,Dir)). + +between_weights(Node, Pairs0, MaxTo) :- + ( Node == 1 -> MaxTo = 0 + ; node_var_low_high(Node, Var, _, _), + phrase(skip_to_var_(Var, _, Pairs0, _), Weights0), + include(<(0), Weights0, Weights), + sum_list(Weights, MaxTo) + ). + +skip_to_var(Var, Weight, Pairs0, Pairs) :- + phrase(skip_to_var_(Var, Weight, Pairs0, Pairs), _). + +skip_to_var_(Var, Weight, [Var0-Weight0|VWs0], VWs) --> + ( { Var == Var0 } -> + { Weight = Weight0, VWs0 = VWs } + ; ( { Weight0 =< 0 } -> [] + ; [Weight0] + ), + skip_to_var_(Var, Weight, VWs0, VWs) + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Projection to residual goals. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + +:- dynamic(clpb_residuals/1). + +attribute_goals(Var) --> + { var_index_root(Var, _, Root) }, + ( { root_get_formula_bdd(Root, Formula, BDD) } -> + { del_bdd(Root) }, + ( { clpb_residuals(bdd) } -> + { bdd_nodes(BDD, Nodes), + phrase(nodes(Nodes), Ns) }, + [clpb:'$clpb_bdd'(Ns)] + ; { prepare_global_variables(BDD), + phrase(sat_ands(Formula), Ands0), + ands_fusion(Ands0, Ands), + maplist(formula_anf, Ands, ANFs0), + sort(ANFs0, ANFs1), + exclude(eq_1, ANFs1, ANFs2), + variables_separation(ANFs2, ANFs) }, + sats(ANFs) + ), + ( { get_attr(Var, clpb_atom, Atom) } -> + [clpb:sat(Var=:=Atom)] + ; [] + ), + % formula variables not occurring in the BDD should be booleans + { bdd_variables(BDD, Vs), + maplist(del_clpb, Vs), + term_variables(Formula, RestVs0), + include(clpb_variable, RestVs0, RestVs) }, + booleans(RestVs) + ; boolean(Var) % the variable may have occurred only in taut/2 + ). + +del_clpb(Var) :- + del_attr(Var, clpb), + del_attr(Var, clpb_hash). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + To make residual projection work with recorded constraints, the + global counters must be adjusted so that new variables and nodes + also get new IDs. Also, clpb_next_id/2 is used to actually create + these counters, because creating them with b_setval/2 would make + them [] on backtracking, which is quite unfortunate in itself. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +b_setval(K, T) :- bb_b_put(K, T). +nb_setval(K, T) :- bb_put(K, T). +b_getval(K, T) :- bb_get(K, T). + +prepare_global_variables(BDD) :- + clpb_next_id('$clpb_next_var', V0), + clpb_next_id('$clpb_next_node', N0), + bdd_nodes(BDD, Nodes), + foldl(max_variable_node, Nodes, V0-N0, MaxV0-MaxN0), + MaxV is MaxV0 + 1, + MaxN is MaxN0 + 1, + b_setval('$clpb_next_var', MaxV), + b_setval('$clpb_next_node', MaxN). + +max_variable_node(Node, V0-N0, V-N) :- + node_id(Node, N1), + node_varindex(Node, V1), + N is max(N0,N1), + V is max(V0,V1). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Fuse formulas that share the same variables into single conjunctions. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +ands_fusion(Ands0, Ands) :- + maplist(with_variables, Ands0, Pairs0), + keysort(Pairs0, Pairs), + group_pairs_by_key(Pairs, Groups), + pairs_values(Groups, Andss), + maplist(list_to_conjunction, Andss, Ands). + +with_variables(F, Vs-F) :- + term_variables(F, Vs0), + variables_in_index_order(Vs0, Vs). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + If possible, separate variables into different sat/1 goals. + A formula F can be split in two if for two of its variables A and B, + taut((A^F)*(B^F) =:= F, 1) holds. In the first conjunct, A does not + occur, and in the second, B does not occur. We separate variables + until that is no longer possible. There may be a better way to do this. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +variables_separation(Fs0, Fs) :- separation_fixpoint(Fs0, [], Fs). + +separation_fixpoint(Fs0, Ds0, Fs) :- + phrase(variables_separation_(Fs0, Ds0, Rest), Fs1), + partition(anf_done, Fs1, Ds1, Fs2), + maplist(arg(1), Ds1, Ds2), + maplist(arg(1), Fs2, Fs3), + append(Ds0, Ds2, Ds3), + append(Rest, Fs3, Fs4), + sort(Fs4, Fs5), + sort(Ds3, Ds4), + ( Fs5 == [] -> Fs = Ds4 + ; separation_fixpoint(Fs5, Ds4, Fs) + ). + +anf_done(done(_)). + +variables_separation_([], _, []) --> []. +variables_separation_([F0|Fs0], Ds, Rest) --> + ( { member(Done, Ds), F0 == Done } -> + variables_separation_(Fs0, Ds, Rest) + ; { sat_rewrite(F0, F), + sat_bdd(F, BDD), + bdd_variables(BDD, Vs0), + exclude(universal_var, Vs0, Vs), + maplist(existential_(BDD), Vs, Nodes), + phrase(pairs(Nodes), Pairs), + group_pairs_by_key(Pairs, Groups), + phrase(groups_separation(Groups, BDD), ANFs) }, + ( { ANFs = [_|_] } -> + list(ANFs), + { Rest = Fs0 } + ; [done(F0)], + variables_separation_(Fs0, Ds, Rest) + ) + ). + + +existential_(BDD, V, Node) :- existential(V, BDD, Node). + +groups_separation([], _) --> []. +groups_separation([BDD1-BDDs|Groups], OrigBDD) --> + { phrase(separate_pairs(BDDs, BDD1, OrigBDD), Nodes) }, + ( { Nodes = [_|_] } -> + nodes_anfs([BDD1|Nodes]) + ; [] + ), + groups_separation(Groups, OrigBDD). + +separate_pairs([], _, _) --> []. +separate_pairs([BDD2|Ps], BDD1, OrigBDD) --> + ( { apply(*, BDD1, BDD2, And), + And == OrigBDD } -> + [BDD2] + ; [] + ), + separate_pairs(Ps, BDD1, OrigBDD). + +nodes_anfs([]) --> []. +nodes_anfs([N|Ns]) --> { node_anf(N, ANF) }, [anf(ANF)], nodes_anfs(Ns). + +pairs([]) --> []. +pairs([V|Vs]) --> pairs_(Vs, V), pairs(Vs). + +pairs_([], _) --> []. +pairs_([B|Bs], A) --> [A-B], pairs_(Bs, A). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Set the Prolog flag clpb_residuals to bdd to obtain the BDD nodes + as residuals. Note that they cannot be used as regular goals. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +nodes([]) --> []. +nodes([Node|Nodes]) --> + { node_var_low_high(Node, Var0, Low, High), + var_or_atom(Var0, Var), + maplist(node_projection, [Node,High,Low], [ID,HID,LID]), + var_index(Var0, VI) }, + [ID-(v(Var,VI) -> HID ; LID)], + nodes(Nodes). + + +node_projection(Node, Projection) :- + node_id(Node, ID), + ( integer(ID) -> Projection = node(ID) + ; Projection = ID + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + By default, residual goals are sat/1 calls of the remaining formulas, + using (mostly) algebraic normal form. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +sats([]) --> []. +sats([A|As]) --> [clpb:sat(A)], sats(As). + +booleans([]) --> []. +booleans([B|Bs]) --> boolean(B), { del_clpb(B) }, booleans(Bs). + +boolean(Var) --> + ( { get_attr(Var, clpb_omit_boolean, true) } -> [] + ; [clpb:sat(Var =:= Var)] + ). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Relate a formula to its algebraic normal form (ANF). +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +formula_anf(Formula0, ANF) :- + parse_sat(Formula0, Formula), + sat_bdd(Formula, Node), + node_anf(Node, ANF). + +node_anf(Node, ANF) :- + node_xors(Node, Xors0), + maplist(maplist(monotonic_variable), Xors0, Xors), + maplist(list_to_conjunction, Xors, Conjs), + ( Conjs = [Var,C|Rest], clpb_var(Var) -> + foldl(xor, Rest, C, RANF), + ANF = (Var =\= RANF) + ; Conjs = [One,Var,C|Rest], One == 1, clpb_var(Var) -> + foldl(xor, Rest, C, RANF), + ANF = (Var =:= RANF) + ; Conjs = [C|Cs], + foldl(xor, Cs, C, ANF) + ). + +monotonic_variable(Var0, Var) :- + ( var(Var0), monotonic -> + Var = v(Var0) + ; Var = Var0 + ). + +clpb_var(Var) :- var(Var), !. +clpb_var(v(_)). + +list_to_conjunction([], 1). +list_to_conjunction([L|Ls], Conj) :- foldl(and, Ls, L, Conj). + +xor(A, B, B # A). + +eq_1(V) :- V == 1. + +node_xors(Node, Xors) :- + phrase(xors(Node), Xors0), + % we remove elements that occur an even number of times (A#A --> 0) + maplist(sort, Xors0, Xors1), + pairs_keys_values(Pairs0, Xors1, _), + keysort(Pairs0, Pairs), + group_pairs_by_key(Pairs, Groups), + exclude(even_occurrences, Groups, Odds), + pairs_keys(Odds, Xors2), + maplist(exclude(eq_1), Xors2, Xors). + +even_occurrences(_-Ls) :- length(Ls, L), L mod 2 =:= 0. + +xors(Node) --> + ( { Node == 0 } -> [] + ; { Node == 1 } -> [[1]] + ; { node_var_low_high(Node, Var0, Low, High), + var_or_atom(Var0, Var), + node_xors(Low, Ls0), + node_xors(High, Hs0), + maplist(with_var(Var), Ls0, Ls), + maplist(with_var(Var), Hs0, Hs) }, + list(Ls0), + list(Ls), + list(Hs) + ). + +list([]) --> []. +list([L|Ls]) --> [L], list(Ls). + +with_var(Var, Ls, [Var|Ls]). + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Global variables for unique node and variable IDs and atoms. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + +make_clpb_var('$clpb_next_var') :- nb_setval('$clpb_next_var', 0). + +make_clpb_var('$clpb_next_node') :- nb_setval('$clpb_next_node', 0). + +make_clpb_var('$clpb_atoms') :- + empty_assoc(E), + nb_setval('$clpb_atoms', E). + +:- initialization((make_clpb_var(_),false;true)). + +clpb_next_id(Var, ID) :- + b_getval(Var, ID), + Next is ID + 1, + b_setval(Var, Next). + +clpb_atom_var(Atom, Var) :- + b_getval('$clpb_atoms', A0), + ( get_assoc(Atom, A0, Var) -> true + ; put_attr(Var, clpb_atom, Atom), + put_attr(Var, clpb_omit_boolean, true), + put_assoc(Atom, A0, Var, A), + b_setval('$clpb_atoms', A) + ). + + + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +Compatibility predicates. +- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + +include(Goal, List, Is) :- + include_(List, Goal, Is). + +include_([], _, []). +include_([X1|Xs1], P, Is) :- + ( call(P, X1) + -> Is = [X1|Is1] + ; Is = Is1 + ), + include_(Xs1, P, Is1). + + +exclude(Goal, List, Is) :- + exclude_(List, Goal, Is). + +exclude_([], _, []). +exclude_([X1|Xs1], P, Is) :- + ( call(P, X1) + -> Is = Is1 + ; Is = [X1|Is1] + ), + exclude_(Xs1, P, Is1). + + +partition(Pred, List, Less, Equal, Greater) :- + partition_(List, Pred, Less, Equal, Greater). + +partition_([], _, [], [], []). +partition_([H|T], Pred, L, E, G) :- + call(Pred, H, Diff), + partition_(Diff, H, Pred, T, L, E, G). + +partition_(<, H, Pred, T, [H|Rest], E, G) :- + partition_(T, Pred, Rest, E, G). +partition_(=, H, Pred, T, L, [H|Rest], G) :- + partition_(T, Pred, L, Rest, G). +partition_(>, H, Pred, T, L, E, [H|Rest]) :- + partition_(T, Pred, L, E, Rest).