Fast Set Bounds Propagation Using a BDD-SAT Hybrid

Binary Decision Diagram (BDD) based set bounds propagation is a powerful approach to solving set-constraint satisfaction problems. However, prior BDD based techniques incur the signiﬁcant overhead of constructing and manipulating graphs during search. We present a set-constraint solver which combines BDD-based set-bounds propagators with the learning abilities of a modern SAT solver. Together with a number of improvements beyond the basic algorithm, this solver is highly competitive with existing propagation based set constraint solvers.


Introduction
It is often convenient to model a constraint satisfaction problem (CSP) using finite set variables and set relationships between them.A common approach to solving finite domain CSPs is using a combination of backtracking search and a constraint propagation algorithm.The propagation algorithm attempts to enforce consistency on the values in the domains of the constraint variables by removing values from the domains of variables that cannot form part of a complete solution to the system of constraints.The most common level of consistency is set bounds consistency (Gervet, 1997) where the solver keeps track for each set of which elements are definitely in or out of the set.Many solvers use set bounds consistency including ECLiPSe (IC-PARC, 2003), Gecode (GECODE, 2008), and ILOG SOLVER (ILOG, 2004).
Set bounds propagation is supported by solvers since stronger notions of propagation such as domain propagation require representing exponentially large domains of possible values.However, Lagoon and Stuckey (2004) demonstrated that it is possible to use reduced ordered binary decision diagrams (BDDs) as a compact representation of both set domains and of set constraints, thus permitting set domain propagation.A domain propagator ensures that every value in the domain of a set variable can be extended to a complete assignment of all of the variables in a constraint.The use of the BDD representation comes with several additional benefits.The ability to easily conjoin and existentially quantify BDDs allows the removal of intermediate variables, thus strengthening propagation, and also makes the construction of propagators for global constraints straightforward.
Given the natural way in which BDDs can be used to model set constraint problems, it is therefore worthwhile utilising BDDs to construct other types of set solver.Indeed it has been previously demonstrated (Hawkins, Lagoon, & Stuckey, 2004, 2005) that set bounds propagation can be efficiently implemented using BDDs to represent constraints and domains of variables.A major benefit of the BDD-based approach is that it frees us from the need to laboriously construct set bounds propagators for each new constraint by hand.Moreover, correctness and optimality of such BDD-based propagators follow by construction.The other advantages of the BDD-based representation identified above still apply, and the resulting solver performs very favourably when compared with existing set bounds solvers.
But set bounds propagation using BDDs still constructs BDDs during propagation, which is a considerable overhead.In this paper we show how we can perform BDD-based set bounds propagation using a marking algorithm that perform linear scans of the BDD representation of the constraint without constructing new BDDs.The resulting set bounds propagators are substantially faster than those using BDDs.
The contributions of this paper are: • Efficient set bounds propagators: No new BDDs are constructed during propagation, so it is very fast.
• Graph reuse: We can reuse a single BDD for multiple copies of the same constraint, and hence handle larger problems.
• Ordering flexibility: We are not restricted to a single global ordering of Booleans for constructing BDDs.
• Filtering: We can keep track of which parts of the set variable can really make a difference, and reduce the amount of propagation.
Pure set-bounds propagation tends to perform badly, however, in problems where a large number of similar regions of the search space must be explored.We therefore embed the set-bounds propagators in MiniSAT (Eén & Sörensson, 2003), to provide SAT-style clause learning.
In the next section, we introduce propagation-based solving for set problems, and briefly discuss SAT solving.In Section 3 we discuss binary decision diagrams (BDDs) and how to implement set bounds propagation using BDDs.Then in Section 4, we present the propagation algorithm used by the hybrid solver, together with a number of variations upon the standard algorithm.In Section 5, we show how to incorporate reason generation with BDD propagation to build a hybrid solver In Section 6 we test the performance of the solver on a variety of set-constraint problems, and compare with other set-constraint solvers.In Section 7 we discuss related work, before concluding in Section 8.

Propagation-based Solving
Propagation based approaches to solving set constraint problems represent the problem using a domain storing the possible values of each set variable, and propagators for each constraint, that remove values from the domain of a variable that are inconsistent with values for other variables.Propagation is combined with backtracking search to find solutions.
A domain D is a complete mapping from the fixed finite set of variables V to finite collections of finite sets of integers.The domain of a variable v is the set D(v).A domain D 1 is said to be stronger than a domain D 2 , written for all variables v ∈ V.A domain D can be interpreted as the constraint v∈V v ∈ D(v).
For set constraints we will often be interested in restricting variables to take on convex domains.A set of sets K is convex if a, b ∈ K and a ⊆ c ⊆ b implies c ∈ K.We use interval notation [a, b] where a ⊆ b to represent the (minimal) convex set K including a and b.For any finite collection of sets K = {a 1 , a 2 , . . ., a n }, we define the convex closure of K: conv (K) = [∩ a∈K a, ∪ a∈K a].We extend the concept of convex closure to domains by defining ran(D) to be the domain such that ran(D)(v) = conv (D(v)) for all v ∈ V.
A valuation θ is a set of mappings from the set of variables V to sets of integer values, written {v 1 → d 1 , . . ., v n → d n }.A valuation can be extended to apply to constraints involving the variables in the obvious way.Let vars be the function that returns the set of variables appearing in an expression, constraint or valuation.In an abuse of notation, we say a valuation is an element of a domain D,

Constraints, Propagators and Propagation
A constraint is a restriction placed on the allowable values for a set of variables.We shall use primitive set constraints such as (membership order) u < v, where u, v, w are set variables, k is an integer.We can also construct more complicated constraints which are (possibly existentially quantified) conjunctions of primitive set constraints.We define the solutions of a constraint c to be the set of valuations θ on vars(c) that make the constraint true.
We associate a propagator with every constraint.A propagator f is a monotonically decreasing function from domains to domains, so A propagation solver solv (F, D) for a set of propagators F and a domain D repeatedly applies the propagators in F starting from the domain D until a fixpoint is reached.
Example 1 A small example of a set-constraint problem would be to, given a universe consisting of the elements {1, 2, 3, 4}, find values for variables x, y, z such that z

Set Bounds Consistency
A domain D is (set) bounds consistent for a constraint c if for every variable v ∈ vars(c) the upper bound of D(v) is the union of the values of v in all solutions of c in D, and the lower bound of D(v) is the intersection of the values of v in all solutions of c in D. We define the set bounds propagator for a constraint c as Then sb(c)(D) is always bounds consistent with c.

Boolean Satisfiability (SAT)
Boolean Satisfiability or SAT solvers are a special case of propagation-based solvers, restricted to Boolean variables and clause constraints.
The Davis-Putnam-Logemann-Loveland (DPLL) algorithm (Davis, Logemann, & Loveland, 1962), on which most modern SAT solvers are based, is also a propagation-based approach to solving SAT problems.It interleaves two phases -search, where an unfixed variable is assigned a value, and propagation (so called unit propagation).
Modern SAT solvers incorporate sophisticated engineering to propagate constraints very fast, to record as nogoods part of the search that lead to failure, and to automate the search by keeping track of how often a variable is part of the reason for causing failure (activity) and concentrating search on variables with high activity.Modern SAT solvers also frequently restart the search from scratch relying on nogoods recording to prevent repeated search, and activity to drive the search into more profitable areas.See e.g., the report by Eén and Sörensson (2003) for a good introduction to modern SAT solving.
A rough architecture of a modern SAT solver is illustrated in Figure 1.Search starts the unit propagation process which interacts with the clause database and may detect failure, which initiates conflict analysis.Unit propagation records for each literal that is made true, the clause that explains why the literal become true.Conflict analysis uses the graph of explanations to construct a nogood which is a resolvent of clauses causing the failure that adds to the strength of unit propagation.This is stored in the clause database and causes search to backjump.It prevents the search revisiting the same set of decisions.Not detailed here are activity counters which record which variables are most responsible for failure, these are the variables chosen for labelling by the search.
Reduced Ordered Binary Decision Diagrams are a well-known method of representing Boolean functions on Boolean variables using directed acyclic graphs with a single root.Every internal node n(v, f, t) in a BDD r is labelled with a Boolean variable v ∈ B, and has two outgoing arcs -the 'false' arc (to BDD f ) and the 'true' arc (to BDD t).Leaf nodes are either F (false) or T (true).Each node represents a single test of the labelled variable; when traversing the tree the appropriate arc is followed depending on the value of ) is shown as a circle labelled v with a dotted arc to the f BDD, and a solid arc to the t BDD.
the variable.Define the size |r| as the number of internal nodes in a BDD r, and VAR(r) as the set of variables v ∈ B appearing in some internal node in r.
Reduced Ordered Binary Decision Diagrams (BDDs) (Bryant, 1986) require that the BDD is: reduced, that is it contains no identical nodes (nodes with the same variable label and identical true and false arcs) and has no redundant tests (no node has both true and false arcs leading to the same node); and ordered, if there is an arc from a node labelled v 1 to a node labelled v 2 then v 1 ≺ v 2 .A BDD has the nice property that the function representation is canonical up to variable reordering.This permits efficient implementations of many Boolean operations.
A Boolean variable v is said to be fixed in a BDD r if either for every node n(v, f, t) ∈ r t is the constant F node, or for every node n(v, f, t) f is the constant F node.Such variables can be identified in a linear time scan over the domain BDD (see e.g., Hawkins et al., 2005).For convenience, if φ is a BDD, we write φ to denote the BDD representing the conjunction of the fixed variables of φ.
Example 3 Figure 2(a) gives an example of a BDD representing the formula v 3 ∧ ¬v 4 ∧ ¬v 5 ∧v 6 ∧v 7 .Figure 2(b) gives an example of a more complex BDD representing the formula x 1 + x 2 + x 3 + x 4 + x 5 ≤ 2 where we interpret the Booleans as 0-1 integers.One can verify that the valuation {x 1 → 1, x 2 → 0, x 3 → 1, x 4 → 0, x 5 → 0} makes the formula true by following the path right, left, right, left, left from the root.

Set Propagation using BDDs
The key step in building set propagation using BDDs is to realize that we can represent a finite set domain using a BDD.

Representing domains
If v is a set variable ranging over subsets of {1, . . ., N }, then we can represent v using the Boolean variables V (v) = {v 1 , . . ., v N } ⊆ B, where v i is true iff i ∈ v.We will order the We can represent a valuation θ using a formula Then a domain of variable v, D(v) can be represented as φ = a∈D(v) R({v → a}).This formula can be represented by a BDD.The set bounds of v can be obtained by extracting the fixed variables from this BDD, φ .
For example the valuation θ of Example 1 is represented by the formula R(θ): And the domain , 6, 7}, {1, 2, 3, 6, 7, 8, 9}] is represented by the BDD in Figure 2(a) since v 3 , v 6 and v 7 are true so 3, 6, 7 are definitely in the set, and v 4 and v 5 are false so 4 and 5 are definitely not in the set.

Representing constraints
We can similarly model any set constraint c as a BDD B(c) using the Boolean variable representation V (v) of its set variables v.By ordering the variables in each BDD carefully we can build small representations of the formulae.The pointwise order of Boolean variables is defined as follows.Given set variables u ≺ v ≺ w ranging over sets from {1, . . ., N } we order the Boolean variables as The representation B(c) is simply ∨ θ∈solns(c) R(θ).For primitive set constraints (using the pointwise order) this size is linear in N .For more details see the work of Hawkins et al. (2005).The BDD representation of |x| ≤ 2 is shown in Figure 2(b), for N = 5.

BDD-based Set Bounds Propagation
We can build a set bounds propagator, more or less from the definition, since we have BDDs to represent domains and constraints.
We simply conjoin the domains to the constraint obtaining φ, then extract the fixed variables from the result, and then project out the relevant part for each variable v.The set bounds propagation can be improved by removing the fixed variables as soon as possible.The improved definition is given by Hawkins et al. (2004).Overall the complexity can be made O(|B(c)|).
The updated set bounds can be used to simplify the BDD representing the propagator.Since fixed variables will never interact further with propagation they can be projected out of B(c), so we can replace B(c) by

Tseitin Transformation
It is possible to convert any Boolean circuit to a pure SAT representation; the method for doing so is generally attributed to Tseitin (1968).Figure 3 gives pseudo code for the translation of a BDD rooted at node, returning a pair of (Boolean variable, set of clauses).The clauses enforce that the Boolean variable takes the truth value of the BDD.Like most BDD algorithms it relies on marking the visited nodes to ensure each node is visited at most once.It assumes the array visit[] is initially all bottom ⊥, and on first visiting a node stores the corresponding Boolean variable in visit[].A more comprehensive discussion of the Tseitin transformation is presented by Eén and Sörensson (2006).
The constraint is enforced by fixing the variable corresponding to the root node to true.An advantage of replacing a BDD by its Tseitin representation is that we can use an unmodified SAT solver to then tackle BDD-based set constraint problems.We shall see in Section 6 that this approach cannot compete with handling the BDDs directly.

Faster Set-bounds Propagation
While set bounds propagation using BDDs is much faster than set domain propagation and often better than set domain propagation (or other variations of propagation for sets) it still creates new BDDs.This is not necessary as long as we are prepared to give up the simplifying of BDDs that is possible in set bounds propagation.
We do not represent domains of variables as BDDs, but rather as arrays of Boolean domains.A domain D is an array where, for variable v ranging over subsets of {1, . . ., N }: The BDD representation of a constraint B(c) is built as before.A significant difference is that since constraints only communicate through the set bounds of variables we do not need them to share a global variable order hence we can if necessary modify the variable order used to construct B(c) for each c, or use automatic variable reordering (which is available in most BDD packages) to construct B(c).Another advantage is that we can reuse the BDD for a constraint c(x) on variables x for the constraint c(ȳ) on variables ȳ (as long as they range over the same initial sets), that is, the same constraint on different variables.Hence we only have to build one such BDD, rather than one for each instance of the constraint.
The set bounds propagator sb(c(x)) for constraint c(x) is now implemented as follows.A generic BDD representation r of the constraint c(ȳ) is constructed.The propagator copies the domain description of the actual parameters x 1 , . . ., x n onto a domain description E for formal parameters y 1 , . . ., y n .It constructs an array E where E[ be the set of Boolean variables occurring in the constraint c(ȳ).The propagator executes the code bddprop(r, V, E) shown in Figures 4 and  5 which returns (r ′ , V ′ , E ′ ).If r ′ = F the propagator returns a false domain, otherwise the propagator copies back the domains of the formal parameters to the actual parameters so We will come back to the V ′ argument in the next subsection.The procedure bddprop(r, V, E) traverses the BDD r as follows.We visit each node n(v, f, t) in the BDD in a top-down memoing manner.We record if, under the current domain, the node can reach the F node, and if it can reach the T node.If the f child can reach the T node we add support for the variable v taking value 0. Similarly if the t child can reach T we add support for the variable v taking 1.If the node can reach both F and T we record that the variable v matters to the computation of the BDD.After the visit we reduce the variable set for the propagator to those that matter, and remove values with no support from the domain.The procedure assumes a global time variable which is incremented between each propagation, which is used to memo the marking phase.The top(n, V ) function returns the variable in the root node of n or the largest variable (under As presented bddprop has time complexity O(|r| × |V |) where |r| is the number of nodes appearing in BDD r.In practice the complexity is O(|r| + |V |) since the |V | factor arises from handling "long arcs", where a node n(v, f, t) has a child node (f or t) are labelled by a Boolean different from that next in the order ≺ after v.For set constraints the length of a long arc is typically bounded by the arity of the set constraint.It is possible to create a version of bddprop which is strictly O(|r|) by careful handling of long arcs.We did so, but in practice it was slower than the form presented here.bddprop has space complexity O(|V | + |r|) the first component for maintaining the domains of variables and the second for memoing the BDD nodes.
Also, no nodes for z 1 are actually visited, and the left node for y 2 only reaches F and the right node only reaches T .Hence matters[z 1 ] and matters[y 2 ] are not marked with the bddprop

Waking up Less Often
In practice a bounds propagation solver does not blindly apply each propagator until fixpoint, but keeps track of which propagators must still be at fixpoint, and only executes those that may not be.For set bounds this is usually managed as follows.To each set variable v is attached a list of propagators c that involve v. Whenever v changes, these propagators are rescheduled for execution.
We can do better than this with the BDD based propagators.The algorithm bddprop collects the set of Boolean variables that matter to the BDD, that is can change the result.If a variable that does not matter becomes fixed, then set bounds propagation cannot learn any new information.We modify the wakeup process as follows.Each propagator stores a list vars of Boolean variables which matter given the current domain.When a Boolean variable x j i becomes fixed we traverse the list of propagators involving x j i and wake those propagators where x j i occurs in vars.On executing a propagator we revise the set vars stored in the propagator.Note the same optimization could be applied to the standard approach, but requires the overhead of computing vars which here is folded into bddprop.It is possible to instead do propagator wake-up on literals, rather than variables.In this case, we observe that fixing a variable v to true matters to a node n(v, f, t) iff T is reachable from f and F is reachable from t -the converse holds for ¬v.In terms of the pseudo-code in Figure 5, the line While this allows for propagators to wake up less frequently, propagator execution is slower due to keeping track of additional reachability information.

Dead Subgraph Memoization and Shortcutting
The algorithm as presented above always explores all reachable parts of the graph in order to determine the set of supported values.However, a number of improvements for Multi-Decision Diagrams (MDDs) were presented by Cheng and Yap (2008) which reduce the portion of the graph which must be traversed in order to enforce consistency.These are dead subgraph memoization, which avoids traversal of subgraphs which cannot provide support for any values, and shortcutting, which recognizes situations where it is only necessary to find one path to T to ensure consistency.These can readily be adapted to a BDD-based set constraint solver.

Dead Subgraph Memoization
The key observation for dead subgraph memoization is that, as search progresses, paths along the graph to T are only ever removed.As such, if T becomes unreachable from a node n, the subgraph incident from n need never again be explored until the solver backtracks.Thus, if the set of dead nodes can be maintained, it is possible to progressively eliminate subgraphs during propagation.We keep for each instance of a constraint c(x) a failure set, fset which records which nodes can not reach T (and hence are equivalent to F).During propagation, once a node n is shown to have no path to T , it is added to the failure set fset.When a node is processed, we first check if it is in fset-if so, we terminate early, otherwise we proceed as normal.The modifications necessary for this are shown on the right in Figure 5.For simplicity the pseudo-code treats fset as a global.
A method for efficiently maintaining the failure sets was presented by Cheng and Yap (2008), which uses sparse-set data structures to provide efficient lookup, insertion and backtracking.The set fset is maintained as a pair of arrays: sparse and dense and a counter  These structures can be improved slightly by the observation that checking membership will occur significantly more often than insertion.Pseudo-code for the modified sparse-set operations are given in Figure 7(b).While insertion operations become more expensive, the overall computation time is reduced.
Example 5 Consider the set illustrated in Figure 8(a).The elements in the set are {1, 7}.We can determine that the element 4 is not in the set S 0 , as sparse[4] is not strictly less than members, indicated by the arrow in Figure 8(a).
To insert an element v using the standard sparse-set operations, we merely overwrite dense[members] with v, and set the value of sparse [v] to members.This is shown in Figure 8

Shortcutting
Shortcutting is an optimization to propagation on the BDD which notices that if all values in the current domains of variables v i , v i+1 , • • • , v N are fully supported, then we do not need to examine the rest of the nodes involving those variables.We keep a high water mark hwater which shows the least variable all of whose values are supported.If we ever reach a node numbered at or below the high water mark we only need to prove that it reaches T , we do not need to fully explore the sub-graph below it.
A modified propagation algorithm taking into account shortcutting (and dead subgraph minimization) is given in Figures 9 and 10.The high water mark hwater is originally larger than the greatest variable appearing in the BDD.
The principle difference of imp bddp is that if we reach a node with variable at or below the high water mark we use the simplified form shortcut bddp which only checks whether the node can reach T .The only other complexity is to update the high water mark hwater when we find all values of v are supported (E ). shortcut bddp has to be careful to mark all variables in nodes visited that reach T as mattering to the propagator.
Example 6 Consider the BDD for the constraint |y ∩ z| = 1 when N = 3 shown in Figure 11(a).As no variables are fixed, we first explore the false paths, and find the T node.This provides complete support for y 2 , x 3 , y 3 , so the high-water mark is updated to y 2 .When searching for support for x 2 false, we no longer need to find support for anything beneath the high-water mark -we need only find a single path to true from the node labelled y 2 .The high water mark then increases to y 1 .Likewise, when finding support for x 1 , everything below that point is already supported, so we explore only the first path to T .The edges explored are shown doubled in 11(b).
Example 6 also illustrates that the impact of shortcutting is highly dependent on the order in which branches are searched, and the structure of the constraint -if we were to shortcut bddp(node,V ,E) { if (in set(fset, node)) return (1,0); switch node { T : return (0,1); n(v, f, t): explore the true branches first, rather than the false branches, we would need to explore all nodes to find support for all variables.Clearly shortcutting does not change the asymptotic time or space complexity of the algorithm.Note that shortcutting for BDDs is more complex than the approach used by Cheng and Yap (2008) since they do not treat "long arcs" in MDDs.

Hybrid SAT Solver
Despite very fast propagation, a pure set bounds-based solver nevertheless suffers from an inability to analyze the reasons for failure, which results in repeated exploration of similar dead subtrees.This limits the performance of the solver on many hard problem instances.
In order to address this, we construct a hybrid solver which embeds BDD-based set bounds propagators within an efficient SAT solver.Search and conflict analysis are per- formed in the SAT solver, and the BDD propagators are used to generate inferences and clauses for the SAT solver to use during propagation.

Efficient Reason Generation
Key to a successful SAT solver is the recording of nogoods, small subsets of the current variable assignments which independently result in failure.This allows similar subtrees to be eliminated from consideration, hence significantly reducing the search space.
In order to construct nogoods, it is necessary to explain the reason why each literal was set. in order to determine the chain of reasoning which resulted in a contradiction.In a pure SAT solver this is easy, as each variable is either a decision variable, or associated with a clause that caused propagation.
BDD-based propagation methods, however, do not automatically provide explanations for inference.The naive approach for generating a reason clause for a BDD inference is to enumerate all the fixed variables which occur in the propagator, and construct a clause from the negations: Unfortunately, this often results in very large reason clauses, particularly in the case of merged propagators or global constraints.As smaller clauses result in stronger nogoods being generated by the SAT solver, it is preferable to determine the minimal set of variables required to cause propagation, and include only those variables in the clause.
A method for constructing such minimal clauses was demonstrated by Hawkins and Stuckey (2006), but this method involves constructing new BDDs, eliminating redundant variables until the minimal BDD is constructed, then reading off the variables remaining in the BDD.Given the propagation algorithm herein avoids expensive BDD operations, we do not wish to use them for explanation.
Given that a set of assignments {l 0 , . . ., l k } entail a literal l with respect to a constraint C, it is also true that As a result, the problem of finding a minimal reason for a given inference from a BDD is equivalent to fixing ¬l and unfixing as many variables as possible without rendering T reachable.
The algorithm presented by Subbarayan (2008) provides a method to do this by traversing a static graph, again avoiding the need to construct intermediate BDDs.The algorithm, given in Figure 12, traverses each node n(v, f, t) in a top-down memoing manner.At each node, it records if, given to the current domain, the T node is reachable.If the variable v has been assigned a value, it also records if T is reachable from the conflicting edge; any such edges must not become relaxed, otherwise the partial assignment is no longer a conflict.
The graph is then traversed a second time, this time in a breadth-first manner.For each variable v, if all nodes which have been reached corresponding to v variable may be relaxed without opening a path to T , the v is unfixed.If the v remains fixed, v is marked as part of the reason, and only the node corresponding to the value of v is marked as reachable.Otherwise, v is not in the minimal reason, and both the f and t nodes are marked as reached.The procedure returns the reason as a clause.The procedure is O(|r|) in time and space complexity, but note this is O(|r|) per new propagation that has to be explained!Example 7 Consider the constraint and assignments obtained in Example 4. It was deter- As such, the naive reason clause to explain 2 ∈ x would be ¬y 1 ∨ ¬z 2 ∨ x 2 ; however, it is possible to construct a smaller clause than this.
In order to construct the minimal reason for E[x 2 ] = {1}, we first set E[x 2 ] = {0}.The corresponding graph is shown in Figure 13(a), with nodes that are consistent with the partial assignment shown doubled.Note that as the solid edge from x 2 is not consistent with the assignment, T is not reachable along a doubled path from the root node.
The algorithm then determines the set of nodes from which T is reachable -these nodes shown doubled in Figure 13(b).These nodes must remain unreachable along the final reason; as such, the nodes which must remain fixed are the x 2 node and the leftmost z 2 node.
Finally, the algorithm progressively unfixes any variables which would not provide a path to T (in this case, y 1 ).The final path is shown in Figure 13(c), the resulting inference being

Lazy Reason Generation
The simplest way to use reason generation is a so called eager generation, where whenever a BDD propagator makes a new inference, a minimal reason clause is generated and added to Figure 12: Pseudo-code for the reason generation algorithm by Subbarayan (2008).Constructs a minimal set of variables required to cause the inference var = sign.the SAT solver.These clauses, however, cannot make any meaningful contribution to search until a conflict is detected -they cannot cause any propagation until the solver backtracks beyond the fixed variable, and no conflict clauses are constructed until there is a conflict.As there is a degree of overhead in adding and maintaining a large set of these clauses in the solver, it may be better to delay constructing these reasons until they are actually required to explain a conflict.
We can instead apply the reason generation only when the SAT conflict analysis asks for the explanation of a literal set by the BDD solver.We call this lazy generation.In order to do so, we must determine the state of the propagator which caused the inference.We implement this by recording the order in which literals become fixed in a propagator.When generating a reason for a variable v becoming fixed, we look at each variable in the propagator, and unfix any variable v ′ such that time(v) ≤ time(v ′ ), then restore them after the reason is constructed.

Hybrid Architecture
The hybrid SAT solver embeds BDD propagators inside the SAT engine.The architecture is illustrated in Figure 14.The usual SAT engine architecture is shown on the left.BDD propagation is added as shown on the right.Unit propagation causes Boolean literals to be fixed which may require that BDD propagators need to be awoken.We attach to each Boolean variable representing part of a set variable x the BDD propagators involving that set variable.When unit propagation reaches a fixpoint, the trail of fixed literals is traversed and each BDD propagator that includes one of these literals is scheduled for execution.If we are using filtering, it is only scheduled if the literal is one which matters to the propagator.Then we execute the scheduled BDD propagators using imp bddp.If the BDD propagator fixes some literals then these are added to the trail of the unit propagation engine.If we are using eager reason generation then we also immediately build a clause explaining the propagation and add it to the clause database and record this clause as the reason for the propagation of the literal.If we are using lazy reason generation, instead we record as the reason simply a pointer to the BDD propagator which causes the literal to be fixed.Then if conflict analysis demands an explanation for the literal, we call the reason generation for the BDD propagator, using the state at the time when the literal was fixed, to build an explaining clause.This is used in conflict analysis.We replace the reason for the literal in the trail by the generated explanation clause and also add the explanation clause to the database.
The implementation inherits almost all features of the underlying SAT solver.Eager reason clauses are added as nogoods, and deleted when the SAT solver decides to eliminate nogoods, lazy reason clauses are only generated on demand during conflict analysis.They are added to the clause database even though this is not necessary, since its makes memoing which explanations have been already performed simpler.The hybrid solver can make use of restarting activity based search, and restarts, although we also extend the search capabilities to allow some simple static searches as these can be preferable for the set problems we tackle.

Experimental Results
We built a hybrid SAT solver implementing the algorithms described above.The solver is based on MiniSAT 2.0 (dated 070721) (Eén & Sörensson, 2003), which has been modified to include the BDD-based propagation engine.BDDs are constructed using the BuDDy BDD package (http://sourceforge.net/projects/buddy/)All BDDs are constructed at the beginning of execution, then converted to the static graph used during propagation.Indeed, for many of the smaller problems solved in Section 6, the majority of the solution time is used in constructing the BDDs.
The BDD propagators are executed at a lower priority level than unit propagation, in order to detect conflict as early as possible.Reason clauses which are generated by the setbounds propagator are added to the SAT solver as learnt clauses, as otherwise the number of clauses added to the solver during propagation of hard problems can overwhelm the solver.
Experiments were conducted on a 3.00GHz Core2 Duo with 2 Gb of RAM running Ubuntu GNU/Linux 8.10.All problems were terminated if not completed within 10 minutes.
We experimented on 3 classes of set benchmarks: social golfers, Steiner systems, and Hamming codes.Unless otherwise specified, the hybrid solver is always executed using lazy reason generation.
We compare with the Gecode 3.1.0set bounds propagation solver since it is acknowledged as one of the fastest solvers available, as well as ECLiPSE 6.0 #100.We also compare with published results of the Cardinal (Azevedo, 2007) and Length-Lex (Yip & Van Hentenryck, 2009) solvers on the same problems.

Social Golfers
A common set benchmark is the "Social Golfers" problem, which consists of arranging N = g × s golfers into g groups of s players for each of w weeks, such that no two players play together more than once.Again, we use the same model as used by Lagoon and Stuckey (2004), using a w × g matrix of set variables v ij where 1 ≤ i ≤ w and 1 ≤ j ≤ g.
The global constraint partition < ensures its arguments are pairwise disjoint and imposes a lexicographic order on its arguments, i.e.
The corresponding propagator is based on a single BDD.We construct BDD propagators for each of the constraint forms The hybrid solver constructs one BDD for each of the 4 terms in the above equation, instantiating constraints accordingly.
Table 1 shows the results using a static search strategy on easy problems.The search fixes the elements of the sets v ij is order v 11 , v 12 , . . ., v 1g , v 21 , . . ., v wg , always trying to first place the least element in the set then excluding it from the set.We compare against the reported results for the original BDD-SAT hybrid solver of Hawkins and Stuckey (2006) versus a number of variations of our hybrid.base is the base solver of Figures 4 and 5, while +f indicates with filtering of Section 4.1 added, +s indicates with dead subgraph memoization and shortcutting added (Section 4.2) using the original sparse set code, +i is these optimizations with the improved sparse set code.We also combine filtering with the other optimizations.The table shows time and number of fails for each variant, where the solvers with identical failure behaviour are grouped together.Note that filtering can change the search by reordering the propagations and hence changing the nogoods that are generated, while the other optimizations cannot except that shortcutting can change the results of filtering (and hence change search).While filtering improves on the base line, dead subgraph memoization and shortcutting do not, although we can see the benefit of the improved sparse set operations.Comparing against the solver of Hawkins and Stuckey (2006), which was run on a (♯) 2.4GHz Pentium 4, we find that, slightly different number of backtracks and slightly faster machine not withstanding, the solver presented here is roughly an order of magnitude faster.
Table 2 shows the results using VSIDS search on easy problems.It compares against the solver of Hawkins and Stuckey (2006) 1, and overall VSIDS is better than static search.The table illustrates some of the difficulty of comparing systems using VSIDS search, since small differences can drastically change the search space.The solver +f is the best except for a bad-performance on 7,5,3.The base solver is around 5 times faster per failure than the solver of Hawkins and Stuckey (2006).The Tseitin decomposition is not competitive, even if we discount the results on 7,5,3.For social golfers, dead-subset memoization and shortcutting provide no advantage (when we discount the drastically different search for 7,5,3 using VSIDS).While the number of nodes processed can be reduced slightly, this is not enough to repay the additional cost of computation at each node.
Table 3 compares the reason generation strategies: eager reasoning which constructs reasons as soon as inference is detected; and lazy reasoning which only those reasons necessary to determine the first UIP or perform conflict clause minimization.
Table 3 compares the base solver with and without filtering (since dead subgraph memoization and shortcutting do not help here) on harder social golfer problems using a static search.It shows time (base) as well as the number of reasons generated and fails in order to find a first solution.For these harder examples filtering is highly beneficial.Here we can see that the number of reasons generated by lazy reasoning is about half of that required by eager reasoning, but it doesn't make that much difference to the computation time, since propagation dominates the time spent in the solver.Interestingly not adding reasons eagerly also seems to generate slightly better nogoods as the search is usually smaller.Table 3: First-solution performance results on harder Social Golfers problems, using a static least-element in set search method.Results are given comparing eager and lazy reason generation.
of lazy reasoning are increased by the use of VSIDS, presumably because the better nogoods are then more useful in driving search.Finally Table 5 compares against a number of different systems.We use the model of social-golfers described in the work of Yip and Van Hentenryck (2009), which in addition fixes the first week, and the first group of the second week to eliminate symmetric solutions.We use the instances reported by Yip and Van Hentenryck (2009).We show results for our base solver with and without filtering.We compare against Gecode 3.1.0and Eclipse 6.0 #100, both which implement a set bounds propagation combined with limited cardinality reasoning, on an identical MiniZinc model of social-golfers running on our 3GHz Core2Duo.Gecode arguably represents the state of the art for set bounds propagation solving.We also compare against the published results of the Cardinal solver (Azevedo, 2007) & Van Hentenryck, 2006;Yip & Van Hentenryck, 2009 for details) running on a ( ‡) C2D-M 2.53GHz machine.The pure set bounds solvers cannot compete with our approach since the search space without using nogood recording is just too big.None of the other systems except Length-lex can solve all of these instances.One can a see a drastic difference between number of failures for Gecode, which uses set bounds propagation without learning, versus our base solver.Gecode can sometimes require less failures on easy problems since it combines cardinality reasoning with bounds reasoning, but on hard problems the advantages of learning to prune similar searches in other parts of the tree dominates completely.The stronger pruning of Length-lex compared to set bounds means it can often improve on fails compared to base but learning is more robust.The hybrid solver is overall around an order of magnitude faster than Length-Lex.

Steiner Systems
Another commonly used benchmark for set constraint solvers is the calculation of small Steiner systems.A Steiner system S(t, k, N ) is a set X of cardinality N and a collection C of subsets of X of cardinality k (called 'blocks'), such that any t elements of X are in exactly one block.Any Steiner system must have exactly m = N t / k t blocks (Theorem 19.2 of van Lint & Wilson, 2001).
We model the Steiner problem similarly to Lagoon and Stuckey (2004) extended for the case of more general Steiner Systems.We model each block as a set variable s 1 , . . ., s m , with the constraints: For comparison with the results of Azevedo (2007) and Yip and Van Hentenryck (2009) A formulation for this problem is: is the symmetric difference.This is similar in structure to the formulation for the Steiner Systems; however, rather than having a fixed number of sets, we find the maximal code by repeatedly adding new sets and the corresponding constraints until no solution can be found.The unsatisfiability of n codewords proves that the maximal code has n − 1 codewords.We create BDD propagators for the constraint form We compare on two different models of the fixed-weight Hamming code problems, one just using the description above, and another where the first two sets are fixed to remove symmetries.We compare against Gecode, the published results of Length-Lex with our hybrid using a static search strategy (the same least element in set strategy as used for Social Golfers), as well as the hybrid solver and a Tseitin decomposition using VSIDS search.For our systems we compare with and without shortcutting and our optimized implementation.Since we are not sure which model was used by Length-Lex we report it results for both models.
Tables 7 to 10 show the results on the 11 hard instances reported by Hawkins et al. (2005).Clearly on these problems the VSIDS hybrid is the most robust.It can solve all but one instance in the basic model, and all with the additional symmetry breaking.This example also clearly shows the potential advantages of shortcutting and our improved data structures: these do not change the search but improve the time by 18% and 21% respectively for the base model, and 24% and 26% respectively for the improved model.Once more Tseitin decomposition is not competitive.

Related Work
Set-constraint problems have been an active area of research in the past decade.Many of the earlier solvers, beginning with PECOS (Puget, 1992), used the set-bounds repre-sentation combined with a fixed set of propagation rules for each constraint.This general approach was also used by Conjunto (Gervet, 1997), ECL i PS e (IC-PARC, 2003), ILOG Solver (ILOG, 2004) and Mozart (Müller, 2001).However, as set-bounds are a relatively weak approximation of the domain of a set variable, a variety of variations have been developed to improve the propagation strength of set-constraint solvers.These include solvers which combine set-bounds representation with either cardinality information, such as that proposed by Azevedo (2002Azevedo ( , 2007)), lexicographic bounds information (Sadler & Gervet, 2004) or both (Gervet & Van Hentenryck, 2006;Yip & Van Hentenryck, 2009).
BDD-based approaches to set-constraint solving, such as that presented by Hawkins et al. (2005) differs greatly from these approaches, as it is possible to perform propagation over arbitrary constraints; Lagoon and Stuckey (2004) also demonstrated the feasibility of a BDD-based solver which maintains a complete domain representation of set variables.These directly BDD-based algorithms were used to construct the earlier hybrid solver presented by Hawkins and Stuckey (2006), which is conceptually similar to the solver presented in this paper.The solver presented here is much more efficient, and includes improvements such as filtering and shortcutting not present in the solver of Hawkins and Stuckey (2006).The solver of Damiano and Kukula (2003) also combines BDD solving and SAT solving, but rather than building BDDs from a high-level problem description and lazily constructing a SAT representation, instead takes a CNF SAT representation and constructs a BDD from a collection of clauses with the primary goal of variable elimination.It is essentially equivalent to the base solver.
The underlying BDD propagation algorithm is similar to propagation of the case constraint of SICStus PRolog (SICS, 2009) and Multi-valued Decision Diagrams (MDDs) (see e.g., Cheng & Yap, 2008).Indeed we have adapted the dead subgraph memoization and shortcutting devices of Cheng and Yap (2008) to BDD propagation.Propagators for case and MDDs do not presently use filtering or generate reasons.
Finally the hybrid set solver we present in this paper is an example of a lazy clause generation solver (Ohrimenko, Stuckey, & Codish, 2007, 2009).The BDD propagators can be understood as lazily creating a clausal representation of the set constraints encoded in the BDD, as search progresses.

Concluding Remarks
In this paper we have improved BDD-based techniques for set-bounds propagation, having demonstrated an approach which avoids the need for expensive BDD construction and manipulation operations.This traversal-based method, when combined with filtering to reduce the number of redundant propagator executions and dead subgraph memoization and shortcutting, is at least an order of magnitude faster than previous techniques which construct BDDs during runtime (Hawkins et al., 2005).
Furthermore, when integrated into a modern SAT solver with clause learning and augmented with a method for generating nogoods, the new hybrid solver is capable of solving hard problem instances several orders of magnitude faster than pure bounds set solvers.Overall the hybrid solver is robust and highly competitive with any other propagation based set-solvers we are aware of.
In many set problems there are significant numbers of symmetries and there is a large body of work solving set problems with symmetry breaking techniques (see e.g., Puget, 2005).It would be interesting to combine symmetry breaking with our hybrid solver.

Acknowledgments
Part of this work was published previously (Gange, Lagoon, & Stuckey, 2008).NICTA is funded by the Australian Government as represented by the Department of Broadband, Communications and the Digital Economy and the Australian Research Council.

Figure 1 :
Figure 1: Architecture for the SAT solver.

Figure 3 :
Figure 3: Pseudo-code for Tseitin transformation of BDD rooted at node where n ′ is the Boolean variable encoding the truth value of node.

Example 4
Consider the BDD for the constraint x = y ∪ z when N = 2 shown in Figure 6(a).Assuming a domain E where E[y 1 ] = {1} (1 ∈ y) and E[z 2 ] = {1} (2 ∈ z), and the remaining variables take value {0, 1}, the algorithm traverses the edges shown with double lines in Figure 6(b).No path from x 1 , or x 2 following the f arc reaches

Figure 5 :
Figure 5: Pseudo-code for processing the constraint graph during propagation.Modifications necessary for using dead-subgraph memoization are shown on the right.

Figure 7 :
Figure 7: Pseudo-code for conventional sparse-set operations, and the corresponding modified versions.
(b), inserting 3 into S 0 .At this point, both sparse[3] and sparse[4] have the value 2. To test if 4 ∈ S ′ 0 , it is not sufficient to determine if sparse[4] < members.One must also check that dense[sparse [4]] = 4.When inserting v using the modified operations, as illustrated in Figure 8(c), we swap the values of sparse[v] and sparse [dense[members]], and likewise switch the values of dense[members] and dense[sparse [v]].This maintains the property that v ∈ S ⇔ sparse[v] < members.
Figure 8: A sparse representation for sets.(a) A possible state of the data structure representing S 0 = {1, 7}.(b) Inserting 3 into the data structure using the standard operations.sparse [3] is updated to point to the next element of dense, and the corresponding entry in dense points back to 3. Notably, both sparse [3] and sparse [4] now point to dense[2].(c) Inserting 3 into the data structure using the modified operations.After the operation, both the sparse and dense arrays are maintained such that ∀v dense[sparse[v]] = v.

Figure 10 :
Figure 10: Pseudo-code for the shortcut phase.

Figure 11 :
Figure 11: (a) The BDD representing |x ∩ y| ≤ 1 where N = 3.A node n(v, f, t) is shown as a circle around v with a dotted arrow to f and full arrow to t.(b) The edges traversed by imp bddp, when E[v] = {0, 1} for all v, are shown doubled.

Figure 13 :
Figure 13: (a) The BDD representing x = y ∪ z where N = 2, with E[y 1 ] = {1}, E[z 2 ] = {1} and E[x 2 ] = {0}.Edges consistent with the partial assignment are shown doubled.(b) Nodes which must remain unreachable in the reason are shown doubled.(c) Edges reachable along the minimal reason are shown doubled, as are nodes which remain fixed.

Figure 14 :
Figure 14: Architecture for the hybrid BDD-SAT solver.

Table 1 :
and a Tseitin decomposition.The results are the Fast Set Bounds Propagation Using a BDD-Sat Hybrid First-solution performance results on the Social Golfers problem using a static, first-element in set ordering.Instances marked with (⋆) are unsatisfiable, entries marked with '-' did not complete within 10 minutes.
same as for Table

Table 2 :
Table4shows the results using VSIDS search on these harder instances.It appears the advantages First-solution performance results on the Social Golfers problem using a VSIDS search strategy.

Table 4 :
Yip and Van Hentenryck (2009)erformance results on harder Social Golfers problems, using a VSIDS search method.Results are given comparing eager and lazy reason generation.morecomplexcardinality reasoning for set solving, using a ( †) Pentium 4 2.4GHz machine, and the recently published results for the Length-Lex solver ofYip and Van Hentenryck (2009), which maintains bounds on sets variables in terms of the length-lex order (see Gervet

Table 7 :
, we construct a dual model with additional variables d 1 , . . ., d N , with additional constraints Results on hard Hamming instances with a static least-element in set search order, with no additional symmetry breaking.

Table 8 :
Results on hard Hamming instances with a VSIDS, with no additional symmetry breaking.

Table 9 :
Results on hard Hamming instances with a sequential least-element in set search order, with fixed first and second sets.

Table 10 :
Results on hard Hamming instances with a VSIDS search strategy, with fixed first and second sets.