% vim:tw=70:
\documentclass[twoside,11pt]{article}
\usepackage[active]{srcltx}

\usepackage{phsty}
\usepackage{theapa}
\usepackage{jair}
\usepackage{amssymb, amsmath, amsthm, stmaryrd, url}
\usepackage{hhline}
\usepackage{multirow}
\usepackage[all]{xy}
%\CompileMatrices
\ifx\pdftexversion\undefined
  \usepackage[dvips]{graphicx}
\else
  \usepackage[pdftex]{graphicx}
\fi



\theoremstyle{plain}
\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}

\theoremstyle{definition}
\newtheorem{example}{Example}[section]

\DeclareMathOperator{\chr}{\mathit{char}}
\DeclareMathOperator{\dom}{\mathit{dom}}
\DeclareMathOperator{\setb}{\mathit{sb}}
\DeclareMathOperator{\sbc}{\mathit{sbc}}
\DeclareMathOperator{\sblex}{\mathit{sblex}}
\DeclareMathOperator{\vars}{\mathit{vars}}
\DeclareMathOperator{\solns}{\mathit{solns}}
\DeclareMathOperator{\solv}{\mathit{solv}}
\DeclareMathOperator{\conv}{\mathit{conv}}
\DeclareMathOperator{\ran}{\mathit{ran}}
\DeclareMathOperator{\ite}{\mathit{ite}}
\DeclareMathOperator{\VAR}{\mathit{VAR}}
\DeclareMathOperator{\plus}{\mathit{plus}}
\DeclareMathOperator{\shl}{\mathit{shl}}
\DeclareMathOperator{\mul}{\mathit{mul}}
\DeclareMathOperator{\bmul}{\mathit{bmul}}
\DeclareMathOperator{\lexlt}{\mathit{lexlt}}
\DeclareMathOperator{\reverse}{\mathit{reverse}}
\DeclareMathOperator{\wsum}{\mathtt{wsum}}
\DeclareMathOperator{\domdc}{\mathit{domdc}}
\DeclareMathOperator{\partition}{\mathtt{partition}}
\DeclareMathOperator{\atmost}{\mathtt{atmost}}
\DeclareMathOperator{\occ}{\mathit{occ}}
\def\True{1}
\def\False{0}



\newcommand{\bexists}{\overline{\exists}}
\newcommand{\eclipse}{ECL$^i$PS$^e$}
\newcommand{\entails}{\vDash}
\newcommand{\essential}[1]{\left \llbracket #1 \right \rrbracket}
\newcommand{\fail}{\bf fail}
\newcommand{\lra}{\leftrightarrow}
\newcommand{\lu}{{L\overline{U}}}
\newcommand{\la}{\rightarrow}
\newcommand{\ra}{\rightarrow}
\newcommand{\card}[2]{|{#1}| = {#2}}
\newcommand{\cardl}[2]{|{#1}| \geq {#2}}
\newcommand{\cardu}[2]{|{#1}| \leq {#2}}
\newcommand{\cardlu}{\mathit{card}}
\newcommand{\range}[2]{\left[#1,#2\right]}
\newcommand{\s}[1]{\left\{#1\right\}} % sets
\newcommand{\n}[1]{\left\langle #1 \right\rangle} % tuples
\newcommand{\llbrace}{\{\!\!\{}
\newcommand{\rrbrace}{\}\!\!\}}

\newcommand{\xyo}[1]{*+[o][F-]{#1}}
\newcommand{\xyoo}[2]{*+[o][F-]{#1}{\n{#2}}}
\newcommand{\xytrue}{*+[F]{\mathrm{\True}}}
\newcommand{\xyfalse}{*+[F]{\mathrm{\False}}}



\newcommand{\fixme}[1]{\textbf{XXX: #1 } }
\newcommand{\ph}[1]{\marginpar{\sc ph}\textbf{FIXME: #1}}
\newcommand{\pjs}[1]{\marginpar{\sc pjs}\textbf{FIXME: #1}}
\newcommand{\vit}[1]{\marginpar{\sc vit}\textbf{FIXME: #1}}

\newcommand{\ignore}[1]{}


%\jairheading{?}{2004}{?}{?}{?}
\ShortHeadings{Solving Set Constraint Satisfaction Problems using ROBDDs}{Hawkins, Lagoon, \& Stuckey}
\firstpageno{109}

\begin{document}

\title{Solving Set Constraint Satisfaction Problems using ROBDDs}

\author{\name Peter Hawkins \email hawkinsp@cs.mu.oz.au \\
  \name Vitaly Lagoon \email lagoon@cs.mu.oz.au \\
  \addr Department of Computer Science and Software Engineering \\
  The University of Melbourne, VIC 3010, Australia
  \AND
  \name Peter J. Stuckey \email pjs@cs.mu.oz.au \\
  \addr NICTA Victoria Laboratory, 
  Department of Computer Science and Software Engineering \\
  The University of Melbourne, VIC 3010, Australia}

\jairheading{24}{2005}{109--156}{10/04}{7/05}

\maketitle

\begin{abstract}
  In this paper we present a new approach to modeling finite set domain
  constraint problems using Reduced Ordered Binary Decision Diagrams
  (ROBDDs). We show that it is possible to construct an efficient
  \emph{set domain} propagator which compactly represents many set
  domains and set constraints using ROBDDs.  We demonstrate that the
  ROBDD-based approach provides unprecedented flexibility in modeling
  constraint satisfaction problems, leading to performance improvements.
  We also show that the ROBDD-based modeling approach can be extended to
  the modeling of integer and multiset constraint problems in a
  straightforward manner. Since domain propagation is not always
  practical, we also show how to incorporate less strict consistency
  notions into the ROBDD framework, such as set bounds, cardinality
  bounds and lexicographic bounds consistency. Finally, we present
  experimental results that demonstrate the ROBDD-based solver performs
  better than various more conventional constraint solvers on several
  standard set constraint problems.

\end{abstract}

\section{Introduction}
\label{sect:intro}

It is often natural to express a constraint satisfaction problem (CSP)
using finite domain variables and relations between those variables,
where the values for each variable are drawn from a finite universe of
possible values.  One of the most common methods of solving finite
domain CSPs is through maintaining and updating a domain for each
variable using the combination of a backtracking search and an
incomplete local propagation algorithm. The local propagation algorithm
attempts to enforce consistency on the values in the variable domains by
removing values that cannot form part of a solution to the system of
constraints. Various levels of consistency can be defined, each with
associated costs and benefits. Most consistency algorithms are
incomplete---that is, they are incapable of solving the problem by
themselves, so they must be combined with a backtracking search
procedure to produce a complete constraint solver.

When attempting to apply this general scheme to the task of solving
constraint satisfaction problems over finite set variables we quickly
run into practical problems. Since the universe of possible values for a
set variable is usually very large, the na\"ive representation of the
domain of a set variable as a set of sets is too unwieldy to solve
realistic problems. For example, if a set variable can take on the value
of any subset of the set $\{1,\ldots,N\}$, then its domain contains
$2^N$ elements, which quickly becomes infeasible to represent as $N$
increases in magnitude. Accordingly, most set constraint solvers to date
have used an approximation to the true domain of a set variable in order
to avoid this combinatorial explosion.

The most common approximation to the true domain of set variable $v$ has
been to maintain an upper bound $U$ and a lower bound $L$ of the domain
under the subset partial ordering relation and to perform set
\emph{bounds} propagation on these bounds. That is, $L$ contains
elements that must be in the set $v$, and $U$ is the complement of the
set of elements that must not be in $v$.  Conventionally, a fixed set of
inference rules specific to each constraint are used to enforce
consistency on these upper and lower bounds.  This basic scheme was
proposed by \citeA{pecos}, and has been implemented by set solvers such
as \textsc{Conjunto} \cite{gervet:conjunto}, the \texttt{fd\_sets} and
\texttt{ic\_sets} libraries of \eclipse{} \cite{eclipse},
\textsc{Mozart} \cite{mueller:mozart}, and \textsc{ILOG Solver}
\cite{ilog}.

Set bounds are a crude approximation to a set domain at best, and thus
various refinements to the bounds propagation scheme have been
proposed so as to more effectively capture the nature of a set domain.
%These include maintaining upper and lower bounds on the cardinality of a set
%domain \cite{azevedo:cardinal,mueller:mozart}, or more recently
%maintaining upper and lower bounds under a total lexicographic
%ordering
%\cite{sadlerg:lexbounds}. While providing more effective propagation
\citeA{azevedo:cardinal} demonstrated that maintaining and performing
inferences upon the upper and lower bounds on the cardinality of a set
domain leads to a significant performance improvement on a variety of
problems. While earlier set solvers such as \textsc{Conjunto} also
maintained cardinality bounds, only partial usage was made of this
information.  More recently, \citeA{sadlerg:lexbounds} showed that
incorporating upper and lower bounds under a lexicographic ordering
leads to significantly stronger propagation, albeit at the cost of a
marked increase in propagation time, leading to only a marginal
performance improvement overall.  While both of these approaches provide
more effective propagation than the simple set bounds scheme, they do
not approach the effectiveness of a true set domain propagator, which
ensures that every value in the domain of a set variable can be extended
to a complete assignment of every variable in any given constraint.

Consequently, we would like to devise a representation for set domains
and constraints on those domains that is tractable enough to permit
domain propagation.  The principal observation that permits the
implementation of a set domain propagator is that a finite integer set
$v$ can be represented by its characteristic function $\chi_v$:
\[
\chi_v:\mathbb{Z} \ra \{0,1\} \text{ where } \chi_v(i) = 1 \text{ iff } i
\in v
\]
Accordingly we can use a set of Boolean variables $v_i$ to represent the
set $v$, which correspond to the propositions $v_i \lra i \in v$.  We
can describe set domains and set constraints in terms of these Boolean
variables. Interestingly, set bounds propagation as described above is
equivalent to performing domain propagation in a na\"ive way on this
Boolean representation.

In this paper we investigate this Boolean modeling approach for modeling
finite domain constraints, and in particular set constraints.  We show
that it is possible to represent the domains of set variables using
Reduced Ordered Binary Decision Diagrams (ROBDDs), a data structure for
representing and manipulating Boolean formul\ae{}.  Such representations
are usually fairly compact, even if they correspond to very large
domains. In addition, it is possible to represent the constraints
themselves as ROBDDs, permitting us to produce efficient \emph{set
  domain} constraint propagators solely using ROBDD operations.

The ROBDD-based representation allows us to easily conjoin constraints
and existentially quantify variables, thus permitting us to remove
intermediate variables and merge constraints for stronger propagation.
The construction of global constraints becomes an almost trivial
exercise, without the requirement to write laboriously new propagators
for every new constraint that we would like to use.

We also demonstrate that only minor changes are needed to the operation
of the set domain propagator in order to implement other, less strict
notions of consistency. In particular, we show how to construct set
bounds, cardinality bounds and lexicographic bounds propagators, and how
to utilize these notions to produce a \emph{split domain} solver that
combines bounds and domain reasoning.

A key theme of this paper is the flexibility of the ROBDD-based modeling
approach. We are not limited only to modeling set variables---for
example ROBDDs can be used to model integer variables and integer
constraints. While the ROBDD-based approach to modeling finite domain
integer variables is in general not as efficient as many existing finite
domain integer solvers, the ability to model integer constraints using
ROBDDs is an essential building block which allows us to construct set
constraints such as cardinality and weighted sum, as well as allowing us
to represent multisets and multiset constraints.

Finally, we present experiments using a variety of standard set
problems to demonstrate the advantages of our modeling approach.

Many of the ideas from this paper have been previously published in two
previous works \cite{cp04-sets,austai04-sets}. This paper contains a
more complete exposition of those ideas, as well as important extensions
of the work that has previously been presented. These include
substantial improvements in the complexity of the propagation algorithm,
as well as cardinality and lexicographic bounds solvers implemented
using ROBDDs.  In addition, we show how to model integer expressions,
allowing us to implement a weighted-sum constraint, and to model
multisets and multiset constraints.  Using this, we present new
experimental results for the Hamming Code and Balanced Academic
Curriculum problems.  Finally, we present results comparing the
ROBDD-based solver against a solver with good cardinality reasoning
(\textsc{Mozart}).
% and Progressive Party problems.

The remainder of this paper is structured as follows.
Section~\ref{sect:prelim} contains essential concepts and definitions
necessary when discussing finite domain solvers and ROBDDs.
Section~\ref{sect:modeling} shows how to model set domains and set
constraints using ROBDDs, and presents a basic outline of an ROBDD-based
set constraint solver. Section~\ref{sect:improvements} demonstrates how
to improve the performance of the ROBDD-based set solver through the
construction of global constraints, the removal of intermediate
variables, and through symmetry-breaking approaches.
Section~\ref{sect:ints} demonstrates how to model integer and multiset
expressions,
as well as how to implement a weighted-sum constraint for set and multiset
variables.
Section~\ref{sect:propagation} shows
how to construct a more efficient domain propagator, as well as set
bounds, set cardinality, and lexicographic bounds propagators. Finally,
in Section~\ref{sect:exp} we present experimental results comparing the
ROBDD-based solver with more conventional set solvers on a variety of
standard problems.

\section{Preliminaries}
\label{sect:prelim}

In this section we define the concepts and notation necessary when
discussing propagation-based constraint solvers. These definitions are
largely identical to those presented by \citeA{cp04-sets} and others.
For simplicity we shall present all of our definitions in the case of
finite set variables; the extensions to multiset and integer variables
are trivial.

\subsection{Lattices, Domains, and Valuations}

Let $\mathcal{L}$ be the powerset lattice $\langle
\mathcal{P}(\mathcal{U}), \subset \rangle$, where $\mathcal{P}(x)$
denotes the powerset of $x$ and the \emph{universe} $\mathcal{U}$ is a
finite subset of $\mathbb{Z}$.  We say a subset $M \subseteq
\mathcal{L}$ is \emph{convex} if for any $a, c \in M$ the relation $a
\subseteq b \subseteq c$ implies $b \in M$ for all $b \in \mathcal{L}$.
The \emph{interval} $[a,b]$ is the set $M = \{ x \in \mathcal{L} ~|~ a
\subseteq x \subseteq b\}$.  Intervals are obviously convex.
%\pjs{interval and convex are the same things (at least for us)?}
Given a subset $K \subseteq \mathcal{L}$, we define the convex closure
of $K$:
\[
\conv(K) = \left [\bigcap_{x \in K} x, \bigcup_{x \in K} x \right ]
\]
The convex closure operation satisfies the properties of extension ($x
\subseteq \conv(x)$), idempotence ($\conv(x) = \conv(\conv(x))$), and
monotonicity (if $x \subseteq y$ then $\conv(x) \subseteq(y)$)
\cite{gervet:conjunto}.

\begin{example}
  The set $X = \{ \{1\}, \{1,3\}, \{1,4\}, \{1,3,4\} \}$ is convex and
  equivalent to the interval $[\{1\}, \{1,3,4\}]$. Conversely, the set
  $Y = \{ \{1\}, \{1,3\}, \{1,3,4\} \}$ is not convex. However, the
  convex closure of $Y$ is precisely the interval $X$, i.e. $\conv(Y) =
  X$.
\end{example}

Let $\mathcal{V}$ denote the fixed finite collection of all set
variables.  Each variable has a domain, which is a finite collection of
possible values from $\mathcal{L}$ (which are themselves sets).  More
generally, we define a \emph{domain} $D$ to be a complete mapping from
the set $\mathcal{V}$ to finite collections of finite sets of integers.
When we speak of the \emph{domain of a variable} $v$, we mean $D(v)$. A
domain $D_1$ is said to be \emph{stronger} than a domain $D_2$, written
$D_1 \sqsubseteq D_2$, if $D_1(v) \subseteq D_2(v)$ for all $v \in
\mathcal{V}$. Two domains are said to be equal, written $D_1 = D_2$, if
$D_1(v) = D_2(v)$ for all $v \in \mathcal{V}$.  We call a domain $D$ a
\emph{range domain} if $D(v)$ is an interval for all $v \in
\mathcal{V}$.  We extend the concept of convex closure to domains by
defining $\ran(D)$ to be the unique (range) domain such that $\ran(D)(v)
= \conv(D(v))$ for all $v \in \mathcal{V}$.

A \emph{valuation} is a function from $\mathcal{V}$ to $\mathcal{L}$,
which we write using the mapping notation $\{ v_1 \mapsto d_1, v_2
\mapsto d_2, \ldots, v_n \mapsto d_n \}$, where $v_i \in \mathcal{V}$
and $d_i \in \mathcal{L}$ for all $i \in \mathbb{N}$.  Clearly a
valuation can be extended to constraints involving the variables in the
obvious way. We define $\vars$ to be the function that returns the set
of variables that are involved in a constraint, expression or valuation.
In an abuse of notation, a valuation $\theta$ is said to be an element
of a domain $D$, written $\theta \in D$, if $\theta(v) \in D(v)$ for all
$v \in \vars(\theta)$.

We say a domain $D$ is a \emph{singleton} or \emph{valuation} domain if
$|D(v)| = 1$ for all $v \in \mathcal{V}$. In this case $D$ corresponds
to a single valuation $\theta$ where $D(v) = \{\theta(v)\}$ for all $v
\in \mathcal{V}$.

\subsection{Constraints, Propagators, and Propagation Solvers}

A \emph{constraint} is a restriction placed upon the allowable values
for a collection of variables. We are interested in the following
\emph{primitive set constraints},
where $k$ is an integer, $d$
is a ground set value, and $u$, $v$ and $w$ are set variables: 
$k \in v$ (membership), $k \notin v$
(non-membership), $u = v$ (equality), $u = d$ (constant equality),
%%PJS not used%%$u \subset v$ (strict subset),
$u \subseteq v$ (non-strict subset),
$u = v \cup w$ (union),
$u = v \cap w$ (intersection),
$u = v \setminus w$ (set difference),
$u = \overline{v}$ (complement),
$u \neq v$ (disequality),
$|u| = k$ (cardinality),
$|u| \geq k$ (lower cardinality bound), and
$|u| \leq k$ (upper cardinality bound). Later
we shall introduce non-primitive set constraints which are formed by
composing primitive set constraints.

We define the \emph{solutions} of a constraint $c$ to be the set of
valuations that make that constraint true, i.e. $\solns(c) = \{ \theta
~|~ \vars(\theta) = \vars(c) \text{ and } \entails \theta(c) \}$.

\begin{example}\label{ex:subset}
  Suppose $v$ and $w$ are set variables, and $D$ is a domain with $D(v)
  = \{ \{1\}, \{1,3\}, \{2,3\} \}$, and $D(w) = \{ \{2\}, \{1, 2\}, \{1,
  3\} \}$.  Let $c$ be the constraint $v \subseteq w$.  Then the
  solutions to $c$ in the domain $D$ are the valuations $\{ v \mapsto
  \{1\}, w \mapsto \{1,2\}\}$, $\{ v \mapsto \{1\}, w \mapsto
  \{1,3\}\}$, and $\{ v \mapsto \{1,3\}, w \mapsto \{1,3\}\}$.
\end{example}

With every constraint we associate a \emph{propagator} $f$, which is a
monotonic decreasing function from domains to domains, so if $D_1
\sqsubseteq D_2$ then $f(D_1) \sqsubseteq f(D_2)$
and $f(D_1) \sqsubseteq D_1$. 
A propagator $f$ is said
to be \emph{correct} for a constraint $c$ if:
\[
    \{ \theta ~|~ \theta \in D \} \cap \solns(c) = 
    \{ \theta ~|~ \theta \in f(D) \} \cap \solns(c)
\]
Correctness is not a strong restriction, since the identity propagator
is correct for all constraints $c$. We usually assume that all
propagators are both correct and \emph{checking}, that is, if $D$ is a
singleton domain formed through propagation then its associated
valuation $\theta$ makes $\theta(c)$ true.

We can use propagators to form the basis for a constraint solver. A
\emph{propagation solver} $\solv(F,D)$ takes a set of propagators $F$
and a current domain $D$, and repeatedly applies propagators from $F$ to
the current domain until a fixpoint is reached. In other words $\solv(F,
D)$ is the weakest domain $D' \sqsubseteq D$ which is a fixpoint (i.e.
$f(D') = D'$) for all $f \in F$. This fixpoint is unique
\cite{apt:propagation}.

\subsection{Local Consistency}

The notion of local consistency is of importance when considering the
solution of constraint satisfaction problems. We can define various
levels of consistency of different strengths and levels of difficulty to
enforce. We describe two here.


A domain $D$ is said to be \emph{domain consistent} for a constraint $c$
if $D$ is the strongest domain which contains all solutions $\theta \in
D$ of $c$; in other words there does not exist a domain $D' \sqsubseteq
D$ such that $\theta \in D$ and $\theta \in \solns(c)$ implies $\theta
\in D'$. A set of propagators maintains domain consistency if
$\solv(F,D)$ is domain consistent for all constraints $c$.

\begin{definition}
\label{def:dom-propagator}
A \emph{domain propagator} for a constraint $c$ is a function $\dom(c)$
mapping domains to domains which satisfies the following identity:
\begin{equation*}
  \dom(c)(D)(v) =
  \begin{cases}
    \{ \theta(v)\ |\ \theta \in D \ \wedge\ \theta \in
    \solns(c) \} & \text{if } v \in \vars(c) \\
    
    D(v) & \text{otherwise}
  \end{cases}
\end{equation*}
\end{definition}

\begin{lemma}
  A domain propagator $\dom(c)$ for a constraint $c$ is correct,
  checking, monotonic, and idempotent.
\end{lemma}
\begin{proof}
    Straightforward from the definitions.
\end{proof}

\begin{example}\label{ex:subset2}
  Consider the constraint $v \subseteq w$ and domain $D$ of
  Example~\ref{ex:subset}.  The domain propagation $D' = dom(c)(D)$
  returns domain $D'$ where $D'(v) = \{ \{1\}, \{1,3\}\}$ and $D'(w) =
  \{\{1,2\}, \{1,3\}\}$.  The missing values are those that did not take
  part in any solution.
\end{example}

Domain consistency may be difficult to achieve for set constraints, so
instead we often need a weaker notion of consistency. The notion of set
bounds consistency is commonly used. A domain $D$ is \emph{bounds
  consistent} for a constraint $c$ if for every variable $v \in
\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$. As
above, a set of propagators $F$ is said to maintain set bounds
consistency for a constraint $c$ if $\solv(F,D)$ is bounds consistent
for all domains $D$.

\begin{definition}
\label{def:setb-propagator}
A \emph{set bounds propagator} for a constraint $c$ is a function
$\setb(c)$ mapping domains to domains satisfying the following identity:
\begin{equation*}
  \setb(c)(D)(v) =
  \begin{cases}
    \conv(\dom(c)(\ran(D))(v)) & \text{if } v \in \vars(c) \\
    D(v) & \text{ otherwise}
  \end{cases}
\end{equation*}
\end{definition}

\begin{lemma}
  A set bounds propagator $\setb(c)$ for a constraint $c$ is correct,
  checking, and idempotent. A propagator $\setb(c)$ is also monotonic on
  range domains.
\end{lemma}
\begin{proof}
  All of these follow trivially from the properties of $\dom(c)$ and
  from the extension and idempotence properties of $\conv$ and $\ran$.
\end{proof}

Clearly for any constraint $c$ and any $v \in \vars(c)$ we have
$\dom(c)(D)(v) \subseteq \setb(c)(D)(v)$ for all domain propagators
$\dom(c)$ and all set bounds propagators $\setb(c)$.

%% We can compare the propagation behaviour of two
%% sets of propagators. We say a set of propagators $F_1$ is
%% \emph{stronger} than a set of propagators $F_2$, written $F_1 \gg
%% F_2$, if $solv(F_1, D) \sqsubset solv(F_2, D)$ for all domains $D$. A
%% set of propagators $F_1$ is \emph{equivalent} to a set of propagators
%% $F_2$, written $F_1 \approx F_2$, if $solv(F_1, D) = solv(F_2,D)$ for
%% all domains $D$.



\subsection{Boolean Formul\ae{} and Binary Decision Diagrams}

We use Boolean formul\ae{} extensively to model sets and set relations.
In particular, we make use of the following Boolean operations: $\wedge$
(conjunction), $\vee$ (disjunction), $\neg$ (negation), $\ra$
(implication), $\leftrightarrow$ (bi-directional implication), $\oplus$
(exclusive OR), and $\exists$ (existential quantification).  We use the
shorthand $\exists_V F$ for $\exists_{x_1} \cdots \exists_{x_n} F$ where
$V = \{x_1, \ldots, x_n\}$, and we use $\bexists_V F$ to mean
$\exists_{V'} F$ where $V' = \vars(F) \setminus V$.

Binary Decision Trees (BDTs) are a well-known method of modeling Boolean
functions on Boolean variables. A Binary Decision Tree is a complete
binary tree, in which each internal node $n(v,t,f)$ is labelled with a
Boolean variable $v$ and each leaf node is labelled with a truth value 0
(false) or 1 (true).  Each internal node corresponds to an if-then-else
test of the labelled variable, and has two outgoing arcs---the ``false''
arc (to BDT $f$) and the ``true'' arc (to BDT $t$).  In order to
evaluate the function represented by the binary tree, the tree is
traversed from the root to a leaf. On reaching any internal node, the
value of the Boolean variable whose label appears on the node is
examined, and the corresponding arc is followed. The traversal stops
upon reaching a leaf node, whereupon the value of the function is taken
to be the label of that node.

A Binary Decision Diagram (BDD) is a variant of a Binary Decision Tree,
formed by relaxing the tree structure requirement, instead representing
functions as directed acyclic graphs. In a Binary Decision Diagram, a
node is permitted to have multiple parents, as opposed to a tree
structure which requires any node to have at most one parent. This
effectively permits common portions of the tree to be shared between
multiple branches, allowing a compact representation of many functions.

The primary disadvantage of both Binary Decision Trees and Binary
Decision Diagrams is that the representation of a given function is not
canonical (i.e. one function may have many representations).  Two
additional canonicity properties allow many operations on BDDs to be
performed efficiently \cite{bryant:graph,bryant:robdds}. A BDD is said
to be \emph{reduced} if it contains no identical nodes (that is, there
are no nodes with the same label and identical ``then'' and ``else''
arcs), and no redundant tests (there are no nodes with both ``then'' and
``else'' arcs to the same node). A BDD is said to be \emph{ordered} if
there is a total order $\prec$ of the variables, such that if there is
an arc from a node labelled $v_1$ to a node labelled $v_2$ then $v_1
\prec v_2$. A \emph{Reduced Ordered BDD} (ROBDD) is canonical function
representation up to reordering, which permits an efficient
implementation of many Boolean function operations. For more details the
reader is referred to the work of \citeA{bryant:robdds}, or the
introduction of \citeA{andersen:bddnotes}.

We define the size $|R|$ to be the number of non-leaf nodes of the ROBDD
$R$, as well as defining $\VAR(R)$ to be the set of ROBDD variables that
appear as labels on the internal nodes of $R$.  We shall be interested
in \emph{stick} ROBDDs, in which every internal node $n(v,t,f)$ has
exactly one of $t$ or $f$ as the constant $0$ node.

\begin{figure}[t]
\begin{center}
\begin{tabular}{ccc}
    $\xymatrix@=5pt{
                                %      \xyo{v_2}\ar@{-->}@/_35pt/[dddddd]\ar@{->}[dr] \\
      & \xyo{v_3}\ar@{-->}@/_25pt/[dddddl]\ar@{->}[dr] \\
      & & \xyo{v_4}\ar@{->}@(dr,r)@/^25pt/[ddddll]\ar@{-->}[dl] \\
      & \xyo{v_5}\ar@{->}@(dr,r)@/^20pt/[dddl]\ar@{-->}[dl] \\
      \xyo{v_6}\ar@{-->}@/_15pt/[dd]\ar@{->}[dr] \\
      & \xyo{v_7}\ar@{-->}[dl]\ar@{->}[dr] \\
      *+[F]{\mathrm{0}}  & & *+[F]{\mathrm{1}}
    }$
    &
    $\xymatrix@=5pt{
      &&& \xyo{v_1}\ar@{-->}[dll]\ar@{->}[drr] \\
      & \xyo{v_2}\ar@{-->}[dl]\ar@{->}[dr] &&&& \xyo{v_2}\ar@{-->}[dl]\ar@{->}[dr] \\
      \xyo{v_8}\ar@{-->}@(dl,l)@/_25pt/[ddrrr]\ar@{->}[dr]
      &&
      \xyo{v_8}\ar@{-->}[dl]\ar@{->}[ddr]
      &&
      \xyo{v_8}\ar@{-->}[ddl]\ar@{->}[dr]
      &&
      \xyo{v_8}\ar@{->}@(dr,r)@/^25pt/[ddlll]\ar@{-->}[dl]
      \\
      & \xyo{v_9}\ar@{-->}@(dl,l)@/_15pt/[drr]\ar@{->}[ddrr]
      &&&&\xyo{v_9}\ar@{->}@(dr,r)@/^15pt/[dll]\ar@{-->}[ddll]
      \\
      &&&  *+[F]{\mathrm{0}} \\
      &&&  *+[F]{\mathrm{1}} \\
    }$
    &
    $\xymatrix@=2pt{
      &&&& \xyo{v_1}\ar@{-->}[dll]\ar@{->}[drr] \\
      && \xyo{v_2}\ar@{-->}[dl]\ar@{->}[dr] &&&&
      \xyo{v_2}\ar@{-->}[dl]\ar@{->}[dr] \\
      & \xyo{v_3}\ar@{->}[dr] && \xyo{v_3}\ar@{->}[dr] &&
      \xyo{v_3}\ar@{->}[dr] && \xyo{v_3}\ar@{->}[dr] \\
      && \xyo{v_4}\ar@{-->}[dl] && \xyo{v_4}\ar@{-->}[dl] && \xyo{v_4}\ar@{-->}[dl] && \xyo{v_4}\ar@{-->}[dl] \\
      & \xyo{v_5}\ar@{-->}[dl] && \xyo{v_5}\ar@{-->}[dl] &&
      \xyo{v_5}\ar@{-->}[dl] && \xyo{v_5}\ar@{-->}[dl]  \\
      \xyo{v_6}\ar@{->}[dr] && \xyo{v_6}\ar@{->}[dr] && \xyo{v_6}\ar@{->}[dr]
      && \xyo{v_6}\ar@{->}[dr] \\
      & \xyo{v_7}\ar@{->}[dr] && \xyo{v_7}\ar@{->}[dr] && \xyo{v_7}\ar@{->}[dr]
      && \xyo{v_7}\ar@{->}[dr]  \\
      && \xyo{v_8}\ar@{->}[dr]
      &&
      \xyo{v_8}\ar@{-->}[dl]
      &&
      \xyo{v_8}\ar@{->}[dr]
      &&
      \xyo{v_8}\ar@{-->}[dl]
      \\
      &&& \xyo{v_9}\ar@{->}[drr]
      &&&&\xyo{v_9}\ar@{-->}[dll]
      \\
      &&&&&  *+[F]{\mathrm{1}}
    }$
    \\
    (a) & (b) & (c)
  \end{tabular}
\end{center}
\caption{ROBDDs for (a)
    $\lu{} = v_3 \wedge \neg v_4 \wedge \neg v_5 \wedge v_6 \wedge v_7$
    (b) $R = \neg (v_1 \leftrightarrow v_9) \wedge \neg (v_2
    \leftrightarrow v_8)$ and (c) $\lu \wedge R$ (omitting the node $0$
    and arcs to it). Solid arcs are ``then'' arcs, dashed arcs are
    ``else'' arcs.  }
  \label{fig:g}
  \rule{\textwidth}{0.4pt}\\
\end{figure}

\begin{example}
  Figure~\ref{fig:g}(a) gives an example of a stick ROBDD $\lu$ representing
  the formula $v_3 \wedge \neg v_4 \wedge \neg v_5 \wedge v_6 \wedge
  v_7$.
  $|\lu| = 5$ and $VAR(\lu) = \{v_3,v_4,v_5,v_6,v_7\}$.  
  Figure~\ref{fig:g}(b) gives an example of a more complex ROBDD
  representing the formula $\neg (v_1 \leftrightarrow v_9) \wedge
  \neg (v_2 \leftrightarrow v_8)$.  
  $|R| = 9$ and $VAR(R) = \{v_1,v_2,v_8,v_9\}$.
  One can verify that the valuation
  $\{v_1 \mapsto 1, v_2 \mapsto 0, v_8 \mapsto 1, v_9 \mapsto 0\}$
  makes the formula true by following the path right, left, right, left
  from the root. 
  \ignore{\ph{Possibly change this example}}
\end{example}

\subsection{ROBDD Operations}


There are efficient algorithms for many Boolean operations applied to
ROBDDs. The complexity of the basic operations for constructing new
ROBDDs is $O(|R_1||R_2|)$ for $R_1 \wedge R_2$, $R_1 \vee R_2$ and $R_1
\leftrightarrow R_2$, $O(|R|)$ for $\neg R$, and $O(|R|^2)$ for $\exists
v ~R$.  Note however that we can test whether two ROBDDs are identical,
whether $R_1 \leftrightarrow R_2$ is equivalent to true (1), in $O(1)$.
 
We give code for conjunction $R_1 \wedge R_2 = \textsf{and}(R_1,R_2)$
and existential quantification (of one variable) $\exists v ~R =
\textsf{exists}(v,R)$. The code for disjunction $R_1 \vee R_2 =
\textsf{or}(R_1,R_2)$ is dual to \textsf{and} and is very similar in
structure.
% ph: This actually isn't true. You need to change some ors to ands as
% well. Probably better to leave things unspecified than labouring the
% point.
%, obtained by replacing
%each 0 by 1 and 1 by 0
The code make use of the auxiliary function
\textsf{node} which builds a new ROBDD node. The \textsf{node}
function returns $t$ if $t = f$,
and is in practice is memoed so any call to \textsf{node} with the same
arguments as a previous call 
returns a reference to the previously created ROBDD node.
\begin{figure}[t]
\begin{quote}
\begin{tabbing}
xx \= xx \= xx \= xx \= xx \= \kill
\textsf{node}($v,t,f$) \textbf{if} ($t = f$) \textbf{return} $t$ \textbf{else}
\textbf{return} $n(v,t,f)$ \\
\\
\textsf{and}($R_1,R_2$) \\
\> \textbf{if} ($R_1 = 0$ \textbf{or} $R_2 = 0$) \textbf{return} 0 \\
\> \textbf{if} ($R_1 = 1)$ \textbf{return} $R_2$ \\
%\> \textbf{if} ($R_2 = 0)$ \textbf{return} 0 \\
\> \textbf{if} ($R_2 = 1)$ \textbf{return} $R_1$ \\
\> $n(v_1,t_1,f_1)$ := $R_1$ \\
\> $n(v_2,t_2,f_2)$ := $R_2$ \\
\> \textbf{if} ($v_1 \preceq v_2$)  \textbf{return} 
      \textsf{node}($v_1$,\textsf{and}($t_1,R_2$),\textsf{and}($f_1,R_2$))\\
\> \textbf{else if} ($v_2 \preceq v_1$)  \textbf{return} 
      \textsf{node}($v_2$,\textsf{and}($t_2,R_1$),\textsf{and}($f_2,R_1$))\\
\> \textbf{else}   \textbf{return} 
      \textsf{node}($v_1$,\textsf{and}($t_1,t_2$),\textsf{and}($f_1,f_2$))\\
%% \\
%% \textsf{or}($R_1,R_2$) \\
%% \> \textbf{if} ($R_1 = 1$ \textbf{or} $R_2 = 1$) \textbf{return} 1 \\
%% \> \textbf{if} ($R_1 = 0)$ \textbf{return} $R_2$ \\
%% \> \textbf{if} ($R_2 = 0)$ \textbf{return} $R_1$ \\
%% \> \textbf{let} $R_1 = n(v_1,t_1,f_1)$ and $R_2 = n(v_2,t_2,f_2)$ \\
%% \> \textbf{if} ($v_1 \preceq v_2$)  \textbf{return} 
%%       \textsf{node}($v_1$,\textsf{or}($t_1,R_2$),\textsf{or}($f_1,R_2$))\\
%% \> \textbf{else if} ($v_2 \preceq v_1$)  \textbf{return} 
%%       \textsf{node}($v_2$,\textsf{or}($t_2,R_1$),\textsf{or}($f_2,R_1$))\\
%% \> \textbf{else}   \textbf{return} 
%%       \textsf{node}($v_1$,\textsf{or}($t_1,t_2$),\textsf{or}($f_1,f_2$))\\
\\
\textsf{exists}($v, R$) \\
\> \textbf{if} ($R = 0)$ \textbf{return} 0 \\
\> \textbf{if} ($R = 1)$ \textbf{return} 1 \\
\> $n(v_r,t,f)$ := $R$ \\
\> \textbf{if} ($v_r \preceq v$)  \textbf{return} 
      \textsf{node}($v_r$,\textsf{exists}($v,t$),\textsf{exists}($v,f$))\\
\> \textbf{else if} ($v \preceq v_r$)  \textbf{return} $R$ \\
\> \textbf{else}   \textbf{return} 
      \textsf{or}($t,f$) 
\end{tabbing}
\end{quote}
\caption{Example ROBDD operations}
\label{fig:ops}
\end{figure}

Modern BDD packages provide many more operations, including specialized
implementations of some operations for improved speed.  Some important
operations for our purposes are existential quantification of multiple
variables $V$
for a formula $R$ $\exists V~ R$, and the combination of conjunction and
existential quantification $\exists V~ R_1 \wedge R_2$.

Although in theory the number of nodes in an ROBDDs can be exponential
in the number of variables in the represented Boolean function, in
practice ROBDDs are often very compact and computationally efficient.
This is due to the fact that ROBDDs exploit to a very high-degree the
symmetry of models of a Boolean formula.


\section{Modeling Set CSPs Using ROBDDs}
\label{sect:modeling}

In this section we discuss how to solve set constraint satisfaction
problems using ROBDDs. There are three parts to this---the modeling of
set domains as ROBDDs, the modeling of set constraints as ROBDDs, and
how to use both to produce a set solver.

\subsection{Modeling Set Domains using ROBDDs}

Suppose the universe $\mathcal{U} = \{1, \ldots, N\}$. We assume that
all set values are subsets of $\mathcal{U}$. Let $x$ be a set
variable, and let $D$ be a domain over $\mathcal{V}$. Given that the size
of the universe is bounded, we can
associate a Boolean variable $x_i$ with each potential element $i
\in \{1, \ldots, N\}$ of $x$. 
%The value of each $x_i$ then corresponds to the state
%of the $i$th bit in the characteristic vector $\chr(\theta(x))$ of a
%valuation $\theta$ of $x$.
%\pjs{never use characteristic vector?}
Hence a set variable $x$ is represented by a vector of Boolean
variables $V(x) = \n{x_1, \ldots, x_N}$.


Take any set $A \in D(x)$. Then we can represent $A$ as a valuation on
the variables $\n{x_1, \ldots, x_N}$ defined as $\theta_A = \{ x_1
\mapsto (1 \in A), \ldots, x_n \mapsto (n \in A) \}$.  We can represent
this valuation $\theta_A$ and consequently the set $A$ \ignore{\pjs{Its
    important that we represent sets by valuations not formul\ae!}}  as
a Boolean formula $B(A)$ which has this valuation as its unique
solution:
\[
B(A) = \bigwedge_{i \in \mathcal{U}} y_i \text{ where } y_i = 
\begin{cases}
x_i & \text{ if } i \in A
\\
\neg x_i & \text{ otherwise}
\end{cases}
\]

Given that we can represent any element of $D(x)$ as a valuation, we
can represent $D(x)$ itself as a Boolean formula $B(D(x))$ 
whose solutions $\theta_A$ correspond to the elements $A \in D(x)$.
That is, $B(D(x))$ is the disjunction of $B(A)$ over all 
possible sets $A \in D(x)$:
\[
B(D(x)) = \bigvee_{A \in D(x)} B(A) \text{ where } B(A) \text{ is defined
above
}
\]

Any solution $\theta(x) \in D(x)$ then corresponds to a satisfying
assignment of the Boolean formula $B(D(x))$.
From now on we will overload the notion of a domain $D$
to equivalently return a set of possible sets $D(x)$ 
for a variable $x$, or its equivalent Boolean representation $B(D(x))$.

\begin{figure}[t]
\centering
  \begin{center}
\begin{tabular}{ccc}
\xymatrix@=10pt{
        &\xyo{v_1}\ar@{-->}[dl]\ar[drr]\\
        \xyo{v_2}\ar[dr]\ar@{-->}@(dl,l)[dddr] &&&
        \xyo{v_2}\ar@{-->}[ddl]\ar@/^2.8pc/[dddll] \\
        &\xyo{v_3}\ar[dr]\ar@{-->}@(dl,ul)[dd]\\
        &&\xytrue\\
        &\xyfalse
      }
&
~~~~~~~~~
    $\xymatrix@=10pt{
      & \xyo{x_1}\ar@{-->}@/_25pt/[dddd]\ar@{->}[dr]\\
      && \xyo{x_2}\ar@{->}@(dr,r)@/^35pt/[dddl]\ar@{-->}[dl] \\
      & \xyo{x_3}\ar@{->}[dr]\ar@{-->}@(dl,l)[dd] \\
      &&\xyo{x_4}\ar@{->}[dr]\ar@{-->}[dl] \\
      & *+[F]{\mathrm{\False}}
      && *+[F]{\mathrm{\True}}
    }$
~~~~~~~~~ & ~~~~
    $\xymatrix@=10pt{
      &&& \xyo{s_1}\ar@{->}[dr]\ar@{-->}[dl]\\
      && \xyo{s_2}\ar@{->}[dr]\ar@{-->}[dl] &&
      \xyo{s_2}\ar@{->}[dr]\ar@{-->}[dl]\\
      &\xyo{s_3}\ar@{->}[dr]
      \ar@{-->}[dl]
      &&\xyo{s_3}\ar@{->}[dr]\ar@{-->}[dl]
      &&\xyo{s_3}\ar@{-->}[dl]\ar@{->}@/^35pt/[ddddlll] \\
      \xyo{s_4}\ar@{->}[dr]\ar@{-->}@/_35pt/[dddrr]
      && \xyo{s_4}\ar@{->}[dr]\ar@{-->}[dl]
      && \xyo{s_4}\ar@{-->}[dl]\ar@{->}@(dr,r)@/^25pt/[dddll]   &&\\
 & \xyo{s_5}\ar@{->}[dr]\ar@{-->}@/_25pt/[ddr]
      &&
\xyo{s_5}\ar@{-->}[dl]\ar@{->}@(dr,r)@/^15pt/[ddl]
&&&&&\\
      && *+[F]{\mathrm{\True}}         &&&   \\
      && *+[F]{\mathrm{\False}} &&&&
    }$

\\
(a) & (b) & (c)
\end{tabular}
  \end{center}
  \caption{ROBDDs for (a) $\{\{1\},\{1,3\},\{2,3\}\}$, 
(b) $\range{\{1,3,4,5\}}{\{1,3,4,5,6,\ldots N\}}$
 and (c) $\{s\subseteq\s{1,2,3,4,5} ~|~ \card{s}{2}\}.$}
  \label{fig:doms}
  \rule{\textwidth}{0.4pt}\\
\end{figure}

\begin{example}
Let $\mathcal{U} = \{1,2,3\}$. Suppose $v$ is a set variable with 
$D(v) = \{\{1\}, \{1,3\},\{2,3\}\}$. 
We associate Boolean variables $\{v_1, v_2, v_3\}$ with $v$. 
We can then represent $D(v)$ as the
Boolean formula $(v_1 \wedge \neg v_2 \wedge \neg v_3) \vee 
(v_1 \wedge \neg v_2 \wedge v_3) \vee (\neg v_1
\wedge v_2 \wedge v_3)$. The three solutions to this formula correspond
to the elements of $D(v)$. The corresponding ROBDD $B(D(v))$ 
is shown in Figure~\ref{fig:doms}(a).
\end{example}

This representation is useful since such Boolean formul\ae{} can be
directly represented as reasonably compact ROBDDs. Given a Boolean
formula representing a domain, we can clearly construct the
corresponding ROBDD in a bottom-up fashion, but in practice we only ever
construct the ROBDD for a domain implicitly through constraint
propagation. As we will show, ROBDDs permit a surprisingly compact
representation of (many) subsets of $\mathcal{P(U)}$.

Since ROBDDs are ordered, we must specify a variable ordering for the
Boolean variables that we use. We arbitrarily order the ROBDD variables
for a single set variable $x$ as follows: $x_1 \prec x_2 \prec \ldots
\prec x_n$. Assuming there is no special relationship between the
elements of the universe, the choice of the ordering of ROBDD variables
that represent a particular set variable is unimportant. However, the
relative ordering of ROBDD variables comprising different set variables
has a drastic effect on the size of the ROBDDs representing constraints,
and is discussed in Section~\ref{sect:modeling-constraints}.

\begin{example}
  An ROBDD representing the $2^{N-5}$ sets in the interval
  $\range{\{1,3,4,5\}}{\{1,3,4,5,6, \ldots N\}}$ is shown in
  Figure~\ref{fig:doms}(b).
\end{example}


The ROBDD representation is flexible enough to be able to represent
compactly a wide variety of domains, even those that might at first seem
difficult to represent. For example, the set of all subsets of
$\{1,2,3,4,5\}$ of cardinality 2 can be represented by the ROBDD shown
in Figure~\ref{fig:doms}(c).

For convenience, we set all of the initial variable domains $D_{init}(x)
= \mathcal{P(U)}$, which corresponds to the constant ``1'' ROBDD.
Restrictions on the initial bounds of a domain can instead be expressed
as unary constraints.

\subsection{Modeling Primitive Set Constraints Using ROBDDs}
\label{sect:modeling-constraints}

A major benefit of using ROBDDs to model set constraint problems is that
ROBDDs can be used to model the constraints themselves, not just the set
domains. Any set constraint $c$ can be converted to a Boolean formula
$B(c)$ on the Boolean variables comprising the set variables $\vars(c)$,
such that $B(c)$ is satisfied if and only if the corresponding set
variable valuations satisfy $c$. As usual, we represent $B(c)$ as an
ROBDD.

\begin{example}
  \label{example:subset}
  Let $v$ and $w$ be set variables over the universe
  $\mathcal{U}=\{1,2,3\}$, and let $c$ denote the constraint $v
  \subseteq w$. Assume that the Boolean variables associated with $v$
  and $w$ are $\n{v_1, v_2, v_3}$ and $\n{w_1, w_2, w_3}$ respectively.
  Then $c$ can be represented by the Boolean formula $(v_1 \ra w_1)
  \wedge (v_2 \ra w_2) \wedge (v_3 \ra w_3)$. This formula can in turn
  be represented by either of the ROBDDs shown in
  Figure~\ref{fig:subset-bdd} (depending on the variable order).
\end{example}

\begin{figure}[t]
  $$
  \begin{array}{ccc}

    \xymatrix@=10pt{
    &&\xyo{v_1}\ar@{-->}[dl]\ar@{->}[dr]\\
    &\xyo{v_2}\ar@{-->}[dl]\ar@{->}[d]&&\xyo{v_2}\ar@{->}[d]\ar@{-->}[dr]\\
    \xyo{v_3}\ar@/_20pt/@{-->}[dddddrr]\ar@/_20pt/[ddddrrr]&\xyo{v_3}\ar@/_10pt/@{-->}[dddr]\ar@/_10pt/@{->}[dddrr]&&\xyo{v_3}\ar@{-->}[dl]\ar@{->}[d]&\xyo{v_3}\ar@{->}[d]\ar@{-->}[dr]\\
    &&\xyo{w_1}\ar@{->}[dd]&\xyo{w_1}\ar@{->}[dd]&\xyo{w_1}\ar@/^10pt/@{->}[dddl]&\xyo{w_1}\ar@/^20pt/@{->}[ddddlll]\\
    \\
    &&\xyo{w_2}\ar@{->}[dd]&\xyo{w_2}\ar@{->}[d]\\
    &&&\xyo{w_3}\ar@{->}[dl]\\
    &&*+[F]{\mathrm{\True}}
    }

\ignore{    \xymatrix@=10pt{
      &&\xyo{v_1}\ar@{-->}[dl]\ar@{->}[dr]\\
      &\xyo{v_2}
      \ar@{-->}@(dl,l)[dddrrr]% `dl^d[dl] `^dr[ddl] [dddr]
      \ar@{->}[ddrrrr]
      && \xyo{v_2}\ar@{-->}[dl]\ar@{->}[dr]\\
      &&\xyo{w_1}\ar@{-->}@(dl,ul)[dd]\ar@{->}[ddrr] &&
      \xyo{w_1}\ar@{-->}[ddll]\ar[dr]\\
      &&&&&\xyo{w_2}\ar@{-->}[dlll]\ar@(dr,r)[dl]\\
      &&*+[F]{\mathrm{\False}} && *+[F]{\mathrm{\True}} }
      }
    &~~~~~~~~&
    \ignore{
    \xymatrix@=10pt{
      \xyo{v_1}\ar[dr]\ar@/_20pt/@{-->}[ddrr] \\
      &\xyo{w_1} \ar[dr]\ar@/_20pt/@{-->}[ddd]\\
      && \xyo{v_2} \ar[dr]\ar@/_20pt/@{-->}[ddr]\\
      &&& \xyo{w2}\ar@/^15pt/[d]\ar@{-->}@(dl,r)[dll]\\
      &*+[F]{\mathrm{\False}} && *+[F]{\mathrm{\True}}
    }
    }
    \xymatrix@=10pt{
    \xyo{v_1}\ar@{->}[dr]\ar@/_20pt/@{-->}[ddrr]\\
    &\xyo{w_1}\ar@{->}[dr]\ar@/_20pt/@{-->}[dddddr]\\
    &&\xyo{v_2}\ar@{->}[dr]\ar@/_20pt/@{-->}[ddrr]\\
    &&&\xyo{w_2}\ar@{->}[dr]\ar@/^10pt/@{-->}[dddl]\\
    &&&&\xyo{v_3}\ar@{->}[dr]\ar@/_20pt/@{-->}[ddr]\\
    &&&&&\xyo{w_3}\ar@/^10pt/@{->}[d]\ar@/^5pt/@{-->}[dlll]\\
    &&*+[F]{\mathrm{\False}} &&& *+[F]{\mathrm{\True}}
    }
    \\
    ~
    \\
    \text{(a) Order: } v_1\prec v_2 \prec v_3\prec w_1\prec w_2 \prec
    w_3&&
    \text{(b) Order: } v_1\prec w_1\prec v_2\prec w_2 \prec v_3 \prec
    w_3\\
    \text{\scriptsize (Arcs to `0' are not shown for clarity)}&&
  \end{array}
  $$
  \caption{Two ROBDDs for $v\subseteq w$}
  \label{fig:subset-bdd}
  \rule{\textwidth}{0.4pt}\\
\end{figure}

Example~\ref{example:varorder} demonstrates that the order of the
variables within the ROBDDs has a very great effect on the size of the
formula representations.

\begin{example}
    \label{example:varorder}
    Consider again the constraint $v \subseteq w$ as described in
    Example~\ref{example:subset}. Figure~\ref{fig:subset-bdd} shows
    the effect of two different variable orderings on the size of the
    resulting BDD. In this case, the variable ordering $v_1 \prec w_1
    \prec v_2 \prec w_2$ gives a much more compact representation of
    constraint than the ordering $v_1 \prec v_2 \prec w_1 \prec w_2$.
    As can be seen from the figure, the latter ordering
    has size exponential in $n$, whereas the former has size linear in
    $n$.

\end{example}

If the variable set $\mathcal{V} = \{v_1, \ldots, v_m\}$ has
corresponding Boolean variables $\n{v_{1,1}, v_{1,2}, \ldots, v_{1,N}}$,
$\n{v_{2,1}, v_{2,2}, \ldots, v_{2,N}}$, \ldots, $\n{v_{m,1}, v_{m,2},
  \ldots, v_{m,N}}$, we choose to order the Boolean variables $v_{1,1}
\prec \cdots \prec v_{m,1} \prec v_{1,2} \prec \cdots \prec v_{m,2}
\prec \cdots \prec v_{1,N} \prec \cdots \prec v_{m,N}$. This ordering
guarantees a linear representation for all of the primitive set
constraints except cardinality.  The reason is that primitive set
constraints except cardinality are defined elementwise, that is an
element $i \in v$ never constrains whether $j \in w$ or $j \notin w$
for $i \neq j$.  For this reason there are no interactions between bits
$v_i$ and $w_j$, $i \neq j$. Placing all the bits for the same element
adjacent in the order means that the interactions of bits $v_i$ and
$w_i$ are captured in a small ROBDD, which is effectively separate from
the ROBDD
describing the interactions of the next 
element's bits.  Table~\ref{tab:size} contains a list of the primitive set
constraints, their corresponding Boolean formul\ae{} and the sizes of
the ROBDDs representing those formul\ae{} under this point-wise ordering.


Cardinality constraints can be represented in a quadratic number of
ROBDD nodes using a simple recursive definition. We define
$\cardlu(\n{v_{i_1},\ldots,v_{i_n}},l,u)$ to be the Boolean formula
which restricts the number of true bits in the vector
$\n{v_{i_1},\ldots,v_{i_n}}$ to between $l$ and $u$ inclusive.
\[
\cardlu(\n{v_{i_1},\ldots,v_{i_n}},l,u) =
\begin{cases}
  \True &  \text{if } l \leq 0 \wedge n \leq u\\
  \False & \text{if } n < l \vee u < 0\\
  \begin{array}{@{}l}
    (\neg v_{i_1} \wedge \cardlu(\n{v_{i_2},\ldots,v_{i_n}},l,u))
    ~~\vee\\
    (v_{i_1} \wedge \cardlu(\n{v_{i_2},\ldots,v_{i_n}},l-1,u-1)
  \end{array}
  & \text{otherwise}
\end{cases}
\]
%% \[
%% \cardlu(v,l,u,n) =
%% \begin{cases}
%%   \True &  \text{if } l \leq 0 \wedge N - n \leq u\\
%%   \False & \text{if } N - n < l \vee u < 0\\
%% %  \neg v_n \wedge \cardlu(v,l,u,n+1)
%% %  & \text{if } n > 0 \wedge u=0 \\
%% %   v_n \wedge \cardlu(v,l-1,u-1,n+1)
%% %  & \text{if } n > 0 \wedge l=n \\
%%   \begin{array}{@{}l}
%%     (\neg v_{n+1} \wedge \cardlu(v,l,u,n+1))
%%     ~~\vee\\
%%     (v_{n+1} \wedge \cardlu(v,l-1,u-1,n+1)
%%   \end{array}
%%   & \text{otherwise}
%% \end{cases}
%% \]
It is clear from the structure 
of $\cardlu(\n{v_{i_1},\ldots,v_{i_n}},l,u)$ 
that 
the resulting ROBDD is $O(n^2)$ in size. A more general
method of
characterising cardinality constraints will be presented in
Section~\ref{sect:wsum}.

\begin{table}[t]
\begin{center}
%\def\arraystretch{1.2}
  \begin{tabular}{l|l@{~~}|c}
    \multicolumn{1}{c|}{$c$} & ~Boolean expression $B(c)$ & ~size of ROBDD \\
    \hline
    $k \in v$ & $v_k$ & $O(1)$ \\
    $k \not\in v$ & $\neg v_k$ & $O(1)$ \\
    $u = d$ & ~$\bigwedge_{i \in d} u_i \wedge \bigwedge_{1 \leq i \leq N, i
    \not\in d} \neg u_i$
    & $O(N)$ \\
    $u = v$ & ~$\bigwedge_{1 \leq i \leq N} (u_i \leftrightarrow v_i$) & $O(N)$
\\
%%PJS    $u \subset v$ & ~$\bigwedge_{1 \leq i \leq N} (u_i \rightarrow v_i)$ & $O(N)$ \\
    $u \subseteq v$ & ~$\bigwedge_{1 \leq i \leq N} (u_i \rightarrow v_i)$ & $O(N)$ \\
    $u = v \cup w~~$ & ~$\bigwedge_{1 \leq i \leq N} (u_i \leftrightarrow (v_i \vee w_i))$
    & $O(N)$ \\
    $u = v \cap w$ & ~$\bigwedge_{1 \leq i \leq N} (u_i \leftrightarrow (v_i \wedge w_i))$
    & $O(N)$ \\
    $u = v - w$ & ~$\bigwedge_{1 \leq i \leq N} (u_i \leftrightarrow (v_i \wedge
 \neg w_i))$
    & $O(N)$ \\
    $u = \overline{v}$ & ~$\bigwedge_{1 \leq i \leq N} \neg (u_i \lra v_i)$ & $O(N)$
\\
    $u \neq v$ & ~$\bigvee_{1 \leq i \leq N} \neg (u_i \leftrightarrow v_i)$ & $
O(N)$ \\
    $\card{u}{k}$ & ~$\cardlu(V(u),k,k))$ & $O(k(N-k))$ \\
    $\cardl{u}{k}$ & ~$\cardlu(V(u),k,N)$ & $O(k(N-k))$ \\
    $\cardu{u}{k}$ & ~$\cardlu(V(u),0,k)$ & $O(k(N-k))$ \\
  \end{tabular}
\end{center}
\caption{Boolean representation of set constraints and the size of the corresponding ROBDD.
    \label{tab:size}}
  \rule{\textwidth}{0.4pt}\\
\end{table}


\subsection{A Basic Set Constraint Solver}
\label{sect:basic-set-solver}

We now show how to construct a simple set domain propagator $\dom(c)$
for a constraint $c$. If $\vars(c) = \{v_1, \dotsc,
v_n\}$, then we define a function $\dom(c)$ mapping domains to domains 
as follows:
\[
\dom(c)(D)(v_i) =
\begin{cases}
    \bexists_{V(v_i)} B(c) \wedge \bigwedge_{i=1}^n D(v_i) & \text{ if
    } v_i \in \vars(c) \\
    D(v_i) & \text{ otherwise}
\end{cases}
\]

In other words, to perform propagation we take the conjunction of the
Boolean representations of the current domains of the variables in
$vars(c)$ with the Boolean representation of the constraint $B(c)$ and
project the result onto the Boolean variables $V(v_i)$ representing
each variable $v_i$. Since $B(c)$ and all $D(v_i)$ are ROBDDs, this
formula can be implemented directly using ROBDD operations.

\begin{figure}[t]
\begin{center}
  \begin{tabular}{cccc}
   $\xymatrix@=10pt{
        &\xyo{w_1}\ar@{-->}[dl]\ar[dr]\\
        \xyo{w_2}\ar[rrrd]\ar@{-->}@(dl,l)[dddrr]
        && \xyo{w_2}\ar@{-->}[dl]\ar[dr]\\
        &\xyo{w_3}\ar[dr]\ar@{-->}@(dl,ul)[ddr] &&
        \xyo{w_3}\ar@{-->}[dl]\ar@(dr,r)[ddl]\\
        && \xytrue\\
        && \xyfalse
      }$
&
   $\xymatrix@=10pt{
        &\xyo{v_1}\ar@{-->}[dddddddl]\ar[dr]\\
        && \xyo{w_1}\ar[rd]\ar@{-->}[ddddddll] \\
        &&& \xyo{v_2}\ar@/^50pt/[dddddlll]\ar@{-->}[dl] \\
        && \xyo{w_2}\ar@{-->}[ddl]\ar[dr]\\
        &&&\xyo{v_3}\ar@{-->}[dl]\ar@(dr,r)[dddlll]\\
        &\xyo{w_3}\ar[ddrr]\ar@{-->}[ddl] &
        \xyo{w_3}\ar@{-->}[ddr]\ar[ddll]\\
        \\
        \xyfalse
        &&& \xytrue
        &&\\
      }$
&
$\xymatrix@=10pt{
    \xyo{v_1}\ar[dr]\ar@{-->}@(dl,l)@/_15pt/[ddrr] \\
    & \xyo{v_2}\ar@{-->}[dl]\ar[dr] \\
    \xytrue && \xyfalse
}$
~~~&~~~
$\xymatrix@=10pt{
  &\xyo{w_1}\ar[dr]\ar@{-->}@(dl,ul)[dddl] \\
&& \xyo{w_2}\ar[dr]\ar@{-->}[dl] \\
        &\xyo{w_3}\ar[dr]\ar@{-->}[dl] &&
        \xyo{w_3}\ar@{-->}[dl]\ar@(dr,r)@/^30pt/[dlll]\\
        \xyfalse
        && \xytrue
      }$
\\
~
\\
(a) & (b) & (c) & (d) 
\end{tabular}
\end{center}
\caption{ROBDDs used in the domain propagation
of $c \equiv v\subseteq w$
(a) $D(w)$, (b) $B(c) \wedge D(v) \wedge D(w)$, (c) $D'(v)$,
and (d) $D'(w)$}
  \label{fig:prop}
  \rule{\textwidth}{0.4pt}\\
\end{figure}

\begin{example}
Consider the domain propagation of constraint $v \subseteq w$
with initial domain 
$D(v) = \{ \{1\}, \{1,3\}, \{2,3\} \}$, 
and $D(w) = \{ \{2\}, \{1, 2\}, \{1, 3\} \}$ 
from Examples~\ref{ex:subset} and~\ref{ex:subset2}.
We conjoin the ROBDD $B(c)$ shown in Figure~\ref{fig:subset-bdd}(b)
with the domain ROBDD $D(v)$ shown in Figure~\ref{fig:doms}(a)
and the ROBDD for $D(w)$ shown in Figure~\ref{fig:prop}(a). 
The result is an ROBDD representing solutions of the formula
$v \subseteq w \wedge v \in D(v) \wedge w \in D(w)$
shown in Figure~\ref{fig:prop}(b).
We project the resulting onto $V(v)$ and $V(w)$, individually
obtaining the ROBDDs $D'(v) = \{\{1\},\{1,3\}\}$ 
and $D'(w) = \{\{1,2\},\{1,3\}\}$ shown in 
Figure~\ref{fig:prop}(c) and (d) respectively.
\end{example}

We need to verify the correctness of the propagator:

\begin{lemma}
    Let $\mathcal{V}$ be the collection of all set variables, and let $c$ be
    a set constraint with $\vars(c) = \{v_1, \dotsc, v_k\} \subseteq
    \mathcal{V}$. Then $\dom(c)$ is a domain propagator for $c$.
\end{lemma}
\begin{proof}
    We need to verify that $\dom(c)$ satisfies the identity of
    Definition~\ref{def:dom-propagator}. Suppose $D$ is a domain over
    $\mathcal{V}$. We need to check for each $v \in
    \vars(c)$ that $\left \{ \chr(\theta(v)) ~|~ \theta \in D \wedge \theta \in
    \solns(c) \right \} = \left \{ u ~|~ u \entails 
    \bexists_{V(v_i)} B(c) \wedge \bigwedge_{i=1}^n D(v_i) \right \}$, where
    $\chr(X)$ denotes the characteristic vector of $X$.
    This is clearly true, since the values $\theta \in \solns(c)$ are by
    definition the satisfying assignments of $B(c)$, and the values
    $\theta \in D$ are by definition the satisfying assignments of
    $\wedge_{i=1}^n D(v_i)$. Hence the equality holds, implying that
    $\dom(c)$ is a domain propagator.
\end{proof}

If $F$ is the set of all
such domain propagators, then we can define a complete propagation
algorithm $\solv(F,D)$:
For simplicity we use the set of constraints $C$ rather than the
corresponding propagators. 

%\begin{figure}[t]
\begin{quote}
\begin{tabbing}
xx \= xx \= xx \= xx \= xx \= \kill
\textsf{solv}($C,D$) \\
\> $Q$ := $C$ \\
\> \textbf{while} ($\exists ~ c \in Q$) \\
%\> \> $Q$ := $Q - \{c\}$ \\
\> \> $D'$ := $\dom(c)(D)$ \\
\> \> \textbf{if} ($|\vars(c)| = 1$) $C$ := $C - \{c\}$ \\ 
\> \> $V$ := $\{ v \in {\cal V} ~|~ D(v) \neq D'(v) \}$ \\
\> \> $Q$ := $(Q \cup \{ c' \in C ~|~ \vars(c') \cap V \neq \emptyset \}) - \{c\}$ \\
\> \> $D$ := $D'$ \\
%\> \textbf{endwhile} \\
\> \textbf{return} $D$
\end{tabbing}
\end{quote}
%\caption{An implementation of a propagation-based constraint solver
%using the set domain propagation operator}
%\label{fig:solver}
%  \rule{\textwidth}{0.4pt}\\
%\end{figure}
 We maintain a queue $Q$ of constraints to be
(re-)propagated, initially all $C$. We select a constraint $c$ from the
queue to propagate, calculating a new domain $D'$.  If the constraint is
unary we remove it from $C$, never to be considered again, since
all information is captured in the domain. We then determine which
variables $V$ have changed domain, and add to the queue constraints $c'
\in C$ involving these variables, with the exception of the current
constraint $c$.
% ph: I don't think this comment matters any more, since we just
% described all of the improvements anyway.
%The implementation of $\solv(F,D)$ shown above
%incorporates several improvements over a basic chaotic iteration scheme.

By combining this algorithm with the modeling techniques described
earlier, we have shown how to construct a simple ROBDD-based set
domain propagator. Various improvements to this basic scheme
will be discussed in Section~\ref{sect:improvements} 
and Section~\ref{sect:propagation}.

The reader may wonder whether in this section we have done anything more than
map set constraints to Boolean constraints and then apply a Boolean
ROBDD solver to them. This is \emph{not} the case. Crucially the
domain propagation is on the original set variables, not on the Boolean
variables that make them up. In fact the set bounds consistency 
approach~\cite{pecos,gervet:conjunto} can be considered as simply a mapping
of set constraints to Boolean constraints.

\section{Effective Modeling of Set Constraints Using ROBDDs}
\label{sect:improvements}

In this section we demonstrate that the ROBDD-based modeling approach
is very flexible, allowing us to produce highly efficient
implementations of many complex constraints. 

\subsection{Combining Constraints and Removing Intermediate Variables}

It is rarely possible to express all of the constraints in a real set
problem directly as primitive set constraints on the variables of the
original problem. Instead, we would usually like to express more
complicated set constraints, which can be decomposed into multiple primitive
set constraints, often requiring the introduction of intermediate variables.

\begin{example}
\label{example:decomposition}

Let $c$ be the constraint $|v \cap w| \leq k$, which requires that $v$
and $w$ have at most $k$ elements in common. This constraint is used
in modeling the problem of finding Steiner systems (see
Section~\ref{sect:exp-steiner}). Since this is not a primitive set
constraint, existing set solvers would usually implicitly decompose it into two
primitive set constraints and introduce an intermediate variable $u$.
The representation of the constraint then becomes $\exists u ~ u = v
\cap w \wedge |u| \leq k$.

\end{example}

In the case of Example~\ref{example:decomposition} the decomposition
does not affect the strength of the resulting propagator. To prove this
fact, we use two results presented by \citeA{choils:propredundancy}.
Choi, Lee, and Stuckey prove these results in the case of a finite
integer domain solver, although the proofs in the set domain case are
identical.  \footnote{The observation of Choi et al. that
  Lemma~\ref{lem:prop-conj-equiv} does not apply if the shared variable
  is a set variable is only true if we are performing set bounds
  propagation, not set domain propagation.}

\begin{lemma}[\citeR{choils:propredundancy}]
\label{lem:prop-conj-stronger}

    Let $c_1$ and $c_2$ be set constraints. Then 
$\solv(\{\dom(c_1 \wedge    c_2)\},D) \sqsubseteq \solv(\{\dom(c_1),
\dom(c_2) \}, D)$ for all domains $D$.

\end{lemma}
\begin{lemma}[\citeR{choils:propredundancy}]
  \label{lem:prop-conj-equiv}
  
  Let $c_1$ and $c_2$ be two set constraints sharing at most one
  variable $x \in \mathcal{V}$. Then $\solv(\{\dom(c_1), \dom(c_2)\},D)
  = \solv(\{\dom(c_1 \wedge c_2)\},D)$ for all domains $D$.

\end{lemma}

Even if the strength of the propagator is unaffected by the
decomposition, splitting a propagator introduces a new variable, thus
slowing down the propagation process. In cases where two set constraints
share more than one variable Lemma~\ref{lem:prop-conj-equiv} does not
apply, so in such a case there can also be a loss of propagation
strength.

The ROBDD representation of the constraints allows us to utilise such
complex constraints directly, thus avoiding the problems associated with
splitting the constraint. We can directly construct an ROBDD for complex
constraints by forming the conjunction of the corresponding primitive
constraints and existentially quantifying away the intermediate
variables.

\begin{figure}[t]
%\begin{center}
\hspace{-1.5cm}
\begin{tabular}{cc}
   $\xymatrix@=10pt{
      &&&&&& \xyo{v_1}\ar@{->}[dr]\ar@{-->}[ddll] \\
      &&&&&&& \xyo{w_1}\ar[dr]\ar@{-->}[dlll] \\
       &&&& \xyo{v_2}\ar@{->}[dr]\ar@{-->}[ddll] &&&&
      \xyo{v_2}\ar@{->}[dr]\ar@{-->}[ddll]\\
      &&&&&  \xyo{w_2}\ar@{->}[dr]\ar@{-->}[dlll] &&&&
      \xyo{w_2}\ar@{->}[dr]\ar@{-->}[dlll]\\
      &&\xyo{v_3}\ar@{-->}[ddddddrr]@(dl,l)\ar@{->}[dr]
      &&&&\xyo{v_3}\ar@{-->}[ddll]\ar@{->}[dr]
      &&&&\xyo{v_3}\ar@{-->}[ddll]\ar@{->}[dr] \\
      &&&\xyo{w_3}\ar@{-->}[dddddr]@(dl,l)\ar@{->}[dr]
      &&&&\xyo{w_3}\ar@{-->}[dlll]\ar@{->}[dr]
      &&&&\xyo{w_3}\ar@{-->}[dlll]\ar@{->}[dddddlll] \\
      &&&& \xyo{v_4}\ar@{->}[dr]\ar@{-->}[dddd]@(dl,ul)
      &&&& \xyo{v_4}\ar@{-->}[ddll]\ar@{->}[dr] \\
      &&&&& \xyo{w_4}\ar@{->}[dr]\ar@{-->}[dddl]@(dl,ul)
      &&&& \xyo{w_4}\ar@{-->}[dlll]\ar@{->}[dddl] \\
      &&&&&& \xyo{v_5}\ar@{-->}[ddll]\ar@{->}[dr] \\
      &&&&&&& \xyo{w_5}\ar@{-->}[dlll]\ar@{->}[dr] \\
      &&&& *+[F]{\mathrm{\True}} 
      &&&& *+[F]{\mathrm{\False}} &&&&&
    }$
&
    $\xymatrix@=10pt{
      &&& \xyo{u_1}\ar@{->}[dr]\ar@{-->}[dl]\\
      && \xyo{u_2}\ar@{->}[dr]\ar@{-->}[dl] &&
      \xyo{u_2}\ar@{->}[dr]\ar@{-->}[dl]\\
      &\xyo{u_3}\ar@{-->}[dddr]@(dl,l)\ar@{->}[dr]
      &&\xyo{u_3}\ar@{-->}[dl]\ar@{->}[dr]
      &&\xyo{u_3}\ar@{-->}[dl]\ar@{->}[dddl] \\
      && \xyo{u_4}\ar@{->}[dr]\ar@{-->}[dd]@(dl,ul)
      && \xyo{u_4}\ar@{-->}[dl]\ar@{->}[dd]   &&\\
 & 
      &&
\xyo{u_5}\ar@{-->}[dl]\ar@{->}[dr]
&&&&&\\
      && *+[F]{\mathrm{\True}}         &&  *+[F]{\mathrm{\False}} &&&&
    }$
\\
(a) & (b) \\
\end{tabular}
%\end{center}
\caption{ROBDDs for 
(a) $|v \cap w| \leq 2$ where $v,w \subseteq \{1,2,3,4,5\}$,
and (b) $|s| \leq 2$
 where $s \subseteq \{1,2,3,4,5\}$.}
  \label{fig:build}
  \rule{\textwidth}{0.4pt}\\
\end{figure}

\begin{example}
  Consider the constraint $c \equiv |v \cap w| \leq k$ as discussed in
  Example~\ref{example:decomposition}. We can build a domain propagator
  $\dom(c)$ directly for $c$ by constructing the ROBDD $\exists V(u) ~ u
  = v \cap w \wedge |u| \leq k$.  In this case the size of the 
resulting ROBDD is
  $O(kN)$.
  
  The ROBDD for $|v \cap w| \leq 2$, that is $\exists V(u) ~ B(u = v
  \cap w) \wedge B(|u| \leq 2)$, where $\mathcal{U} = \{1,2,3,4,5\}$ is
  shown in Figure~\ref{fig:build}(a). The ROBDD for $|u| \leq 2$ is
  shown in Figure~\ref{fig:build}(b) for comparison. Note how each $u_i$
  node in the Figure~\ref{fig:build}(b) 
  is replaced by $v_i \wedge w_i$ (the formula
  for $u_i$) in Figure~\ref{fig:build}(a).
\end{example}

\subsection{Modeling Global Constraints}
\label{sect:global}

As the previous section demonstrates, it is possible to join the ROBDDs
representing primitive set constraints into a single ROBDD representing
the conjunction of the constraints. We can use this to join large
numbers of primitive constraints to form \emph{global constraints},
which in many cases will improve performance due to stronger
propagation.

It is one of the strengths of the ROBDD-based modeling approach that it
is trivial to construct global constraints simply using ROBDD operations
on primitive set constraints, without laboriously writing code to
perform propagation on the global constraint. This approach is very
powerful, but is not feasible for all combinations of primitive
constraints. As we shall see, some global constraints that we might
desire to construct lead to ROBDDs that are exponentially large.

\begin{figure}[t]
 \[
 \begin{array}{ccc}
    \xymatrix@=10pt{
        &\xyo{x_1}\ar@{-->}[dl]\ar@{->}[dr]\\
        \xyo{y_1}\ar@/_15pt/@{-->}[d]\ar@{->}[drr]&&\xyo{y_1}\ar@/^15pt/@{-->}[d]\\
        \xyo{z_1}\ar@{->}[dr]&&\xyo{z_1}\ar@{-->}[dl]\\
        &\xyo{x_2}\ar@{-->}[dl]\ar@[->][dr]&\\
        \xyo{y_2}\ar@/_15pt/@{-->}[d]\ar@{->}[drr]&&\xyo{y_2}\ar@/^15pt/@{-->}[d]\\
        \xyo{z_2}\ar@{->}[dr]&&\xyo{z_2}\ar@{-->}[dl]\\
%        \vdots\ar@`{[d]}[dr]&\vdots&\vdots\ar@`{[d]}@{-->}[dl]\\
%        \ar[dr]&&\ar@{-->}[dl]\\
        &\xyo{x_3}\ar@{-->}[dl]\ar@{->}[dr]&\\
        \xyo{y_1}\ar@/_15pt/@{-->}[d]\ar@{->}[drr]&&\xyo{y_3}\ar@/^15pt/@{-->}[d]\\
        \xyo{z_3}\ar@{->}[dr]&&\xyo{z_3}\ar@{-->}[dl]\\
        & *+[F]{\mathrm{\True}}
    }
    &~~~~~~~~~&
%%     \xymatrix@=10pt{
%%         &&\xyo{x_1}\ar@{-->}[dl]\ar[dr]\\
%%         &\xyo{y_1}\ar@{-->}[dr]\ar@`{[dl],[ddl],[ddddl],[ddddddl]}[dddddd]&&\xyo{y_1}\ar[dl]\ar@{-->}@`{[dr],[ddr],[ddddr],[ddddddr]}[dddddd]\\
%%         &&\xyo{x_2}\ar@{-->}[dl]\ar[dr]&&\\
%%         &\xyo{y_2}\ar@{-->}[dr]\ar@/_13pt/[dddd]&&\xyo{y_2}\ar[dl]\ar@{-->}@/^13pt/[dddd]&\\
%%         &&\xyo{x_3}\ar@{-->}[dl]\ar[dr]&&\\
%%         &\xyo{y_3}\ar@{->}[dd]\ar@{-->}[ddrr]&&\xyo{y_3}\ar[ddll]\ar@{-->}[dd]&\\
%%         \\
%%         & *+[F]{\mathrm{\True}} && 
%%         *+[F]{\mathrm{\False}}& 
%%     }
   \xymatrix@=10pt{
        &&\xyo{x_1}\ar@{-->}[dl]\ar[dr]\\
        &\xyo{y_1}\ar@{-->}[dr]\ar@`{[dl],[ddl],[ddddl],[ddddddl]}[dddddd]&&\xyo{y_1}\ar[dl]\ar@{-->}@`{[dr],[ddr],[ddddddr]}[dddddd]\\
        &&\xyo{x_2}\ar@{-->}[dl]\ar[dr]&&\\
        &\xyo{y_2}\ar@{-->}[dr]\ar@/_13pt/[dddd]&&\xyo{y_2}\ar[dl]\ar@{-->}@/^13pt/[dddd]&\\
        &&\xyo{x_3}\ar@{-->}[dl]\ar[dddr]&&\\
        &\xyo{y_3}\ar@{->}[dd]\ar@{-->}[ddrr]\\
        \\
        & *+[F]{\mathrm{\True}} && 
        *+[F]{\mathrm{\False}}& 
    }
\\
\\
\partition(x,y,z)&&\lexlt(x,y) \\
 \text{\scriptsize (Arcs to `0' are omitted for clarity)}&&
\end{array}
\]

\caption{ROBDDs for the constraints $\partition(x,y,z)$ and
$\lexlt(x,y)$, where $x$, $y$
and $z$ are set variables taking values from the universe $\mathcal{U}
= \{1,2,3\}$.}
\label{fig:partition-and-lexlt}
\rule{\textwidth}{0.4pt}
\end{figure}
%\pjs{Suspicious that the lexlt bdd is for lexlteq??}

%\begin{example}
A useful global constraint is the constraint
$\partition(v_1,\dotsc,v_n)$, which requires that the sets $v_1, v_2,
\dotsc, v_n$ form a partition of the universe $\mathcal{U} =
\{1,\dotsc,N\}$. We can easily construct this constraint from primitive
set constraints as follows:
\begin{align*}
  \partition(v_1,\dotsc,v_n) =& \bigwedge_{i=1}^{n-1}
  \bigwedge_{j=i+1}^n \exists_{u_{ij}} (u_{ij}
  = v_i \cap v_j \wedge u_{ij} = \phi) \\
  \wedge \quad& \exists_{w_0} \cdots \exists_{w_n} ( w_0 = \phi \wedge
  (\bigwedge_{i=1}^n w_i = w_{i-1} \cup v_i) \wedge w_n = \mathcal{U})
\end{align*}

The propagator $\dom(\partition(\cdots))$ is stronger than a domain
propagator on the decomposition. For example, consider the constraint $c
\equiv \partition(x,y,z)$ as depicted in
Figure~\ref{fig:partition-and-lexlt}, and the domain $D$ where $D(x) =
D(y) = \{\{1\},\{2\}\}$ and $D(z) = \{\{2\}, \{3\}\}$.  Domain
propagation using the decomposition of the constraint will not alter
$D$, but domain propagation using the global constraint will give
$\dom(c)(D)(z) = \{\{3\}\}$. Hence stronger propagation can be gained
from global propagation using this constraint.
%\end{example}
 
Unfortunately not all global constraints can be modelled efficiently
using this approach. In particular, there is a risk that the ROBDD
representation of a global constraint could be extremely large, making
it infeasible to construct and use for propagation.

For example, consider the constraint $\atmost(\n{v_1,\dotsc,v_n}, k)$
proposed by \citeA{sadlerg:global}, which requires that each of the sets
$v_1, \dotsc, v_n$ has cardinality $k$ and the intersection of any pair
of the sets is at most 1 element in size. This constraint models the
Steiner triple systems of Section~\ref{sect:exp-steiner}.
\citeA{bessierehhw:globalsetcomplexity} proved that enforcing bounds
consistency on this constraint is NP-hard, so it follows that enforcing
domain consistency on this constraint is at least NP-hard as well.
Theoretically we can still construct an ROBDD representing this
constraint from primitive constraints as follows:
\[
\bigwedge_{i=1}^n |v_i| = n ~ \wedge ~ 
    \bigwedge_{i=1}^{n-1} \bigwedge_{j=i+1}^n \exists_{u_{ij}} (u_{ij}
    = v_i \cap v_j \wedge |u_{ij}| \leq 1) 
\]
Unfortunately, the resulting ROBDD turns out to be exponential in
size, making it impractical for use as a global propagator. This is
not surprising in light of NP-hardness of the problem in 
general---in fact it would be surprising if the resulting ROBDD was \emph{not}
exponential in size! 


\subsection{Avoiding Symmetry Through Ordering Constraints}
\label{sect:lexorder}

It is important in modeling a constraint satisfaction problem to 
minimize symmetries in the model of the problem. A model
that contains symmetrical solutions often has a greatly enlarged
search space, leading to large amounts of time being spent
in searching sections of a search tree that are identical up to
symmetric rearrangement. It is
therefore highly desirable to remove whatever symmetries exist in a
problem.

One approach to symmetry-breaking is the introduction of additional
ordering constraints between the variables of the problem. A convenient
ordering to use on sets is a \emph{lexicographic order} on the
characteristic bit vectors of the sets. In other words, if $v$ and $w$
are set variables, then $v < w$ in the lexicographic ordering if and
only if the list of bits $V(v)$ is lexicographically smaller than
$V(w)$. We can model this lexicographic ordering constraint as
$\lexlt(v,w,1)$, which is defined recursively in the following manner:
\[
\lexlt(v,w,n) =
        \begin{cases}
            0 & \text{if } n > N \\
            (\neg v_n \wedge w_n) \vee ((v_n \lra w_n) \wedge
            \lexlt(v,w,n+1) & \text{otherwise}
        \end{cases}
\]
An example of such an ROBDD is depicted in
Figure~\ref{fig:partition-and-lexlt}.

We shall make use of the lexicographic ordering constraint extensively
in our experiments in Section~\ref{sect:exp}.

\section{Modeling Integers, Multisets, and Weighted Sum Constraints}
\label{sect:ints}

In this section we show how to model integer variables and integer
constraints using ROBDDs. In general such representations can be very
large, which limits the usefulness of ROBDDs as the basis for a general
purpose finite-domain constraint solver. Despite this, the ability to
represent integers is extremely useful as a component of a set solver.
We propose two major uses of the integer representation.

%\pjs{reorder uses}
Firstly, we can use the integer representation to model values such as
the weighted sum of the elements of a set.  Constraints on the weighted
sum of elements of a set have been shown to be useful in practical
applications \cite{mailharro:configuration}.

Secondly, we can model finite multisets using ROBDDs by replacing the
individual ROBDD variables of the set representation with bundles of
ROBDD variables, each bundle corresponding to a binary integer. Multiset
operations can then be constructed by composing integer operations on
the variable bundles.

In addition, the integer representation as described here could be used
as an interface between an ROBDD-based set solver and a more
conventional integer finite-domain solver. Such an interface could
easily be implemented using channeling constraints between ROBDD integer
and finite-domain versions of the same variable.

\subsection{Representing Integer Values using ROBDDs}

In order to model integers and integer operations we must choose an
appropriate representation in terms of Boolean formul\ae{}.  In general
we are free to use any encoding of an integer as a binary sequence, such
as unary, unsigned binary and twos-complement encodings, but for
simplicity we choose to represent all of our integers in unsigned binary
form.

\ignore{
(PH) I'm choosing my version of this cvs conflict:

ROBDDs model Boolean formul\ae{}, so in order to model
integers and integer operations we must describe them in terms of
Boolean formulae. First, we must decide how to encode our integers as
sequences of bits. 

In general we are free to use any encoding of an integer as a binary
sequence. In particular, unary and binary encodings suggest themselves
for this application due to their simplicity. Unary is undesirable in
general, since it places
severe limits on the sizes of the integers that we are able to
represent. We concentrate our efforts here on modeling integers in
binary. Fortunately this is well-understood---the 
representation of integers
and integer operations within a computer is binary, so our ROBDDs
can mirror the design of the hardware logic that handles integers. We
will assume that all of our integers are unsigned, but the method
described below could be extended to use mechanisms such as
twos-complement arithmetic with relative ease.
}

We can represent an arbitrary integer expression $e$ by a list of
Boolean formul\ae{}.  Each formula in the list corresponds to a single
bit of the unsigned binary value of the expression. We will denote such
a list by $\n{e_{n-1}, e_{n-2}, \ldots, e_1, e_0}$, where each $e_i$ is
a Boolean formula, in order from the most significant bit to the least
significant bit. We interpret a formula as a ``1'' bit if the formula is
logically true, and as a ``0'' bit otherwise; for simplicity we will
just call the formul\ae{} the ``bits of the expression''. We will also
use the expression $e$ and its list of constituent bits interchangeably.
As usual, we will represent the Boolean formul\ae{} as ROBDDs. As we
shall see, this notation is flexible enough to represent arbitrary
integer expressions.

\begin{example}
  Consider the integer constant $k=25$, which is $11001$ in unsigned
  binary. We will represent $k$ as the list $\n{1,1,0,0,1}$.
\end{example}

In order to represent an integer variable $x$, we associate with $x$ a
fixed set of Boolean variables $\{x_{k-1}, \ldots, x_0\}$. The value of
$x$ is then taken to be the value of the unsigned binary integer
$\n{x_{k-1}, \ldots, x_0}$. By varying the value of the $x_i$ variables,
the value of $x$ can range from $0$ to $2^k - 1$ inclusive.  As always
the ordering of the Boolean variables has a significant effect on the
sizes of the ROBDD representations of formul\ae{}. If $x = \n{x_{k-1},
  \dotsc, x_0}$, $y = \n{y_{k-1}, \dotsc, y_0}$, and $z = \n{z_{k-1},
  \dotsc, z_0}$, then we choose to order the corresponding ROBDD
variables in a interleaved most-significant-bit-first order, i.e.
$x_{k-1} \prec y_{k-1} \prec z_{k-1} \prec x_{k-2} \prec \cdots \prec
x_0 \prec y_0 \prec z_0$.

Since we permit arbitrary Boolean formul\ae{} as the bits of an
expression, we can also model arbitrary integer expressions. For
example, suppose $x$ and $y$ are integer variables with bits $\n{x_2,
  x_1, x_0}$ and $\n{y_2,y_1,y_0}$ respectively. Then, for example, we
can represent the expression $x \wedge y$ (where the $\wedge$ denotes
the bitwise AND operator) as the list $\n{x_2 \wedge y_2, x_1 \wedge
  y_1, x_0 \wedge y_0}$. We are not limited to logical operations, as we
shall see in the next section.

\ignore{
\begin{example}
  Suppose $x$ is an integer variable, where $x=3$ or $x=5$. $3_{10} =
  110_2$, so $x=3$ can be represented by the formula $x_1 \wedge x_2
  \wedge \neg x_3$. Similarly $5_{10} = 101_2$, so $x=5$ can be
  represented by the formula $x_1 \wedge \neg x_2 \wedge x_3$. The
  constraint that $x=3$ or $x=5$ can be represented by the disjunction
  $(x_1 \wedge x_2 \wedge \neg x_3) \vee (x_1 \wedge \neg x_2 \wedge
  x_3)$, which simplifies to just $x_1 \wedge (x_2 \lra \neg x_3)$. We
  can represent this formula as an ROBDD in the usual way.
\end{example}
}

\subsection{Representing Integer Operations using ROBDDs}
\label{sect:int-ops}

We can also use ROBDDs to model arithmetic operations such as addition,
by analogy with the design of the corresponding logic circuits. For the
purpose of the set and multiset solvers, we only require implementations
of the operations of addition, left shift, minimum, maximum,
multiplication by a constant, and multiplication by a single variable
bit. We do not require a general implementation of multiplication using
ROBDDs.

It is convenient to assume that all integer expressions have the same
number of bits. We may assume this without loss of generality since we
can freely pad the left of the shorter of any pair of expressions with
``0'' bits.

To model addition, we simulate the
operation of a full binary adder. Suppose $x$ and $y$ are integer
expressions, with bit representations $\n{x_{l-1}, \ldots, x_0}$ and
$\n{y_{l-1}, \ldots, y_0}$. 
We can use ROBDDs to compute
the output bits $\plus(x,y)$ of the operation $x+y$
as follows (here
$c_i$ denotes a carry bit and $s_i$ denotes a sum bit):
\begin{align*}
    c_{-1} &= 0 &\\
    s_i &= x_i \oplus y_i \oplus c_{i-1} & 0 \leq i <
    l\\
    c_i &= (\neg c_{i-1} \wedge x_i \wedge y_i) \vee
          (c_{i-1} \wedge (x_i \vee y_i)) &0 \leq i < l\\
    \plus(x,y) &= \langle c_{l-1}, s_{l-1}, \ldots, s_1, s_0 \rangle &
\end{align*}
Note we avoid overflow by extending the size of the result by one bit.
\begin{example}
  
  Suppose $x$ is an integer variable, with bits $\n{x_1, x_0}$. We can
  represent the expression $x + 3$ by the bits $\n{ (\neg x_0 \wedge
    x_1) \vee x_0), x_1 \oplus 1 \oplus x_0, x_0 \oplus 1 } = \n{x_0
    \vee x_1, x_0 \lra x_1, \neg x_0}$.

\end{example}

The left shift operation is trivial to implement. If $x$ is an integer
expression represented by $\n{x_{l-1}, \ldots, x_0}$ and $k$ is a
non-negative integer, then we can represent the left shift of $x$ by $k$
bits $\shl(x,k)$ by the following formula:
\begin{equation*}
    \shl(\n{x_{l-1}, \ldots, x_0},k) = \langle
    x_{l-1}, \dotsc, x_0,
    \underbrace{0, \dotsc,
    0}_{\text{k bits}} \rangle
\end{equation*}


We can implement the operation of multiplication by a constant using the
$\plus$ and $\shl$ operators. If $x$ is an integer expression with bits
$\n{x_{l-1}, \dotsc, x_0}$ and $k$ is a non-negative integer, then $x
\times k$ corresponds to $\mul(x,k)$ in the following formula:
\begin{align}
  \label{eqn:mul-expr}
    \mul(x,k) &= 
        \begin{cases}
            \n{} & \text{ if } k = 0\\
            \mul(\shl(x, 1), \frac{k}{2}) & \text{ if $k$ is even and
            $k > 0$} \\
            \plus(x, \mul(\shl(x, 1), \lfloor \frac{k}{2} \rfloor)) &
                \text{ if $k$ is odd and $k > 0$} 
        \end{cases}
\end{align}

\begin{example}
  Let $x = \n{x_1,x_0}$, and consider the expression $e = \mul(x, 3)$.
  By applying Equation~\eqref{eqn:mul-expr}, we obtain $e = \n{x_1
    \wedge (x_1 \wedge x_0), x_1 \oplus (x_0 \wedge x_1), x_0 \oplus
    x_1, x_0}$ which can be simplified to $\n{ x_0 \wedge x_1, \neg x_0
    \wedge x_1, x_0 \oplus x_1, x_0 }$. The corresponding ROBDD
  representations are shown in Figure~\ref{fig:mul-exprs}.
\end{example}

\begin{figure}
\[
  \begin{array}{ccccccc}
    e_3 && e_2 && e_1 && e_0 
    \\
    \\
    \xymatrix@=10pt{
      &\xyo{x_1}\ar@/_5pt/@{-->}[ddl]\ar[d]\\
      &\xyo{x_0}\ar@/_5pt/@{-->}[dl]\ar[d]\\
      \xyfalse & \xytrue
    }
    &~~~~&
    \xymatrix@=10pt{
      \xyo{x_1}\ar@/^5pt/[dr]\ar@/_5pt/@{-->}[dd]\\
      &\xyo{x_0}\ar@/^15pt/@{-->}[d]\ar[dl]\\
      \xyfalse & \xytrue
    }
    &~~~~&
    \xymatrix@=10pt{
      &\xyo{x_1}\ar@/_5pt/@{-->}[dl]\ar@/^5pt/[dr]\\
      \xyo{x_0}\ar@/_15pt/@{-->}[d]\ar[drr] && \xyo{x_0}\ar@/^15pt/@{-->}[d]\ar[dll]\\
      \xyfalse && \xytrue
    }
    &~~~~&
    \xymatrix@=10pt{
      &\xyo{x_0}\ar@/_5pt/@{-->}[dl]\ar@/^5pt/[dr]\\
      \xyfalse && \xytrue
    }
    \\
    \\
    (x_0 \wedge x_1) &&
    (\neg x_0 \wedge x_1) &&
    (x_0 \oplus x_1) &&
    (x_0)
    \\
  \end{array}
\]
\caption{ROBDDs representing the bits $\n{e_3, e_2, e_1, e_0}$ of the
expression $e = \mul(\n{x_1, x_0},3)$, together with
the simplified Boolean expressions for each $e_i$.}
\label{fig:mul-exprs}
\rule{\textwidth}{0.4pt}\\
\end{figure}

%%We will also require the operation of multiplication of an expression
%%$x$ by a single bit $b$, which corresponds to $\bmul(x,b)$ in the
%%following equation:
%%\[
%%    \bmul(\n{x_{l-1},\dotsc,x_1,x_0}, b) = \n{x_{l-1} \wedge b, \dotsc,
%%    x_1 \wedge b, x_0 \wedge b}
%%\]

\subsection{Integer Constraints using ROBDDs}
\label{sect:int-cons}

We can also express constraints on integer expressions using ROBDDs. In
particular, we show how to implement equality and inequality
constraints.  As usual, we assume any two expressions $x$ and $y$ have
equal lengths of $l$ bits; if not, we pad the shorter expression with
``0'' bits on the right.

Equality of two integer expressions $x$ and $y$ is easy to represent as
an ROBDD---the corresponding bits of $x$ and $y$ must be equal.  Hence
we can represent the equality constraint $x = y$ as the ROBDD $B(x =
y)$:
\[
B(x = y) = \bigwedge_{i=0}^{l-1} x_i \lra y_i
\]
Note that this is identical to the implementation of equality for two
set expressions with the addition of zero-padding.

It turns out that we already have an implementation for inequality
constraints, albeit in a different guise. Inequalities of the binary
integers correspond to inequalities on the lexicographic ordering of the
bit representations, so we can implement, for example, the strict
less-than constraint $x < y$ for two integer variables $x$ and $y$ using
the $\lexlt$ operation from Section~\ref{sect:lexorder}.
\[
B(x < y) = \lexlt(x, y)
\]
We can then use this to construct the reverse inequality by swapping
the order of the operands, and the non-strict inequality by negating
the formula (reversed as necessary).

The implementation of inequalities also leads us to an implementation of
the minimum and maximum expressions. Consider the problem of finding the
smaller of two integer expressions $x$ and $y$. If $x$ and $y$ have bit
vectors $\n{x_{l-1}, \dotsc, x_0}$ and $\n{y_{l-1}, \dotsc, y_0}$
respectively, we recursively define $\min(x,y)$ as follows:
\begin{align*}
L_{l} & = R_{l} = 0 \\
L_i & = L_{i+1} \vee (\neg L_{i+1} \wedge \neg R_{i+1} \wedge \neg x_i \wedge y_i) & 1
\leq i < l\\
R_i & = R_{i+1} \vee (\neg L_{i+1} \wedge \neg R_{i+1} \wedge x_i \wedge \neg y_i) &
1 \leq i < l \\
m_i & = (L_{i+1} \wedge x_i) \vee (R_{i+1} \wedge
y_i) \vee (\neg L_{i+1} \wedge \neg R_{i+1} \wedge x_i \wedge y_i) &
0 \leq i < l\\
\min(x,y) &= \n{m_{l-1}, \dotsc, m_0}
\end{align*}
In the equation above, the $L_i$ and $R_i$ values are flag bits which
state whether the higher order bits have already allowed us to
conclude that one of the two
values is the minimum. The maximum operation is defined similarly.

\subsection{Modeling Multisets and Multiset Constraints}
\label{sect:multisets}

Various authors have suggested that multisets are a valuable addition
to the modeling abilities of a set constraint solver
\cite{kiziltanw:multisets}. 
\ignore{
\pjs{omit inflamatory statement}
This is a questionable assertion---constraints 
on multisets are equivalent to global constraints on
integers, so the expressive power of such a solver would not seem to
be greatly enlarged. Nonetheless,}
In this section we briefly show how multisets and
multiset constraints can be modelled using ROBDDs by making use of the
integer building blocks described above.

A multiset $m$ is a unordered list of elements $\llbrace m_0, \dotsc,
m_n \rrbrace$ drawn from the universe $\mathcal{U}$, in which (unlike a set)
repetition of elements is
permitted. Most set operations have parallel operations on
multisets, although the multiset operations are not strict
generalisations of the set operations. Let $\occ(i,m)$ denote the
number of occurrences of an element $i$ in a multiset $m$. Suppose $m$
and $n$ are multisets, and $k$ is an integer constant. We define
the following multiset relations and operations by their actions on the number
of occurrences of each element $i$ in the universe:
\begin{itemize}
\item Equality: $m=n$ iff $\occ(i,m) = \occ(i,n)$
for all $i \in \mathcal{U}$.
\item Subset:  $m \subseteq n$ iff $\occ(i,m) \leq \occ(i,n)$ for
all $i \in \mathcal{U}$.
\item Union: $\occ(i,m\cup n) = \occ(i,m) + \occ(i,n)$ for all $i
\in \mathcal{U}$.
\item Intersection: $\occ(i,m\cap n) = \min\{\occ(i,m),
\occ(i,n)\}$ for all $i \in \mathcal{U}$.
\item Difference: $\occ(i,m \setminus n) = \max\{0, \occ(i,m) -
\occ(i,n)\}$ for all $i \in \mathcal{U}$.
\item Cardinality: $|m| = \sum_{i\in\mathcal{U}} \occ(i,m)$ for all
$i \in \mathcal{U}$.
\end{itemize}
%\pjs{replaced strict subset by subset, since def was wrong!}

To represent a set variable $x$ we associated a vector of Boolean
variables $\n{x_1,\dotsc,x_n}$ with the bits of the characteristic
vector of a valuation of $x$. In the case of a multiset, the characteristic
vector of a multiset $m$ is a vector of
integers, and so we need to associate an integer value $m_i$ with each
potential element $i$ of the multiset $m$. We can model each such
integer value using the approach described above.

If $m$ is a multiset variable, then we associate a \emph{bundle} of
ROBDD variables with every $i \in \mathcal{U}$, the contents of which
comprise the
bits of the corresponding integer expression $m_i$.
In order that we can represent multisets in a finite
number of bits, we assume that the number of occurrences of any element
in a multiset variable is bounded above by a reasonably small $M$,
allowing us to use only
$k = \lceil \log_2 M \rceil$ Boolean variables per bundle. We then 
write $m$ as a list of bundles $\n{m_1, \dotsc, m_n}$, where each
$m_i$ is in turn a list of bits $\n{m_{i,k-1}, \dotsc, m_{i,0}}$.

Given a representation of multiset variables, we now turn our
attention to the implementation of multiset expressions and
constraints.
Multiset expressions
can be implemented in the obvious way as sequences of integer
expressions. For example, suppose $x$ and $y$ are multiset variables
with associated bundles $\n{x_1,\dotsc,x_N}$ and $\n{y_1,\dotsc,y_N}$
respectively. Then the bundles corresponding to the expression $x \cup
y$ are $\n{\plus(x_1, y_1), \plus(x_2, y_2), \dotsc, \plus(x_N,
y_N)}$. Similarly, the bundles corresponding to the expression $x \cap
y$ are $\n{\min(x_1,y_1), \dotsc, \min(x_N,y_N)}$, and so on for 
other expressions.  We show how to implement cardinality and
weighted sum constraints in Section~\ref{sect:wsum}.

Multiset constraints are also trivial to implement---for example, 
two multisets $x$
and $y$ are equal if and only if all of their constituent bundles are
equal, and so multiset equality can be modelled by a conjunction of integer
equalities.
Relations
such as subset correspond to a conjunction of integer inequalities on the
constituent
bundles, the implementation of which was described in
Section~\ref{sect:int-cons}.

Up until this point we have left the ordering of the ROBDD variables that
make up a multiset variable unspecified. Unfortunately, unlike the set
case, there is no single optimal variable ordering that is guaranteed
to produce compact descriptions of all the primitive multiset constraints.
For example, a subset constraint can be compactly represented by a bundle-major
bit ordering but not a bit-major ordering (see
Figure~\ref{fig:multiset-bit-orders}), since a subset constraint
consists of a series of integer inequalities between the corresponding
bundles of two multiset variables, and so a bundle-major ordering gives an
interleaving of variables as above.
However, the opposite is true for a cardinality constraint, which
consists of a sum of the values of the bundles within a variable,
the variables of which are interleaved under a bit-major ordering but not a
bundle-major ordering. These two orderings are mutually
exclusive, and hence we can conclude that in general there need not be
an optimal variable ordering for modeling multiset constraints.
\begin{figure}[t]
\begin{center}
\begin{tabular}{c@{~~~~~}c}
\vspace{-\ht\strutbox}
%\vspace{0pt}
\includegraphics[height=2.5cm]{bitmajor}
&
\vspace{-\ht\strutbox}
\includegraphics[height=2.3cm]{bundlemajor} \\
\\
\\
Bit major & Bundle major
\end{tabular}
\end{center}
\caption{Bit and bundle major orderings for two multiset variables}
\label{fig:multiset-bit-orders}
\rule{\textwidth}{0.4pt}
\end{figure}

\subsection{Weighted Sum and Cardinality Constraints}
\label{sect:wsum}

In many practical applications we are interested in placing
constraints on a weighted sum of the elements of a set variable. For
example, in the Balanced Academic Curriculum
problem (problem \texttt{prob030} of CSPLib), every
course has an associated weight corresponding to its academic load, and
there is a limit to the total academic load that can be undertaken during
any given time period.
In the set model of the problem proposed by
\citeA{hnichkw:bacp}, the limits on the academic load in a period are made
through the use of a weighted sum constraint.
In addition the cardinality
constraint which we have already described is a special case of the
weighted sum constraint where all the weights are set to ``1''. It
therefore seems essential to implement this constraint in the ROBDD
framework.

Suppose $x$ is a set or multiset expression with bit bundles
$\n{x_1, \dotsc, x_n}$. (If $x$ is a set expression, then each bit
bundle has size $1$). We can use the integer operations described
earlier to produce an integer expression 
$\wsum(x, w)$ corresponding to the weighted sum $\sum_{i=1}^n x_i w_i$,
where $w$ is a vector of integers $\n{w_1, \dotsc, w_n}$:
\begin{align*}
\phi_0 &= \n{} \\
\phi_i &= \plus(\phi_{i-1}, \mul(x_i, w_i)) \\
\wsum(x,w) &= \phi_n 
\end{align*}


Expressions involving the cardinality of set and multiset variables can then be
expressed as a special case of the weighted sum expression where $w_i
= 1$ for all $1 \leq i \leq n$. We have already seen one method of
constructing ROBDDs for the constraints $|x| = k$, $|x| \leq k$ and
$|x| \geq k$ in the set variable case; this method is more general,
since it permits us to directly model constraints such as $|x| + 5 \leq |y| +
|z|$.
This is of great practical value---for example it is needed to
implement the Hamming Code experiment of
Section~\ref{sect:exp-hamming}.  In the cases already discussed in
Section~\ref{sect:modeling-constraints} the ROBDDs produced by the
two methods are identical since ROBDDs are a canonical 
representation.



\section{Efficient Constraint Propagation Using ROBDDs}
\label{sect:propagation}

In this section we discuss improvements to and variants on the basic
domain propagation scheme presented in
Section~\ref{sect:basic-set-solver}. 

The implementation of the domain propagator $\dom(c)$ has a
substantial effect on the performance of the solver. The definition of
$\dom(c)$ was given in its purest mathematical form for simplicity. 
Section~\ref{sect:prop-domain} discusses several implementation details that
can lead to greatly improved efficiency when performing domain propagation.

In general, the inferences that can be obtained using the domain 
propagator may be too costly to be practical, and so in some
circumstances it may be desirable to enforce a weaker form of
consistency. Weaker consistency may be substantially cheaper to
enforce, permitting more time to be spent in searching the solution
space.

As always there is a compromise between propagation time 
and search time, and for certain problems it may be more productive to
spend more time searching than performing more accurate propagation, while
for other problems the converse may be true. The
ROBDD based representation allows this to be taken to extremes---in
theory it is possible to form a single ROBDD representing the
solutions to a constraint satisfaction problem by forming the
ROBDD conjunction of the constraints. A solution could then be
trivially read from the ROBDD as a satisfying assignment. 
Usually such an ROBDD would be prohibitively expensive
to construct in both time and space, forcing us to maintain less
strict consistency levels.

Accordingly, we show how to implement weaker levels of consistency
using the ROBDD representation by combining the domain propagator
with an approximation operation. The approximation operation simplifies
the ROBDD representing a domain to its bounds closure under a suitable
definition of bounds. The ROBDD representing the bounds of a
domain is almost always smaller than the domain itself, leading to
better performance for all future
operations on that domain.
As  we shall see, this can lead to substantial performance
improvements for the overall solver.

\subsection{Domain Propagation}
\label{sect:prop-domain}

Let $c$ be a constraint, with $\vars(c) = \{v_1, \dotsc, v_n\}$. 
Section~\ref{sect:basic-set-solver} gave the following definition of
a domain propagator:
\begin{equation}
\label{eqn:domain-propagator}
    \dom(c)(D)(v_i) = \bexists_{V(v_i)} (B(c) \wedge \bigwedge_{j=1}^n
    D(v_j))
\end{equation}


Since $B(c)$ and $D(v_j)$ are ROBDDs, we can directly implement
Equation~\eqref{eqn:domain-propagator} using ROBDD operations. In
practice it is more efficient to perform the existential
quantification as early as possible to limit the size of the
intermediate ROBDDs. Here we make use of
the efficient combined
conjunction and existential quantification operation,
which we'll call \emph{and-abstraction},
provided by most ROBDD packages.

This leads to the following implementation:
\begin{equation}
\label{eqn:domain-impl}
\begin{split}
  \phi_i^0 &= B(c) \\
  \phi_i^j &=
  \begin{cases}
    \exists_{V(v_j)} (D(v_j) \wedge \phi_i^{j-1}) & 1 \leq i, j \leq n,
    \ i \neq j ~\text{(and-abstraction)}\\
    \phi_i^{i-1} & i = j \\
  \end{cases} \\
  \dom(c)(D)(v_i) & = D(v_i) \wedge \phi_i^n\\
\end{split}
\end{equation}
The worst case complexity is still $O(|B(c)| \times \prod_{j=1}^{n} |D(v_j)|)$
for each $v_j$. 
Clearly some of the computation can be shared between propagation of
$c$ for different variables since $\phi_i^j = \phi_{i'}^j$ when $j < i$
and $j < i'$. Even with this improvement the algorithm of
Equation~\eqref{eqn:domain-impl} uses $O(n^2)$ and-abstraction
operations (which experimentally have been shown to occupy the
majority of the execution time of the set solver).

The domain propagator implementation of Equation~\eqref{eqn:domain-impl}
can be significantly improved by observing that in the case of an
$n$-variable constraint ($n \geq 3$) many similar sub-formul\ae{} are computed.
Due to the need to perform the existential quantification operations
as early as possible, we do not have complete freedom to rearrange the
order of evaluation as we see fit. However, a simple
divide-and-conquer strategy for calculating $\dom(c)(D)(v_i)$ allows us
to perform domain propagation using just $O(n \log n )$ and-abstraction
operations. We define $dom(c)(D) = \mathsf{dom\_divide\_conquer}(D,
B(c), \vars(c))$, where \textsf{dom\_divide\_conquer} is defined 
in Figure~\ref{fig:dom-divide-conquer}.

\begin{figure}[t]
\begin{quote}
\begin{tabbing}
xx \= xx \= xx \= xx \= xx \= \kill
\textsf{dom\_divide\_conquer}($D$, $\phi$, $V$): \\
\> \textbf{if} ($|V| = 0$) \textbf{return} $D$ \\
\> \textbf{else if} ($V = \{v_i\})$ \\
\> \> $D(v_i)$ := $D(v_i) \wedge \phi$ \\
\> \> \textbf{return} $D$ \\
\> \textbf{else} \\
\> \> $\{v_1, \dotsc, v_k\}$ := $V$ \\
\> \> $h$ := $\lfloor \frac{k}{2} \rfloor$ \\
\> \> $R$ := $\exists_{V(v_1)} D(v_1) \wedge \exists_{V(v_2)} D(v_2) \wedge \cdots
\exists_{V(v_h)} D(v_h) \wedge \phi$ \\
\> \> $L$ := $\exists_{V(v_{h+1})} D(v_{h+1}) \wedge \exists_{V(v_{h+2})}
D(v_{h+2}) \wedge \cdots \exists_{V(v_k)} D(v_k) \wedge \phi$ 
 \\
\> \> $D_1$ := \textsf{dom\_divide\_conquer}($D$, $L$, $\{v_1,\dotsc,v_h\}$) \\
\> \> $D_2$ := \textsf{dom\_divide\_conquer}($D_1$, $R$, $\{v_{h+1},\dotsc,v_k\}$) \\
\> \> \textbf{return} $D_2$
\end{tabbing}
\end{quote}
\caption{A divide and conquer algorithm for domain propagation}
\label{fig:dom-divide-conquer}
  \rule{\textwidth}{0.4pt}
\end{figure}

\ignore{
\begin{example}
Consider the evaluation of a constraint $c$ where $var(c) =
\{v_1,v_2,v_3,v_4\}$.  If we are to simultaneously compute
$\dom(c)(D)(v_i)$ for $1 \leq i \leq 4$ then using the
implementation of Equation~\ref{eqn:domain-impl} we compute
the following subformul\ae{}
$$
\begin{array}{l|llll}
\phi_i^j    &j = 1& j = 2 & j =3 & j = 4 \\ 
\hline
i=1 & B(c) & \exists_{V(v_2)} D(v_2) \wedge \phi_1^1
&\exists_{V(v_3)} D(v_3) \wedge \phi_1^2 &
\exists_{V(v_4)} D(v_4) \wedge \phi_1^3  \\
i = 2 & \exists_{V(v_1)} D(v_1) \wedge B(c) &
 \phi_2^1 
&\exists_{V(v_3)} D(v_3) \wedge \phi_2^2 &
\exists_{V(v_4)} D(v_4) \wedge \phi_2^3  \\
i = 3 & \phi_2^1 & 
\exists_{V(v_2)} D(v_2) \wedge \phi_3^1 &
\phi_3^2 &
\exists_{V(v_4)} D(v_4) \wedge \phi_3^3  \\
i = 4 & \phi_2^1 & 
\phi_3^2 &
\exists_{V(v_3)} D(v_3) \wedge \phi_4^2 &
\phi_4^3  \\
\end{array}
$$
Overall there are 9 and-abstraction operations, saving 3 out of
the potential 12.

Using \textsf{dom\_divide\_conquer} we compute
\textsf{dom\_divide\_conquer}($D$, $B(c)$, $\{v_1,v_2,v_3,v_4\}$)
which computes
$R_{1,2} = \exists_{V(v_1)} D(v_1) \wedge \exists_{V(v_2)} D(v_2) \wedge
B(c)$
and
$L_{3,4} =  \exists_{V(v_3)} D(v_3) \wedge \exists_{V(v_4)} D(v_4) \wedge
B(c)$
calls
\textsf{dom\_divide\_conquer}($D$, $L_{3,4}$, $\{v_1,v_2\}$)
and
\textsf{dom\_divide\_conquer}($D$, $R_{1,2}$, $\{v_3,v_4\}$).
The first call computes
$R_{1,3,4} = \exists_{V(v_1)} D(v_1) \wedge L_{3,4}$
and
$L_{2,3,4} = \exists_{V(v_2)} D(v_2) \wedge L_{3,4}$
and calls 
\textsf{dom\_divide\_conquer}($D$, $L_{2,3,4}$, $\{v_1\}$)
and
\textsf{dom\_divide\_conquer}($D$, $R_{1,3,4}$, $\{v_2\}$)
which both immediately update the domains of $v_1$ and
$v_2$. The other call  
\textsf{dom\_divide\_conquer}($D$, $R_{1,2}$, $\{v_3,v_4\}$).
computes
$R_{1,2,3} = \exists_{V(v_3)} D(v_3) \wedge R_{1,2}$
and
$L_{1,2,4} = \exists_{V(v_4)} D(v_4) \wedge R_{1,2}$
and calls 
\textsf{dom\_divide\_conquer}($D$, $L_{1,2,4}$, $\{v_3\}$)
and
\textsf{dom\_divide\_conquer}($D$, $R_{1,2,3}$, $\{v_4\}$)
which immediately update domains.
Overall there are 8 and-abstractions.  In this case the improvement 
over Equation~\ref{eqn:domain-impl} may
seem minor, but as $n$ grows the advantage is significant. 
\end{example}
}


%\pjs{probably omit above example, boring!}
%\pjs{Add explanation and example of divide and conquer in action!}

\subsection{Set Bounds Propagation}
\label{sect:set-bounds}

As domain propagation may be prohibitively expensive to enforce for some
problems, it is useful to investigate less strict notions of
consistency. In such cases, we can speed propagation by simplifying the
domains through approximation. Since set bounds under the subset partial
ordering relation are one of the most commonly used approximations to
a set domain, it seems natural to implement a set bounds propagator in
the ROBDD framework. Only relatively minor changes are needed to the
domain propagator to turn it into a set bounds propagator.

Given the ROBDD representation of a set domain, we can easily identify
the corresponding set bounds. In an ROBDD-based set domain representation,
the set bounds on a domain correspond to the \emph{fixed} variables
of the ROBDD representing the domain. We say an ROBDD variable $v$ is
\emph{fixed} if either for all nodes $n(v,t,e)$ $t$ is the constant
$0$ node, or for all nodes $n(v,t,e)$ $e$ is the constant $0$ node,
and such a node appears in every path from the root of the diagram to
the ``1'' node.

Such nodes can be identified in a single pass over the domain ROBDD,
in time proportional to its size. If $\phi$ is an ROBDD, we will write
$\essential{\phi}$ to denote the ROBDD representing the conjunction of
the fixed variables of $\phi$. If $\phi$ represents a set of sets $S$,
then $\essential{\phi}$ represents $\conv(S)$. 
An ROBDD $\phi$ where $\phi = \essential{\phi}$ is a stick ROBDD by definition.

\begin{example}
    Let $\phi$ be the ROBDD depicted in Figure~\ref{fig:g}(c). Then
    $\essential{\phi}$ is the ROBDD of Figure~\ref{fig:g}(a).
\end{example}

Using this operation we can convert our domain propagator 
into a set bounds propagator by discarding all of the non-fixed
variables from the domain ROBDDs after each propagation step. Suppose
that $D(v)$ is a stick ROBDD for each $v \in \mathcal{V}$. If $c$ is a
constraint, with $\vars(c) = \{v_1, \dotsc, v_n \}$, we define a set
bounds propagator $\setb(c)$ thus:
\begin{equation}
\label{eqn:bounds-propagator}
 \setb(c)(D)(v_i) = \Big\llbracket \bexists_V(v_i) (B(c) \wedge
 \bigwedge_{j=1}^n D(v_j)) \Big\rrbracket
\end{equation}

Despite only relatively minor differences between the set bounds
propagator and the domain propagator, the set bounds propagator is
usually significantly faster than a domain propagator for two reasons.
Firstly, as the domains $D(v)$ are all sticks, all of the ROBDD
operations are cheap, compared to operations on the possibly very large
ROBDDs representing arbitrary domains. The entire propagator can be
implemented with $O(|B(c)|)$ complexity, since all of the other ROBDDs
are sticks.  Secondly, we can use the updated set bounds to simplify the
propagator ROBDD $B(c)$. Since domains are monotonic decreasing in size,
fixed variables will remain fixed up to backtracking, and so we can
project them out of $B(c)$, thus reducing the size of the propagator
ROBDD in future propagation steps. This leads us to the following
implementation of the propagator:
%\pjs{Confused, isnt the correct implementation just}
\begin{equation}
\label{eqn:bounds-propagator-impl}
\begin{array}{r@{\,}lr}
\phi_0 &= B(c) & \\
\phi_j &= \exists_{VAR(D(v_j))} D(v_j) \wedge  \phi_{j-1} & 1 \leq j
\leq n \\
sb(c)(D)(v_i) &= D(v_i) \wedge 
                \bar{\exists}_{V(v_i)} \llbracket \phi_n \rrbracket &
                1 \leq i \leq n
\end{array}
\end{equation}
After this propagation step we can replace the representation
of the constraint $B(c)$ by $\phi_n$ since the fixed variables
will no longer have any new impact.

\begin{example}\label{ex:bp}
Consider bounds propagation for the constraint $c \equiv v \subseteq w$
where $N = 3$. The ROBDD representation $B(c)$ is given
in Figure~\ref{fig:subset-bdd}. Assume the domains of $v$ and $w$
are respectively $[\{1\},\{1,2,3\}]$, represented by the formula $v_1$,
and  $[\emptyset,\{1,2\}]$, represented by the formula $\neg w_3$.
The ROBDD 
$\phi_n \equiv \exists v_1 \exists w_3 ~B(c) \wedge v_1 \wedge \neg w_3$ 
is shown in
Figure~\ref{fig:bp}(a).  We have that
$\llbracket \phi_n \rrbracket = w_1 \wedge \neg v_3$.
We can project out the
fixed variables $v_1, w_1, v_3, w_3$ from $B(c)$ to get a new simplified
form of the constraint $v_2 \rightarrow w_2$ shown in Figure~\ref{fig:bp}(b).
\end{example}
\begin{figure}
$$
\begin{array}{ccc}
    \xymatrix@=10pt{
%     \xyo{v_1}\ar@{->}[dr]\ar@/_30pt/@{-->}[dddddrrrrr]\\
    &\xyo{w_1}\ar@{->}[dr]\ar@/_25pt/@{-->}[ddddrrrr]\\
    &&\xyo{v_2}\ar@{->}[dr]\ar@/_20pt/@{-->}[ddrr]\\
    &&&\xyo{w_2}\ar@{->}[dr]\ar@/_15pt/@{-->}[ddrr]\\
    &&&&\xyo{v_3}\ar@{->}[dr]\ar@{-->}[dll]\\
    &&*+[F]{\mathrm{\True}} &&& *+[F]{\mathrm{\False}}
    }
&
~~
&
    \xymatrix@=10pt{
    &&\xyo{v_2}\ar@{->}[dr]\ar@/_15pt/@{-->}[ddrr]\\
    &&&\xyo{w_2}\ar@{->}[dr]\ar@{-->}[dl]\\
    &&*+[F]{\mathrm{\False}} && *+[F]{\mathrm{\True}}
    }
\\
(a) & & (b)
  \end{array}
$$
\caption{Set bounds propagation on the constraint $v \subseteq w$
showing (a) resulting ROBDD after conjunction with domains,
and (b) simplified constraint ROBDD after removing fixed variables.}
\label{fig:bp}
\rule{\textwidth}{0.4pt}
\end{figure}  

This set bounds solver retains all of the modeling advantages of the
domain solver, including the ability to easily conjoin and
existentially quantify constraints, to remove intermediate variables
and to form global constraints. In some cases this permits a 
substantial performance
improvement over more traditional set bounds solvers.

Experimentally it appears that a direct implementation of
Equation~\eqref{eqn:bounds-propagator}, written to use a divide and conquer
approach to calculate the and-abstractions, is faster than an
implementation of Equation~\eqref{eqn:bounds-propagator-impl}, even if
a divide and conquer approach is used in calculating the existential
quantification. The former approach calculates fewer intermediate
results, which leads to a faster propagator overall.
%\pjs{What exactly should be implemeneted using divide and conquer}
Experimental results for the bounds propagator are given in
Section~\ref{sect:exp}.
\ignore{

the
following implementation:
\ph{fix this}
\begin{quote}
\begin{tabbing}
xx \= xx \= xx \= xx \= xx \= \kill
\textsf{setb\_divide\_conquer}($D$, $\phi$, $V$): \\
\> \textbf{if} ($|V| = 0$) \\
\> \> \textbf{return} $D$ \\
\> \textbf{else if} ($V = \{v_i\})$ \\
\> \> $D(v_i)$ := $D(v_i) \wedge \essential{\phi}$ \\
\> \> \textbf{return} $D$ \\
\> \textbf{else} \\
\> \> $\{v_1, \dotsc, v_k\}$ := $V$ \\
\> \> $h$ := $\lfloor \frac{k}{2} \rfloor$ \\
\> \> $R$ := $\exists_{V(v_1)} v_1 \wedge \exists_{V(v_2)} v_2 \wedge \cdots
\exists_{V(v_h)} v_h \wedge \phi$ \\
\> \> $L$ := $\exists_{V(v_{h+1})} v_{h+1} \wedge \exists_{V(v_{h+2})}
v_{h+2} \wedge \cdots \exists_{V(v_k)} v_k \wedge \phi$ \\
\> \> $D$ := \textsf{setb\_divide\_conquer}($D$, $L$, $\{v_1,\dotsc,v_h\}$) \\
\> \> $D$ := \textsf{setb\_divide\_conquer}($D$, $R$, $\{v_{h+1},\dotsc,v_k\}$) \\
\> \> \textbf{return} $D$
\end{tabbing}
\end{quote}

}


\subsection{Split Domain Propagation}
\label{sect:split}

We can combine the set bounds propagator with the domain propagator to
produce a more space efficient split domain propagator. By separating the
domain representation into fixed and unfixed parts, we can reduce the
total size of the representation, also hopefully speeding propagation.

One of the unfortunate characteristics of ROBDDs is that the size of
the ROBDD representing a domain is highly dependent on the variable
ordering.  Consider an ROBDD representing a set domain which
contains several fixed variables. If these variables do not
appear at the beginning of the variable ordering, then the ROBDD
representing the domain will in effect contain several copies of the
sticks representing the fixed variables. For example,
Figure~\ref{fig:g}(c) contains several copies of the stick in
Figure~\ref{fig:g}(a). Since many of our ROBDD operations take time
proportional to the product of the number of ROBDD nodes of their
arguments, this overly large representation has a performance cost. We
can solve this problem in two ways---either by reordering the
ROBDD variables or by splitting up the domain representation.

Variable reordering is capable of eliminating redundancy in the
representation of any individual domain, but in general cannot
eliminate redundancy across a set of domains.  By reordering the ROBDD
variables, we can reduce the size of a domain by placing
the fixed variables at the beginning of the variable order, thus
removing the unnecessary duplication for that domain. Unfortunately, the
variable order is a global property of all ROBDDs in existence,
whereas the fixed variables of a domain are a local property specific
to a particular domain, so there may not be a variable
ordering that is optimal for all of the domains in a problem.

In the context of applying ROBDDs to the groundness analysis of logic
programs, \citeA{bagnara:96reactive} demonstrated that the performance
of an ROBDD-based program analyzer could be improved by splitting up
ROBDDs into their fixed and non-fixed parts. We can apply the same
technique here.

We split the ROBDD representing a domain $D(v)$ into a pair of ROBDDs
$(\lu, R)$.  $\lu$ is a stick ROBDD representing the lower and upper set
bounds on $D(v)$, and $R$ is a remainder ROBDD representing the
information on the unfixed part of the domain. Logically $D = \lu \wedge
R$. We will write $\lu(D(v))$ and $R(D(v))$ to denote the $\lu$ and $R$
parts of $D(v)$ respectively.

The following results provide an upper bound of the size
of the split domain representation:
\begin{lemma}
Let $G$ be an ROBDD, and let $v$ be a fixed variable of $G$. Then
$|\exists_v G| < |G|$.
\end{lemma}
\begin{proof}
    Since $v$ is a fixed variable, either for every node $n(v,t,f)$ in $G$
    $t$ is the constant 0 node, or for every such node 
    $f$ is the constant 0 node. Since a node $n(v,t,f)$ corresponds to
    the proposition $(v \wedge t) \vee (\neg v \wedge f)$, it is clear
    that $\exists_v n(v,t,f)$ corresponds to simply $t \vee f$, and
    moreover since $v$ is a fixed node one of $t$ or $f$ is zero.
    Hence the existential quantification of a fixed variable $v$
    simply removes all nodes labelled $v$ from $D$. Since there is at
    least one such node, the result follows.
\end{proof}

\begin{lemma}
  \label{prop:split-sizes}
  Let $G$ be an ROBDD, $\lu = \essential{G}$, and $R = \exists_{\VAR(\lu)}
  G$. Then $G \lra \lu \wedge R$ and $|\lu| + |R| \leq
  |G|$.
\end{lemma}
\begin{proof}
    The result $G \lra \lu \wedge R$ is straightforward, so we
    only prove the result on sizes.
    
    Suppose that $\VAR(\lu) = \{v_1, \dotsc, v_n\}$
    is the set of fixed variables of $G$. Then, since
    $|\exists_{v_1} G| < |G|$, $1 + |\exists_{v_1} G| \leq |G|$. By
    repeating this operation for each $v_i$, we obtain $n +
    |\exists_{\VAR(\lu)} G| \leq |G|$. But as $\lu$ is a stick,
    trivially $|\lu| = n$, and as $R = \exists_{\VAR(\lu)} G$ by
    definition this is the required inequality.
\end{proof}
    
\ignore{    Suppose $\VAR(\lu) = \{x_1, x_2, \ldots, x_n\}$, where $x_1 \prec x_2
    \prec \cdots \prec x_n$. Define $\lu_k$ and $G_k$ to be $\lu$ and
    $G$ respectively with the variables $x_i$ for $i > k$
    existentially quantified out.
    We use induction on the sequence of variables $\VAR(\lu)$ to show
    that $|\lu_k| + |R| \leq |G_k| + 1$ for all $k \leq n$.

    If $k = 0$ the result is trivially true. As $\lu_0$ is a
    constant ROBDD, $|\lu_0| = 0$, and as $G_0$ has all variables in
    $\VAR(\lu)$ quantified out, by definition $R = G_0$. Therefore $|\lu_0| +
    |R| = 1 + |R| \leq |G_0|$.

    Now assume the result is true for $k \leq m$. By definition we have
    $\lu_{m+1} = \lu_m \wedge (x_{m+1})$ or $\lu_{m+1} = \lu_m \wedge
    (\neg x_{m+1})$. Without loss of generality we assume the former.
    Then $G_{m+1} \ra x_{m+1}$ and hence every node labelled by
    $x_{m+1}$ in $G_{m+1}$ has $0$ as its ``else'' child, and there must
    be at least one node so labelled. We then have $|G_m| + 1 \leq
    |G_{m+1}|$, since we simply omit the nodes labelled $x_{m+1}$ from
    $G_{m+1}$ to obtain $G_m$. Hence $|\lu_{m+1}| + |R| = |\lu_m| + |R|
    + 1 \leq |G_m| + 1 \leq |G_{m+1}|$. The result follows by induction.
\end{proof}
    }

Note that $|D|$ can be $O(|\lu| \times |R|)$.  For example, considering
the ROBDDs in Figure~\ref{fig:g} where
$\lu$ is shown in (a), $R$ in (b), and $D = \lu \wedge R$ in (c),
we have that $|\lu| = 5$ and $|R| = 9$ but $|D| = 9 + 4 \times 5 = 29$.

We now show how to construct a propagator on split domains.
Firstly, we eliminate any fixed variables (as in the bounds propagator)
and then apply domain propagation on the remainder domains $R$.
The propagator produces a new pair $(\lu, R)$ consisting of new fixed
variables and a new remainder. This process is shown below:
\begin{equation}
  \label{eq:split}
  \begin{split}
    \phi_0 &= B(c) \\
    \phi_j &= \exists_{\VAR(\lu(D(v_j)))} (\lu(D(v_j)) \wedge \phi_{j-1}) \\
    \delta_i &= \bexists_{V(v_i)}\,
    (\phi_n \wedge
    \bigwedge_{j=1}^n R(D(v_j))) \\
    \beta_i &= \lu(D(v_i)) \wedge \essential{\delta_i} \\
    \dom(c)(D)(v_i) &= (\beta_i, \exists_{\VAR(\beta_i)} \delta_i) \\
  \end{split}
\end{equation}

For efficiency the $\delta_i$ components can be calculated using the
divide-and-conquer approach described for the domain propagator.

The split domain representation has three main advantages. 
Proposition~\ref{prop:split-sizes} tells us that the split domain
representation is no larger than the original domain representation.
However, often the split representation is substantially smaller,
which can lead to improvements in propagation performance. The split solver
can also use the propagator simplification technique from the
bounds solver by abstracting the fixed variables out of the
propagator ROBDDs.
Finally using the split solver we can mix the usage
of domain and bounds propagators in the same problem.


Experimental results for the split domain propagator are given in
Section~\ref{sect:exp}.

\subsection{Cardinality Bounds Propagation}

Given that we are able to model a set bounds propagator using ROBDDs,
it is also appropriate to consider how we might model other
levels of consistency for set constraint problems. One level of consistency
that is commonly used
\cite{azevedo:cardinal,mueller:mozart} is that of combined set bounds and
cardinality consistency, in which upper and lower bounds 
on the cardinality of each domain are maintained in addition to the
bounds under the subset partial ordering. This hybrid approach allows a more
accurate representation
of some domains, particularly those with constrained cardinality,
which are common in set problems.

\begin{figure}
\centerline{\includegraphics[height=6cm]{cardinality} }
\caption{The set interval $[\emptyset, \{1,2,3,4\}]$
with upper and lower cardinality bounds $u=3$ and $l=2$
respectively}
\label{fig:bounds-and-cardinality}
\rule{\textwidth}{0.4pt}
\end{figure}

\begin{example}
    Figure~\ref{fig:bounds-and-cardinality} depicts the set interval
    $[\emptyset, \{1,2,3,4\}]$, together with lower and upper cardinality
    bounds 2 and 3 respectively. In general an interval consists of
    a large number of sets, making it a crude approximation to a set
    domain. Cardinality bounds permit a more fine
    grained representation by effectively allowing us to select a subset
    of the rows of the lattice diagram.
\end{example}

Just as we were able to implement a set bounds propagator using the
domain propagator together with a 
function which extracts the set bounds of a domain, so too
can we create a combined set bounds and set cardinality propagator.
Here we will extend the split domain solver by simplifying the
``remainder'' component of a split domain to an ROBDD representing its 
cardinality bounds.

As before, we need a method for extracting an ROBDD representing the
cardinality bounds of an arbitrary domain ROBDD. We perform this
operation in two stages.  Firstly we define a function
\textsf{bdd\_count\_cardinality} which takes an ROBDD representing a
set domain and returns upper and lower bounds on its cardinality.  We
can represent these bounds in ROBDD form by constructing a new
cardinality constraint ROBDD as described in
Section~\ref{sect:modeling-constraints}. 

\begin{figure}[t]
\begin{quote}
\begin{tabbing}
xx \= xx \= xx \= xx \= xx \= \kill
\textsf{bdd\_count\_cardinality}($D$, $V$): \\
\> \textbf{if} ($D = 0$) \textbf{return} $\n{\infty, -\infty}$ \\
\> \textbf{else if} ($D = 1$) \\
\> \> \textbf{if} ($|V| = 0$) \textbf{return} $\n{0, 0}$ \\
\> \> \textbf{else} \\
\> \> \> $\n{v_1, v_2, \dotsc, v_n}$ := $V$ \\
\> \> \> \textbf{return} $\n{0, n}$ \\
%%PJS\> \> \> $\n{l,u}$ := \textsf{bdd\_count\_cardinality}($D$, $\n{v_2,\dotsc, v_n}$) \\ 
%%PJS\> \> \> \textbf{return} $\n{l, u+1}$ \\
\> \textbf{else} \\
\> \> \textbf{if} ($|V| = 0$) \textbf{error} \\
\> \> $n(v,t,e)$ := $D$ \\
\> \> $\n{v_1, v_2, \dotsc, v_n}$ := $V$ \\
\> \> \textbf{if} ($v_1 \succ v$) \textbf{error} \\
\> \> \textbf{else if} ($v = v_1$) \\
\> \> \> $\n{l_t, u_t}$ := \textsf{bdd\_count\_cardinality}($t$, $\n{v_2,
\dotsc, v_n}$) \\
\> \> \> $\n{l_e, u_e}$ := \textsf{bdd\_count\_cardinality}($e$, $\n{v_2,
\dotsc, v_n}$) \\
\> \> \> \textbf{return} $\n{\min(l_t + 1, l_e), \max(u_t + 1, u_e)}$ \\
\> \> \textbf{else} \\
\> \> \> $\n{l,u}$ := \textsf{bdd\_count\_cardinality}($D$, $\n{v_2,
\dotsc, v_n}$) \\ 
\> \> \> \textbf{return} $\n{l, u+1}$ \\
\\
\textsf{bdd\_card\_bounds}($D$, $V$): \\
\> $\n{l,u}$ := \textsf{bdd\_count\_cardinality}($D$, $V$) \\
\> \textbf{if} ($l = \infty$ \textbf{ or } $u = -\infty$) \textbf{return} 0 \\
\> \textbf{return} $\cardlu(V, l, u)$ 
\end{tabbing}
\end{quote}
\caption{An algorithm to determine the cardinality bounds on the 
domain of a set variable represented by an ROBDD $D$,
where $V$ is the vector of bits in
the Boolean representation of the set variable}
\label{fig:alg-card-count}
\rule{\textwidth}{0.4pt}
\end{figure}


An implementation of \textsf{bdd\_count\_cardinality} is shown in
Figure~\ref{fig:alg-card-count}. 
This function can be implemented to run in $O(|D||V|)$ time if dynamic
programming/caching is used to save the results of the
intermediate recursive calls. In practice since $D$ and $V$ are
highly interrelated, it is $O(|D|)$. In our implementation we utilise
the global cache mechanism of the ROBDD library, which 
also permits caching of partial results between multiple calls to
\textsf{bdd\_count\_cardinality}.

We can now use this function to construct a function
\textsf{bdd\_card\_bounds} which takes an ROBDD $D$ and a list of
Boolean variables $V$ and returns a new
ROBDD representing the bounds on the cardinality of the solutions of
$D$. The implementation of \textsf{bdd\_card\_bounds} is shown in
Figure~\ref{fig:alg-card-count}.

The algorithm for a split set bounds and set cardinality propagator for a
constraint $c$ is given by $\sbc(c)$ in the following equation:
\begin{equation}
  \label{eqn:card-bounds-prop}
  \begin{split}
    \phi_0 &= B(c) \\
    \phi_j &= \exists_{\VAR(\lu(D(v_j)))} (\lu(D(v_j)) \wedge \phi_{j-1}) \\
    \delta_i &= \bexists_{V(v_i)}\,
    (\phi_n \wedge
    \bigwedge_{j=1}^n R(D(v_j))) \\
    \beta_i &= \lu(D(v_i)) \wedge \essential{\delta_i} \\
    \sbc(c)(D)(v_i) &= (\beta_i,
    \textsf{bdd\_card\_bounds}(\exists_{\VAR(\beta_i)} \delta_i,
    V(v_i) \setminus \VAR(\beta_i)))
  \end{split}
\end{equation}

Note that we only keep the cardinality bounds on the remaining non-fixed
Boolean variables for a set variable $v$
rather than on all the original variables $V(v)$, 
since we do not need to consider fixed variables again, and it leads to
a slightly smaller cardinality ROBDD.
As usual Equation~\ref{eqn:card-bounds-prop} should be implemented using a divide-and-conquer
approach for efficiency. Experimental results for this propagator are
shown in Section~\ref{sect:exp}.

\begin{example}
We illustrate set bounds and cardinality propagation
on the constraint $c \equiv lexlt(x,y)$ whose
ROBDD $B(c)$ is shown in Figure~\ref{fig:partition-and-lexlt}(b).
Assume the original domains for $x$ and $y$ are
universal, so $D(x) = D(y) = [\emptyset,\{1,2,3\}]$, represented by
the
ROBDD $\True$.  The ROBDD for $\delta_x = \bar{\exists}_{V(x)} B(c)$
is shown in Figure~\ref{fig:cp}(a), and 
$\beta_x = \llbracket \delta_x \rrbracket = \True$.
The tree of calculations for
\textsf{bdd\_card\_bounds}($\delta_x,\n{x_1,x_2,x_3}$) is shown in
Figure~\ref{fig:cp}(b). Overall the cardinality of $x$ is
determined to be in the range $[0,2]$.
\end{example}
\begin{figure}
$$
\begin{array}{ccc}
     \xymatrix@=10pt{
    \xyo{x_1}\ar@{->}[dr]\ar@/_15pt/@{-->}[dddr]\\
    &\xyo{x_2}\ar@{->}[dr]\ar@{-->}[dd]\\
    &&\xyo{x_3}\ar@{->}[dr]\ar@{-->}[dl]\\
    &*+[F]{\mathrm{\True}} && *+[F]{\mathrm{\False}}
    }
&
~
&
     \xymatrix@=10pt{
    & \n{0,2}\ar@{->}[dr]\ar@{-->}[dl]\\
      \n{0,2}&&\n{0,1}\ar@{->}[dr]\ar@{-->}[dl]\\
    & \n{0,1}&&\n{0,0}\ar@{->}[dr]\ar@{-->}[dl]\\
    &&\n{0,0} && \n{\infty,-\infty}
    }
\\
~\\
(a) && (b)
\end{array}
$$
\caption{Cardinality propagation example showing (a) resulting
ROBDD projected onto $x$, and (b) calculation of cardinality bounds}
\label{fig:cp}
\rule{\textwidth}{0.4pt}
\end{figure}
%\pjs{Give an example of propagating cardinality!}

\subsection{Lexicographic Bounds Propagation}
\label{sect:lexbounds}

An alternative form of set consistency proposed by
\citeA{sadlerg:lexbounds} is to maintain bounds under a 
lexicographic ordering in addition to set bounds. The lexicographic
ordering is a total ordering on sets which embeds the subset partial
ordering. Bounds under the lexicographic ordering alone are not sufficient
to express the effects of many constraints (in particular the
inclusion of a single element), so Sadler and Gervet constructed a
hybrid solver which combines lexicographic bounds with traditional
set bounds. They demonstrated that a set solver based
upon lexicographic bounds consistency techniques produced stronger
propagation than a traditional set bounds solver, although this came
at a substantial cost in propagation performance. However, given that
the use of the ROBDD representation leads to a performance improvement
in the case of set bounds propagation, it is worth investigating the
performance of an ROBDD-based lexicographic bounds propagator.

The lexicographic bounds of a domain can be very compactly represented
as an ROBDD. 
Like set bounds, an ROBDD representing the upper or lower
lexicographic bounds of a domain is an ROBDD of size $O(N)$,
and so is the combination.
Because these ROBDDs are compact this hopefully leads to 
fast propagation. 
Moreover, given an ROBDD domain, it is
very easy to extract the lexicographic bounds on that domain in a
single pass. By a process analogous to the construction of the bounds
and cardinality propagator, we can use a split domain propagator
combined with a function which determines the lexicographic bounds of a domain
to construct a highly efficient
lexicographic bounds solver.

\ignore{\ph{Diagram of lex bounds robdd}}

We define two functions \textsf{bdd\_lex\_lower\_bound} and
\textsf{bdd\_lex\_upper\_bound}. These functions, given an ROBDD representing
the
domain $D$ of a variable $V$ together with a list of Boolean variables
$B$ corresponding to the bits of $V$, return the lower and upper 
lexicographic bounds (respectively) of $D$.
\textsf{bdd\_lex\_lower\_bound} is implemented as shown in
Figure~\ref{fig:bdd-lex-lower-bound} (the
implementation of \textsf{bdd\_lex\_upper\_bound} is very similar).

\begin{figure}[t]
\begin{quote}
\begin{tabbing}
xx \= xx \= xx \= xx \= xx \= \kill
\textsf{bdd\_lex\_lower\_bound}($D$, $B$): \\
\> \textbf{if} ($|B| = 0$ \textbf{or} $D = 1$) \textbf{return} 1 \\
\> \textbf{if} ($D = 0$) \textbf{error} \\
\> $n(v,t,e)$ := $D$ \\
\> $\n{b_1, \dotsc, b_n}$ := $B$ \\
\> \textbf{if} ($b_1 \succ v$) \textbf{error} \\
\> \textbf{else if} ($b_1 = v$ \textbf{and} $e = 0$) \\
\> \> \textbf{return} $b_1 \wedge \textsf{bdd\_lex\_lower\_bound}(t,
\n{b_2,\dotsc,b_n})$ \\
\> \textbf{else} \\
\> \> \textbf{if} ($b_1 = v$) $r$ := \textsf{bdd\_lex\_lower\_bound}($e$,
$\n{b_2,\dotsc,b_n}$) \\
\> \> \textbf{else} $r$ := \textsf{bdd\_lex\_lower\_bound}($D$,
$\n{b_2,\dotsc,b_n}$) \\
\> \> \textbf{return} $b_1 \vee (\neg b_1 \wedge r)$ 
\end{tabbing}
\end{quote}
\caption{An algorithm to extract the lower lexicographic bounds on the
domain of a set variable represented by an ROBDD $D$, where $B$ is the
vector of (non-fixed) bits  in the Boolean representation of the set variable}
\label{fig:bdd-lex-lower-bound}
\rule{\textwidth}{0.4pt}
\end{figure}

We define:
\[
    \textsf{bdd\_lex\_bounds}(D, B) =
    \textsf{bdd\_lex\_lower\_bound}(D, B) \wedge
    \textsf{bdd\_lex\_upper\_bound}(D, B)
\]

The split set and lexicographic bounds propagator can then be
implemented exactly as in Equation~\eqref{eqn:card-bounds-prop}, using
\textsf{bdd\_lex\_bounds} in place of \textsf{bdd\_card\_bounds}.
\ignore{
\begin{equation*}
  \begin{split}
    \phi_0 &= B(c) \\
    \phi_j &= \exists_{\VAR(\lu(D(v_j)))} (\lu(D(v_j)) \wedge \phi_{j-1}) \\
    \delta_i &= \bexists_{V(v_i)}\,
    (\phi_n \wedge
    \bigwedge_{j=1}^n R(D(v_j))) \\
    \beta_i &= \lu(D(v_i)) \wedge \essential{\delta_i} \\
    \sblex(c)(D)(v_i) &= (\beta_i,
    \textsf{bdd\_lex\_bounds}(\exists_{\VAR(\beta_i)} \delta_i,
    \VAR(D(v_i)) \setminus \VAR(\beta_i)) \\
  \end{split}
\end{equation*}
}

\begin{example}
Consider lexicographic bounds propagation on the constraint
$c \equiv |s| = 2$ where $s \subseteq \{1,2,3,4,5\}$.
The ROBDD $B(c)$ is shown in Figure~\ref{fig:doms}(b)
and initial domain $D(s) = [\emptyset,\{1,2,3,4,5\}]$.  Then
$\delta_s = B(c)$ and $\beta_s = \True$ so the call
to \textsf{bdd\_lex\_bounds} calculates the lower bounds
lexicographic ROBDD shown in Figure~\ref{fig:lp}(a),
the upper bounds lexicographic ROBDD shown in Figure~\ref{fig:lp}(b),
and final answer the conjunction shown in Figure~\ref{fig:lp}(c).
Note that we have lost some information relative to the original cardinality
ROBDD.
\end{example}
\begin{figure}
$$
\begin{array}{ccccc}
    \xymatrix@=10pt{
      &&& \xyo{s_1}\ar@{->}@/^25pt/[dddddl]\ar@{-->}[dl]\\
      && \xyo{s_2}\ar@{->}@/^15pt/[dddd]\ar@{-->}[dl] \\
      &\xyo{s_3}\ar@{->}[dddr]\ar@{-->}[dl] \\
      \xyo{s_4}\ar@{->}[dr]\ar@{-->}@/_35pt/[dddrr] \\
 & \xyo{s_5}\ar@{->}[dr]\ar@{-->}@/_25pt/[ddr] \\
      && *+[F]{\mathrm{\True}}            \\
      && *+[F]{\mathrm{\False}} 
    }
&&
    \xymatrix@=10pt{
      & \xyo{s_1}\ar@{->}[dr]\ar@{-->}@/_15pt/[dddddl]\\
       &&
      \xyo{s_2}\ar@{->}[dr]\ar@{-->}[ddddll]\\
      &
      &&\xyo{s_3}\ar@{-->}[dl]\ar@{->}@/^35pt/[ddddlll] \\     
      && \xyo{s_4}\ar@{-->}[dl]\ar@{->}@(dr,r)@/^25pt/[dddll]   &&\\
 & 
\xyo{s_5}\ar@{-->}[dl]\ar@{->}@(dr,r)@/^15pt/[ddl]
&&&&&\\
       *+[F]{\mathrm{\True}}         &&&   \\
       *+[F]{\mathrm{\False}} &&&&
    }
&&
    \xymatrix@=10pt{
      &&& \xyo{s_1}\ar@{->}[dr]\ar@{-->}[dl]\\
      && \xyo{s_2}\ar@{->}@/^15pt/[dddd]\ar@{-->}[dl] &&
      \xyo{s_2}\ar@{->}[dr]\ar@{-->}[ddddll]\\
      &\xyo{s_3}\ar@{->}[dddr]
      \ar@{-->}[dl]
      &&
      &&\xyo{s_3}\ar@{-->}[dl]\ar@{->}@/^35pt/[ddddlll] \\
      \xyo{s_4}\ar@{->}[dr]\ar@{-->}@/_35pt/[dddrr]
      &&
      && \xyo{s_4}\ar@{-->}[dl]\ar@{->}@(dr,r)@/^25pt/[dddll]   &&\\
 & \xyo{s_5}\ar@{->}[dr]\ar@{-->}@/_25pt/[ddr]
      &&
\xyo{s_5}\ar@{-->}[dl]\ar@{->}@(dr,r)@/^15pt/[ddl]
&&&&&\\
      && *+[F]{\mathrm{\True}}         &&&   \\
      && *+[F]{\mathrm{\False}} &&&&
    } \\
~\\
(a) && (b) && (c) 
\end{array}
$$
\caption{Calculating lexicographic (a) lower and (b) upper bounds
and (c) their conjunction.}
\label{fig:lp}
\rule{\textwidth}{0.4pt}
\end{figure}

%\pjs{Give an example of propagating lex bounds!}


As observed in Section~\ref{sect:int-cons}, the lexicographic ordering
for set variables actually corresponds to a numeric ordering on integer variables,
so a pure lexicographic bounds propagator would also be coincidentally an
integer bounds propagator.

Experimental results for the lexicographic bounds propagator are given
in Section~\ref{sect:exp}.

\section{Experiments}
\label{sect:exp}

We have implemented a set solver using the ideas described in this
paper. Our implementation is written primarily in \textsc{Mercury}
\cite{mercury}, using an interface to the C language ROBDD library
\textsc{CUDD} as a platform for ROBDD manipulations \cite{cudd}.
\ignore{It should be noted that the ROBDDs as used by the CUDD library do not
correspond precisely to the definition of ROBDDs as presented in this
paper. The ROBDDs implemented by CUDD incorporate an extension known
as complement arcs, which allow a more concise expression of the
some ROBDDs.
This does not substantially affect the nature 
of our solver since ROBDDs are effectively treated as a `black-box'.}
The ROBDD library is effectively treated as a black box.

We used our solver to implement a series of standard constraint benchmarks as
described below. Many of these problems are in the CSPLib library of
constraint satisfaction problems \cite{csplib}. For the purposes of
comparison, we also implemented our benchmarks using the
\texttt{ic\_sets} library of \eclipse{} v5.6\footnote{Very
recently a new sets library \texttt{Cardinal} has been added to \eclipse{}
which supports better cardinality reasoning. Unfortunately
we cannot directly compare with it on these benchmarks since it does
not support lexicographic ordering constraints. 
Testing without lexicographic orderings showed that it
was about 2--3 times slower than LU + Card.} 
and the finite sets
library of \textsc{Mozart} v1.3.0.  We conducted all of our experiments on a
cluster of 8 identical 2.8GHz Pentium 4 machines, each with 1 Gb of
RAM and running Debian GNU/Linux 3.1. All three solvers were limited
to 1 Gb of memory to minimise swapping. All experiments were repeated
three times, and the lowest time out of the three runs was taken as
the result. 

\subsection{Steiner Systems}
\label{sect:exp-steiner}

%\pjs{It would be nice to show with an  experiment or at least
%comment on the speed of the three different implementations
%of dom(c): naive, linear, logarithmic}

A 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. Steiner
systems have been extensively studied in combinatorial mathematics. If
$t=2$ and $k=3$, we have the so-called Steiner triple systems, which
are often used as benchmarks
\cite{gervet:conjunto,azevedo:cardinal,mueller:mozart}. Any Steiner
system must have exactly $m={N \choose t}/{k \choose t}$ blocks
(Theorem 19.2 of \citeR{vanlintw:combinatorics}).

It is natural to model a Steiner system using set variables $s_1,
\ldots, s_m$, where each set variable corresponds to a single block,
subject to the following constraints:
\begin{align}
    \label{eqn-steiner-unary}
    \bigwedge_{i=1}^m (|s_i| = k) \\
    \wedge
    \label{eqn-steiner-binary}
    \bigwedge_{i=1}^{m-1} \bigwedge_{j=i+1}^{m}
    (\exists u_{ij} ~ u_{ij} = s_i \cap s_j \wedge |u_{ij}| \leq (t-1))
    \wedge (s_i < s_j) 
\end{align}
The lexicographic ordering constraint $s_i < s_j$ has been added to 
remove symmetries from the problem formed by permuting the blocks.

A necessary condition for the existence of a Steiner system is that
${{N-i} \choose {t-i}}/{{k-i} \choose {t-i}}$ is an integer for all $i
\in \{0,1,\ldots,t-1\}$ \cite{vanlintw:combinatorics}; we say a set of
parameters $(t,k,N)$ is admissible if it satisfies this condition. In
order to choose test cases, we ran each solver on every admissible set
of $(t,k,N)$ values for $N < 32$. Results are shown for every test
case that at least one solver was able to solve within a time limit of
10 minutes. ``$\times$'' denotes abnormal termination due to exceeding
an arbitrary limit on the maximum number of ROBDD variables imposed 
by the CUDD package, 
while ``---'' denotes failure to
complete a testcase within the time limit.

In all cases we use a sequential variable-ordering
heuristic, and a value-ordering heuristic that chooses the largest
unfixed value within a variable's domain. Once a variable and an
unfixed value
have been chosen for labeling, either that value is a member of the set
represented
by the variable or it is not. The order in
which we try the two alternatives has a significant effect on the
performance of the solver. In this case we elect to choose the
element-not-in-set option first.

\begin{table}[p]
\caption{Results for Steiner systems, with split
constraints and intermediate variables.}
%$\times$ denotes abnormal
%termination on a test-case.}
\label{tab:steiner-split}
\begin{tabular}{||l||r|r||r|r||r|r||r|r||r|r||r|r||r|r||}
\hhline{|t:=:t:==:t:==:t:==:t:==:t:==:t:==:t:==:t|}
%Testcase&\eclipse&&\textsc{Mozart}&&Bounds&&Domain&&LU+R&&LU+Lex&&LU+Card& \\
Testcase &
\multicolumn{2}{c||}{\eclipse} &
\multicolumn{2}{c||}{\sc Mozart} &
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{2}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card}  \\
\hhline{||~||--||--||--||--||--||--||--||}
&Time&Fails&Time&Fails&Time&Fails&Time&Fails&Time&Fails&Time&Fails&Time&Fails
\\[-3pt]
&
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} & & 
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||-||--||--||--||--||--||--||--||}
S(2,3,7)&0.3&10&0.1&21&\bf $<$0.1&10&0.1&\bf 0&0.1&\bf 0&0.1&4&\bf $<$0.1&2 \\
S(3,4,8)&0.5&21&\bf 0.1&52&\bf 0.1&21&0.4&\bf 0&0.4&\bf 0&0.4&4&\bf 0.1&4 \\
S(2,3,9)&7.7&1394&\bf 1.0&5102&1.6&1394&\bf 1.0&\bf 100&1.3&\bf 100&3.3&421&2.4&1072
\\
S(2,4,13)&1.8&313&\bf 0.4&1685&0.6&313&1.7&\bf 32&1.5&\bf 32&2.1&127&0.7&157 \\
S(2,3,15)&3.6&65&\bf 0.5&354&2.2&65&20.2&\bf 0&19.6&\bf 0&20.4&127&3.3&41 \\
S(3,4,16)&\bf 67.5&\bf 289& --- & --- & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$  \\
S(2,5,21)&3.2&421&\bf 0.4&668&2.3&421&110.2&\bf 0&59.8&\bf 0&21.5&139&2.6&124 \\
S(3,6,22)&\bf 49.7& \bf 1619& --- & ---& $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$ & $\times$  \\
\hhline{|b:=:b:==:b:==:b:==:b:==:b:==:b:==:b:==:b|}
\end{tabular}
\end{table}

\begin{table}[p]
\caption{Results for Steiner systems, with merged
constraints and no intermediate variables.}
%`---' denotes
%failure to complete a testcase within 10 minutes.}
\label{tab:steiner-merged}
\begin{tabular}{||l||r|r||r|r||r|r||r|r||r|r||}
\hhline{|t:=:t:==:t:==:t:==:t:==:t:==:t|}
%Testcase&Eclipse&&Mozart&&Bounds&&Domain&&LU+R&&LU+Lex&&LU+Card& \\
Testcase &
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{2}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card} \\
\hhline{||~||--||--||--||--||--||}
&Time&Fails&Time&Fails&Time&Fails&Time&Fails&Time&Fails
\\[-3pt]
&
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} & & 
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} & &
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||-||--||--||--||--||--||}
S(2,3,7)&\bf $<$0.1&8&\bf $<$0.1&\bf 0&\bf $<$0.1&\bf 0&\bf $<$0.1&\bf 0&\bf $<$0.1&\bf 0 \\
S(3,4,8)&\bf 0.1&18&\bf 0.1&\bf 0&\bf 0.1&\bf 0&\bf 0.1&\bf 0&\bf 0.1&\bf 0 \\
S(2,3,9)&0.2&325&\bf 0.1&\bf 9&\bf 0.1&\bf 9&\bf 0.1&11&0.2&113 \\
S(2,3,13)&---&---&\bf 109.2&\bf 24723&144.6&\bf 24723&518.1&30338&---&--- \\
S(2,4,13)&\bf 0.1&157&\bf 0.1&\bf 0&\bf 0.1&\bf 0&0.4&11&\bf 0.1&27 \\
S(2,3,15)&\bf 0.4&56&1.3&\bf 0&1.4&\bf 0&2.9&\bf 0&0.7&32 \\
S(2,4,16)&421.4&522706&\bf 0.6&\bf 15&\bf 0.6&\bf 15&2.5&16&577.0&209799 \\
S(2,6,16)&---&---&\bf 80.7&\bf 15205&82.7&\bf 15205&---&---&---&--- \\
S(3,4,16)&\bf 9.7&274&548.7&\bf 0&485.3&\bf 0&428.9&\bf 0&18.4&162 \\
S(2,5,21)&\bf 0.5&413&1.4&\bf 0&1.4&\bf 0&32.9&\bf 0&0.6&116 \\
S(3,6,22)&\bf 8.3&1608&---&---&---&---&---&---&12.7&\bf 381 \\
S(2,3,31)&\bf 23.3&280&---&---&---&---&---&---&48.6&\bf 224 \\
\hhline{|b:=:b:==:b:==:b:==:b:==:b:==:b|}
\end{tabular}
\end{table}

\begin{table}[p]
\caption{All-solutions results on Steiner systems. ``---''
denotes failure to complete a test case within one hour}
\label{tab:steiner-allsolns}
\begin{tabular}{||l||r||r|r||r|r||r|r||r|r||r|r||r|r||}
\hhline{|t:=:t:=:t:==:t:==:t:==:t:==:t:==:t:==:t|}
\multirow{3}{33pt}{Problem} &
\multirow{3}{28pt}{Solns.} &
\multicolumn{2}{c||}{\eclipse} &
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{2}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card}
\\
\hhline{||~||~||--||--||--||--||--||--||}
&&
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails}
\\[-3pt]
&&
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} &&
\multicolumn{1}{c|}{\scriptsize /s} &&
\multicolumn{1}{c|}{\scriptsize /s} & 
\\
\hhline{||-||-||--||--||--||--||--||--||}
S(2,3,7) & 30 & 16.3 & 3,015 & 0.2 & 537 & \bf 0.1 & \bf 47 & \bf 0.1
& \bf 47 & 0.2 & 76 & 0.3 & 267 \\
S(3,4,8) & 30 & --- & --- & 726.5 & 610271 & \bf 1.4 & \bf 492 & 2.1
& \bf 492 & 22.4 & 3248 & 951.2 & 431801 \\
S(2,3,9) & 840 & --- & --- & 398.1 & 391691 & \bf 23.4 & \bf 16794 &
37.8 & \bf 16794 & 110.4 & 29133 & 593.3 & 224131 \\
S(2,6,16) & 0 & --- & --- & --- & --- & \bf 80.9 & \bf 15205 & 83.0 & \bf 15205 & --- & --- & --- & --- \\
\hhline{|b:=:b:=:b:==:b:==:b:==:b:==:b:==:b:==:b|}
\end{tabular}

\end{table}


In order to compare the raw performance of the various solvers,
irrespective of any modeling advantages of the ROBDD-based solver, we
performed experiments using a model of the problem which contains only
primitive constraints and makes use of intermediate variables. This
``split'' model contains $m$ unary constraints corresponding to
Equation~\eqref{eqn-steiner-unary} and $3m(m-1)/2$ binary
constraints corresponding to Equation~\eqref{eqn-steiner-binary}
(containing intermediate variables $u_{ij}$). The same model was used
in the \eclipse{}, \textsc{Mozart} and ROBDD-based solvers, permitting a direct
comparison of propagation performance. Results for this model are
shown in Table~\ref{tab:steiner-split}. In particular, observe that
with this model \eclipse{} and the ROBDD-based bounds solver produce the
same number of failures, demonstrating that the search spaces explored
by the two solvers are identical.

One of the main strengths of the ROBDD-based modeling approach is
that it gives us the freedom to merge arbitrary constraints and
existentially quantify away intermediate variables, allowing us to
model set constraint problems more efficiently. In the case of
Steiner systems, this allows us to model the problem as $m$ unary
constraints corresponding to Equation~\eqref{eqn-steiner-unary}, and
just $m(m-1)/2$ binary constraints $\psi_{ij}$ of the form
$\psi_{ij} = (|s_i \cap s_j| \leq (t-1)) \wedge (s_i < s_j)$. The
binary constraints do not contain the intermediate variables 
$u_{ij}$---they are not required since they can be existentially quantified
out of the ROBDD representation.  Experimental results for this
revised model are shown in Table~\ref{tab:steiner-merged}. 

In all
cases the revised model propagates much more strongly than the
original model, leading to a substantial decrease in solution time. In
addition, the decrease in the number of set variables required permits
the solution of larger test cases. Clearly it is beneficial to remove
intermediate variables and merge constraints.


Despite weaker propagation the ROBDD bounds solver is often the
fastest method of finding a single solution to a Steiner System. In
order to determine whether this was due to the efficiency of the
solver, or whether the solver was just ``lucky'' in finding a solution
quickly, we also ran experiments to find all solutions of the Steiner
systems up to reordering of the blocks.
The results for all test cases that at least one of the solvers
was able to solve within a time limit of one hour are shown in
Table~\ref{tab:steiner-allsolns}. In all cases
the reduction both of time and number of fails demonstrate the
superiority of the propagation approaches based on domain consistency.

\subsection{Social Golfers}
\label{sect:exp-golfers}

Another problem often used as a set benchmark is the ``Social
Golfers'' problem (problem \texttt{prob010} of CSPLib). The aim of
this problem is to arrange  $N = g \times s$ golfers into $g$ groups
of $s$ players for each of $w$ weeks, such that no two players play
together more than once. We can model this problem as a set constraint
problem using a $w \times g$ matrix of set variables $v_{ij}$, where
$1 \leq i \leq w$ is the week index and $1 \leq j \leq g$ is the group
index.

We use the following model of the problem:
\begin{equation}
\begin{split}
\bigwedge_{i=1}^w (\partition^{<}(v_{i1}, \ldots, v_{ig})) 
\wedge \bigwedge_{i=1}^w \bigwedge_{j=1}^g |v_{ij}| = s  \\
\wedge \bigwedge_{\substack{i,j \in \{1,\ldots,w\} \\  i \neq j}} \, 
       \bigwedge_{k,l \in \{1,\ldots,g\}}
       |v_{ik}\cap v_{jl}| \leq 1 
\wedge \bigwedge_{i=1}^{w-1} \bigwedge_{j=i+1}^w v_{i1} \leq v_{j1}
\end{split}
\end{equation}

The $\texttt{partition}^{<}$ global constraint is a combined
partitioning and lexicographical ordering constraint, formed by merging
the $\texttt{partition}$ constraint of Section~\ref{sect:global} with
constraints imposing a lexicographic order on the variables. This
constraint is trivial
to construct using ROBDDs, but is not available in either \eclipse{} or 
\textsc{Mozart}.\footnote{It would, of course, be possible to implement such a
constraint in both \eclipse{} and \textsc{Mozart}, but such an implementation
would be a fairly laborious process. A strength of the
ROBDD-based modeling approach 
is that we can construct global constraints with no extra code.}

\begin{table}[t!]
\caption{First-solution performance results on the Social Golfers
problem, using a sequential ``smallest-element-in-set'' labeling
strategy. Time and number of failures are given for all solvers.
``---'' denotes failure to complete a test case within 10 minutes.
The cases 5-4-3, 6-4-3, and 7-5-5 have no solutions}
\label{tab:golfer}
\begin{tabular}{||l||r|r||r|r||r|r||r|r||r||r|r||r|r||}
\hhline{|t:=:t:==:t:==:t:==:t:==:t:=:t:==:t:==:t|}
\multirow{2}{41pt}{Problem} &
\multicolumn{2}{c||}{\eclipse}&
\multicolumn{2}{c||}{\sc Mozart}&
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{1}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card} 
\\
\hhline{||~||--||--||--||--||-||--||--||}
& 
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c||}{time} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} 
\\[-3pt]
$w$-$g$-$s$&
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c||}{\scriptsize /s} &
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||-||--||--||--||--||-||--||--||}
2-5-4&7.6&10468&1.0&7638&\bf 0.1&30&\bf 0.1&\bf 0&\bf 0.1&\bf 0.1&3&\bf 0.1&5 \\
2-6-4&49.2&64308&6.4&42346&0.6&2036&0.2&\bf 0&\bf 0.1&0.7&194&0.4&326 \\
2-7-4&95.1&114818&10.9&66637&1.7&4447&\bf 0.4&\bf 0&\bf 0.4&2.3&692&1.6&1608 \\
2-8-5&---&---&---&---&---&---&2.0&\bf 0&\bf 1.6&---&---&---&--- \\
3-5-4&12.5&14092&2.5&10311&\bf 0.1&30&0.3&\bf 0&0.3&0.4&3&0.2&5 \\
3-6-4&76.3&83815&14.0&51134&1.6&2039&1.6&\bf 0&\bf 1.4&2.2&194&1.9&328 \\
3-7-4&146.8&146419&27.3&88394&\bf 4.6&4492&8.9&\bf 0&8.4&7.6&695&5.1&1629 \\
4-5-4&14.1&14369&3.9&10715&\bf 0.2&30&0.8&\bf 0&0.6&0.7&3&0.4&5 \\
4-6-5&---&---&---&---&21.9&12747&118.6&\bf 0&80.7&\bf 19.3&499&32.0&2122 \\
4-7-4&169.3&149767&46.4&90712&\bf 8.7&4498&---&---&481.6&14.1&\bf 696&10.5&1632 \\
4-9-4&27.3&19065&7.8&12489&\bf 2.6&71&---&---&---&22.2&\bf 8&8.9&33 \\
5-4-3&---&---&---&---&113.9&63642&\bf 28.6&\bf 5165&32.1&88.9&10210&202.7&50542 \\
5-5-4&350.6&199632&217.7&416889&7.0&2686&3.8&\bf 41&\bf 2.3&12&313&20.8&1584 \\
5-7-4&---&---&---&---&\bf 14.6&4583&---&---&---&23.7&\bf 700&17.5&1683 \\
5-8-3&5.0&2229&\bf 0.9&1820&1.1&14&9.2&\bf 0&7.8&3.2&3&2.1&4 \\
6-4-3&---&---&---&---&158.2&61770&\bf 20.3&\bf 2132&23.0&60.8&4506&293.5&49966 \\
6-5-3&458.9&240296&287.4&471485&4.1&1455&1.8&\bf 82&\bf 1.5&4.1&202&8.3&1078 \\
6-6-3&3.3&1462&1.0&1462&\bf 0.5&5&1.8&\bf 0&1.4&1.2&0&0.7&1 \\
%7-5-3&---&---&---&---&---&---&---&---&---&---&---&---&--- \\
7-5-5&---&---&---&---&---&---&0.5&\bf 0&\bf 0.4&2.4&0&1.3&22 \\
\hhline{|b:=:b:==:b:==:b:==:b:==:b:=:b:==:b:==:b|}

\end{tabular}
\end{table}

\begin{table}[t!]
\caption{First-solution performance results on the Social Golfers
problem, using a first-fail ``smallest-element-in-set'' labeling
strategy. Time and number of failures are given for all solvers.
The cases 5-4-3, 6-4-3, and 7-5-5 have no solutions}
\label{tab:golferff}
\begin{tabular}{||l||r|r||r|r||r|r||r|r||r||r|r||r|r||}
\hhline{|t:=:t:==:t:==:t:==:t:==:t:=:t:==:t:==:t|}
\multirow{2}{41pt}{Problem} &
\multicolumn{2}{c||}{\eclipse}&
\multicolumn{2}{c||}{\sc Mozart}&
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{1}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card} 
\\
\hhline{||~||--||--||--||--||-||--||--||}
& 
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c||}{time} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} 
\\[-3pt]
$w$-$g$-$s$&
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c||}{\scriptsize /s} &
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||-||--||--||--||--||-||--||--||}
2-5-4 & 7.9 & 10468 & 1.1 & 7638 & \bf 0.1 & 30 & \bf 0.1 & \bf 0 &
\bf 0.1 & \bf 0.1 & 3 & 0.4 & 447 \\
2-6-4 & 51.3 & 64308 & 6.5 & 42346 & 0.6 & 2036 & 0.2 & \bf 0 & \bf 0.1 & 0.7 & 186 & 3.8 & 3820 \\
2-7-4 & 99.9 & 114818 & 11.1 & 66637 & 1.7 & 4447 & \bf 0.4 & \bf 0
& \bf 0.4 & 2.0 & 390 & 4.6 & 6424 \\
2-8-5 & --- & --- & --- & --- & --- & --- & 1.9 & \bf 0 & \bf 1.6 & --- & --- & --- & --- \\
3-5-4 & 14.5 & 14092 & 2.6 & 10311 & \bf 0.1 & 44 & 0.3 & \bf 0 & 0.3 & 0.4 & 5 & 1.3 & 481 \\
3-6-4 & 91.8 & 83815 & 15.1 & 51134 & 1.9 & 2361 & 1.1 & \bf 0 &
\bf 0.9 & 2.4 & 209 & 13.7 & 3722 \\
3-7-4 & 183.0 & 146419 & 28.7 & 88394 & 5.2 & 5140 & 1.9 & \bf 0 &
\bf 1.7 & 6.2 & 512 & 15.6 & 6067 \\
4-5-4 & 18.0 & 14369 & 4.1 & 10715 & \bf 0.3 & 47 & 0.8 & \bf 0 & 0.6 & 0.7 & 7 & 2.6 & 494 \\
4-6-5 & --- & --- & --- & --- & \bf 38.5 & 19376 & 62.5 & \bf 0 & 40.1 & 24.6 & 607 & --- & --- \\
4-7-4 & 243.9 & 149767 & 49.6 & 90712 & 9.9 & 5149 & 6.5 & \bf 0 &
\bf 5.1 & 10.6 & 405 & 29.5 & 5717 \\
4-9-4 & 40.9 & 19065 & 8.3 & 12489 & \bf 2.7 & 143 & 152.0 & \bf 0 & 107.4 & 22.6 & 13 & 12.5 & 516 \\
5-4-3 & --- & --- & --- & --- & 187.5 & 103972 & \bf 23.2 & \bf 3812 & 26.0 & 91.4 & 10422 & 309.1 & 72669 \\
5-5-4 & 394.9 & 199632 & 224.5 & 416889 & 4.5 & 2388 & 2.5 & \bf 18
& \bf 1.7 & 23.4 & 776 & 41.4 & 4730 \\
5-7-4 & --- & --- & --- & --- & 17.6 & 5494 & 18.2 & \bf 0 & \bf 12.8 & 16.9 & 447 & 47.1 & 5473 \\
5-8-3 & 5.4 & 2229 & \bf 0.9 & 1820 & 1.1 & 19 & 4.5 & \bf 0 & 3.9 & 3.4 & 2 & 2.7 & 83 \\
6-4-3 & --- & --- & --- & --- & 234.2 & 90428 & \bf 14.6 & \bf 1504 & 15.2 & 54.0 & 4013 & 349.9 & 59805 \\
6-5-3 & 501.0 & 240296 & 294.2 & 471485 & 1.6 & 495 & 1.2 & \bf 34 & \bf 1.0 & 58.2 & 3787 & 473.6 & 68673 \\
6-6-3 & 3.6 & 1462 & \bf 1.0 & 1462 & --- & --- & 1.6 & \bf 7 & 1.3 &
5.1 & 292 & \bf 1.0 & 8 \\
7-5-3 & --- & --- & --- & --- & --- & --- & 16.9 & \bf 528 & \bf 13.1 & 288.8 & 9829 & --- & --- \\
7-5-5&---&---&---&---&---&---&0.5&\bf 0&\bf 0.4&2.4&0&1.3&22 \\
\hhline{|b:=:b:==:b:==:b:==:b:==:b:=:b:==:b:==:b|}
\end{tabular}
\end{table}

Results for each of the solvers are shown in
Table~\ref{tab:golfer} and Table~\ref{tab:golferff}. In the former
table,
a sequential
smallest-element-in-set labeling strategy was used to enable a fair
comparison of propagation performance, whereas in the latter table a
first-fail labeling strategy was used in order to give a measure of
the peak performance of each solver. For every test case in both tables,
with a single exception (5-8-3), at least one of the
ROBDD-based solvers performs equal or better to both \eclipse{} and
\textsc{Mozart}. 
It should also be observed that when using a first-fail labeling
strategy, the domain and split domain solvers are the only solvers
able to solve every test case.

There are several other features of the results that are worth noting.
As in the case of Steiner systems, the ROBDD-based set bounds solver
is often the fastest, despite weak propagation. Amongst the solvers
with stronger propagation, the split domain solver is almost always
faster than the original domain solver due to smaller domain sizes. It
is, however, slower than the original domain solver in the presence of
backtracking (due to the
requirement to trail more values---in particular the propagator
ROBDDs). The lexicographic bounds solver is almost as effective as the
domain solvers in restricting search space, although it is usually
outperformed by the domain and bounds solvers.

\subsection{Weighted Hamming Codes}
\label{sect:exp-hamming}

The problem of finding maximal Hamming Codes can be modelled as a set
constraint problem.

We define an $l$-bit codeword to be a bit-string (or vector
of Boolean values) of length $l$. Given two $l$-bit codewords $A$
and $B$, we
define the \emph{Hamming distance} $d(A,B)$ between $A$ and $B$ to be
the number of positions at which the two bit-strings differ. 
An \emph{$(l,d)$-Hamming Code} is a set of $l$-bit codewords
such that the Hamming distance between any two codewords in the set 
is at least $d$.

Given a codeword length $l$ and the minimum Hamming distance $d$,
the problem is to construct a Hamming code with the largest possible
number of codewords. A variant of this problem, used as a benchmark by
\citeA{sadlerg:lexbounds}, has the additional requirement that each
codeword have a fixed weight $w$, where the weight of a codeword is
defined to be the number of ``$1$'' bits that codeword contains. We will
denote an instance of this problem by $H(l,d,w)$.

As proposed by \citeA{muellerm:fsetsoz}, we can model this problem for
$n$ codewords using $n$ set variables $S_i$, where $1 \leq i \leq n$
. A codeword $C_i$ corresponds to the characteristic
function of the set $S_i$, i.e. bit $j$ is set in codeword $C_i$ if and
only if $j \in S_i$. The Hamming distance $d(C_i, C_j)$ between two
codewords $C_i$
and $C_j$ can be
calculated from the associated sets $S_i$ and $S_j$ thus:
\[
d(C_i, C_j) = l - |S_i \cap S_j| - |\{1, \ldots, l\} \setminus (S_i
\cup S_j)|
\]

We can remove symmetries created by permuting the codewords by
introducing lexicographic ordering constraints $S_i < S_j$, for all $1
\leq i < j \leq n$. The complete model of the problem is:
\begin{align}
\bigwedge_{i=1}^{n} |S_i| &= w \\
\label{eqn-hamming-binary}
\wedge ~ \bigwedge_{i=1}^{n-1} \bigwedge_{j=i+1}^{n}
    (|S_i \cap S_j| + |\overline{(S_i \cup S_j)}| \leq l - d) \wedge
    (S_i &< S_j)
\end{align}

The constraints described by Equation~\eqref{eqn-hamming-binary} are
implemented using a single ROBDD for each pair of $i$ and $j$ values.
This is possible since we can model the integer addition and comparison
operations
as ROBDDs as described in Section~\ref{sect:ints},
using the representation of the cardinality of a set variable
as an integer expression as described in Section~\ref{sect:wsum}.

In order to find an optimal solution, we initially set $n$ to
$1$, and repeatedly solve instances of the problem, progressively
incrementing $n$
to find larger and larger codes. We prove optimality of a solution for
$n = k$ by failing to solve the problem for $n = k + 1$; we then know that the
optimal value for $n$ was $k$.

In order to assist with timing comparisons, we use the same set of
problem instances as \citeA{sadlerg:lexbounds}. We consider sets of
values $H(l,d,w)$ where $l \in \{6,7,8,9,10\}$, $d \in
\{4,6,8,10,12\}$, and $w \in \{3,4,5,6,7,8\}$, with $d < l$ and $w \leq
l$ (trivially there is at most one solution if $d \geq l$ and none if
$w > l$). There are $62$ such testcases; some of these are
almost identical---in particular testcases $H(l,d,w)$ and
$H(l,d,l-w)$ have solutions that differ only through complementation
of the bits. The ROBDD-based solver can solve all but seven test cases
(namely $H(9,4,4),
H(9,4,5), H(10,4,3), H(10,4,4), H(10,4,5), H(10,4,6), H(10,4,7)$), which in
reality contains three pairs of mirror image testcases.

Since there are too many results to list each testcase individually, 
performance statistics for the testcases that all the solvers were able to
solve are shown in Table~\ref{tab:hamming-stats}. Those cases that
were solved by at least one but not all of the ROBDD-based solvers are
shown in Table~\ref{tab:hamming-tricky-cases}.

\begin{table}[t]
\caption{Statistics for the 51 Weighted Hamming Code testcases that
were solved by all of the solvers.}
\label{tab:hamming-stats}
\begin{tabular}{||l||r|r||r|r||r|r||r|r||r|r||}
\hhline{|t:=:t:==:t:==:t:==:t:==:t:==:t|}
&
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{2}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card} 
\\
\hhline{||~||--||--||--||--||--||}
& 
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} 
\\[-3pt]
&
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||-||--||--||--||--||--||}
Mean & 17.7 & 110034 & \bf 0.2 & \bf 210.7 & 0.3 & \bf 210.7 & 4.4 & 1886 & 3.6
& 6604 \\
\hhline{||-||--||--||--||--||--||}
Total & 903.3 & & \bf 11.4 & & 14.2 & & 222.6 & & 184.6 & \\
\hhline{||-||--||--||--||--||--||}
Minimum &   \bf 0.03 &  \bf 0 & \bf 0.03  & \bf 0 & \bf 0.03  & \bf 0
& \bf 0.03 &  \bf 0 & \bf 0.03  & \bf 0 \\
25th Percentile  &  \bf 0.03  & 19   &\bf 0.03  & \bf 0 & \bf 0.03 &
\bf 0 & \bf 0.03  & 3 & \bf 0.03 &  \bf 0 \\
Median   &  0.05  & 196&\bf 0.04 &  \bf 2 & \bf 0.04  & \bf 2 & 0.04
& 20 &\bf 0.04  & \bf 2 \\
75th Percentile  &  0.67  & 5604 &  \bf 0.06 &  \bf 25 &\bf 0.06 &
\bf 25& 0.18 &  143.5 & 0.09  & 112 \\
Maximum  &  415.55 &3021057&\bf 4.16 &  \bf 3740 &  4.69 &  \bf 3740 &  113.9
& 45667 & 92.09 & 211677 \\
\hhline{|b:=:b:==:b:==:b:==:b:==:b:==:b|}
\end{tabular}
\end{table}

\begin{table}
\caption{Weighted Hamming Code testcases that were solved by at least
one but not all solvers}
\label{tab:hamming-tricky-cases}
\begin{tabular}{||l||r|r||r|r||r|r||r|r||r|r||}
\hhline{|t:=:t:==:t:==:t:==:t:==:t:==:t|}
&
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{2}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card} 
\\
\hhline{||~||--||--||--||--||--||}
& 
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} 
\\[-3pt]
&
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||-||--||--||--||--||--||}
H(8,4,4)  & ---&---&\bf 1.6  & \bf 224&1.9 &  \bf 224&8.5 &  986&---&---\\
H(9,4,3)  & ---&---&\bf 11.3 & \bf 5615  & 20.6 &  \bf 5615  & ---&---&---&---\\
H(9,4,6)  & ---&---&\bf 25.4 & \bf 16554 & 45.7 &  \bf 16554 & 321.9 & 56599 &
---&--- \\
H(10,6,5) & ---&---&\bf 26.7 & \bf 16635 & 29.6 & \bf 16635 & 528.8 & 169457 &
428.3 &762775  \\
\hhline{|b:=:b:==:b:==:b:==:b:==:b:==:b|}
\end{tabular}
\end{table}

\citeA{sadlerg:lexbounds} report results for this problem using set
bounds and lexicographic bounds solvers implemented in \eclipse{}.
Their solvers were able to solve 50 of
the testcases with a time limit of 240 seconds for each testcase.
Some individual testcases took upwards of 100 seconds. By contrast,
the ROBDD-based domain solver is capable of solving 55 testcases in
76.4 seconds in total. Clearly in this case enforcing domain consistency
brings a considerable reduction in search space, leading to a highly
efficient solver.

Moreover, the set bounds and lexicographic bounds solvers implemented
using the ROBDD platform appear to perform considerably better than
those of Sadler and Gervet. Due to the graphical nature of their
presentation it is difficult to quantify this performance difference; 
exact figures are presented for only two testcases---namely
$H(9,4,7)$ and $H(10,6,7)$. For the testcase $H(9,4,7)$ the ROBDD
set bounds and lexicographic bounds solvers each found and proved an
optimal solution in 0.6 and 0.3 seconds
respectively, compared with the
$>240$ and $167.1$ seconds respectively quoted for the bounds and
lexicographic solvers of Sadler and Gervet. Similarly for $H(10,6,7)$
the ROBDD set bounds and lexicographic bounds solvers found and proved
an optimal solution in 1.9 and 1.1
seconds respectively, as compared with $>240$ and 98.5 seconds.
Given this dramatic performance difference, it appears the
additional modeling flexibility of the ROBDD-based solver provides a
substantial performance gain.

It should be noted that the performance results reported by Sadler and
Gervet were run on a slightly slower machine (a 2 Ghz Pentium 4
machine). Nonetheless, the contribution to the performance difference
due to machine speed is dwarfed by the performance gap between the two
solvers.

\subsection{Balanced Academic Curricula}
\label{sect:exp-bacp}

The Balanced Academic Curriculum problem (problem \texttt{prob030} of
CSPLib) involves planning an academic curriculum over a sequence of
academic periods in order to provide a balanced load in each period.

A curriculum consists of $m$ courses ($1 \leq i \leq m$) and $n$
academic periods. Each course $i$ has a set of prerequisites and an
associated academic load $t_i$.  Every course must be assigned to
exactly one period. In any given period the total number of courses
must be at least some minimum number $c$ and at most a maximum
number $d$. In addition, within any given period the total academic
load of the courses must be at least some minimum load $a$ and
at most a maximum load $b$.  Let $R$ be the set of prequisite pairs
$\n{i,j}$, where $i$ and $j$ are courses. Prerequisite relationships
must be observed, so for each pair of courses $\n{i,j} \in R$, if
course $i$ is scheduled in period $p$, then course $j$ must be
scheduled in a period strictly prior to $p$.


A model of this problem using set variables and set constraints was
proposed by \citeA{hnichkw:bacp}, although no experimental results
were presented for that model.  In this ``primal'' model, we use one set
variable
$S_i$ per academic period. Each $S_i$ represents the set of courses
assigned to academic period $i$, so $k \in S_i$ if and only if course
$k$ is assigned to period $i$. We can model the problem using the
following constraints:

\begin{itemize}
\item (S1) Every course is taken exactly once: 
    $\partition(S_1, \dotsc, S_n)$
\item (S2) The number of courses in a period is between $c$ and
$d$:\\
    $\quad\forall_{i=1}^n ((c \leq |S_i|) \wedge (|S_i| \leq d))$
\item (S3) The total academic load in a period is between $a$ and
$b$: \\
    $\quad\forall_{i=1}^n ((a \leq \wsum(S_i, \n{t_1, \dotsc,
    t_m})) \wedge (\wsum(S_i,\n{t_1, \dotsc, t_m}) \leq b))$
\item (S4) Prerequisites are respected:\\
    $\quad\forall_{i=1}^{n} \forall_{j=1}^i (\bigwedge_{\n{c,p} \in
    R} (p \in S_i) \ra (c \notin S_j))$
\end{itemize}



In general the $\partition$ and $\wsum$ constraints of the primal set
model can be very large, making domain propagation impractical.

In addition to presenting results for the primal model, 
\citeA{choils:propredundancy} demonstrated that a substantial
performance improvement can be obtained for this problem through the
use of redundant models together with channelling constraints.

We can obtain better results from a ``dual'' set model, where additional
set variables are introduced to model each course. We define set
variables $X_i$ ($1 \leq i \leq m$) representing each course, such
that if $k \in X_i$ then course $i$ is assigned to period $k$. We can
then define the following constraints:
\begin{itemize}
\item (CX) Channelling constraints:
    $\forall_{i=1}^m \forall_{j=1}^n (i \in S_j) \lra (j \in X_i)$
\item (X1) Each course may be assigned to at most one period:
    $\forall_{i=1}^m |X_i| = 1$
\item (X2) Prerequisites are respected: $\forall \n{i,j} \in R ~
    \lexlt(X_j, X_i)$
\end{itemize}

We define the dual model to consist of the constraints
\{S2, S3, CX, X1, X2\}. Constraints S1 and S4 are no longer required, since
they are propagation redundant.

Unfortunately, while the dual model performs better than the primal
set model, it is still incapable of proving the optimality of the
solutions to any of the problems. We can obtain stronger propagation
by introducing redundant global constraints on the total academic load and
number of courses of all academic periods, as suggested 
by~\citeA{choils:propredundancy}.

For each academic period $j$ define two integer variables $l_j$,
representing the academic load in period $j$ and $q_j$, representing
the number of courses in period $j$. We can then define the following
constraints on these new integer variables:
\begin{itemize}
\item (CI1) Channeling constraints on $l_i$:
    $\forall_{i=1}^n \wsum(S_i, \n{t_1, \dotsc, t_m}) = l_i$
\item (CI2) Channeling constraints on $q_i$:
    $\forall_{i=1}^n |S_i| = q_i$
\item (I1) Range constraints on $l_i$:
    $\forall_{i=1}^n (a \leq l_i) \wedge (l_i \leq b)$
\item (I2) Range constraints on $q_i$:
    $\forall_{i=1}^n (c \leq q_i) \wedge (q_i \leq d)$
\item (I3) All loads should be undertaken:
    $\sum_{i=1}^n l_j = \sum_{i=1}^n t_i$
\item (I4) All courses must be taken:
    $\sum_{i=1}^n q_j = m$
\end{itemize}

We then define the ``hybrid dual'' model consisting of the constraints
\{CX, X1, X2, CI1, CI2, I1, I2, I3, I4\}. Note that none of the
original constraints S1--S4 remain. For completeness it is also
possible to define a ``hybrid primal'' model consisting of constraints
\{S1, S4, CI1, CI2, I1, I2, I3, I4\}, although as before the large
ROBDD sizes make domain propagation impractical.

\begin{table}
\caption{Performance results for the Balanced Academic Curriculum
Problem for the 8, 10 and 12 period problems. The first column
contains the number of periods, the second column contains the model
type (HD=hybrid dual, D=dual, P=primal, HP=hybrid primal), 
and the third column contains the maximum load per period $b$. In all cases a
sequential ``smallest-element-in-set'' labeling method was used
}
\label{tab:bcur-results}
\small
\begin{center}
\begin{tabular}{||l|l|l||r|r||r|r||r||r|r||r|r||}
\hhline{|t:===:t:==:t:==:t:=:t:==:t:==:t|}
\multicolumn{3}{||l||}{} &
\multicolumn{2}{c||}{Bounds} &
\multicolumn{2}{c||}{Domain} &
\multicolumn{1}{c||}{LU+R} &
\multicolumn{2}{c||}{LU+Lex} &
\multicolumn{2}{c||}{LU+Card} 
\\
\hhline{||~~~||--||--||-||--||--||}
\multicolumn{3}{||l||}{Problem} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c||}{time} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} &
\multicolumn{1}{c|}{time} &
\multicolumn{1}{c||}{fails} 
\\[-3pt]
\multicolumn{3}{||l||}{} &
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c||}{\scriptsize /s} &
\multicolumn{1}{c|}{\scriptsize /s} && 
\multicolumn{1}{c|}{\scriptsize /s} &
\\
\hhline{||---||--||--||-||--||--||}
\multirow{12}{16pt}{8} & \multirow{3}{16pt}{HD} & 16 & --- & --- &\bf  0.3 &\bf  0 &\bf  0.3 & \bf 0.3 &\bf  0 & --- & --- \\
& & 17 & 8.1 & 5586 & 1.7 & \bf 3 & 2.7 &\bf  1.0 &\bf  3 & 2.3
&\bf  3  \\
& & 18 & 35.7 & 30430 & 1.7 &\bf  5 & 3.0 &\bf  1.0 & \bf 5 & 2.3
& \bf 5  \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{D} & 16 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 17 & 7.7 & 5348 & 1.1 & \bf 3 & 1.2 & \bf 0.6 & \bf 3 & 1.7 & \bf 3  \\
& & 18 & 32.8 & 29781 & 1.2 & \bf 5 & 1.3 & \bf 0.7 & \bf 5 & 1.6
& \bf 5  \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{P} & 16 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 17 & \bf 6.0 & 5353 & --- & --- & --- & 13.9 & \bf 2046 & 32.7 &
4757 \\
& & 18 & \bf 16.3 & 29781 & --- & --- & --- & 53.1 & \bf 26431 & 61.2 &
29628 \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{HP} & 16 & --- & --- & --- & --- & --- & \bf 0.3 & \bf 0 & --- & --- \\
& & 17 & \bf 6.8 & 5586 & --- & --- & --- & 12.8 & \bf 1510 & 40.5 &
5410 \\
& & 18 & \bf 19.3 & 30431 & --- & --- & --- & 69.8 & \bf 26145 & 70.5 &
30141 \\
\hhline{||---||--||--||-||--||--||}
\multirow{12}{16pt}{10} & \multirow{3}{16pt}{HD} & 13 & --- & --- & \bf 0.3 & \bf 0 & \bf 0.3 & \bf 0.3 & \bf 0 & --- & --- \\
& & 14 & 17.2 & 11439 & 1.3 & \bf 2 & 1.5 & \bf 0.9 & \bf 2 & 2.6
& \bf 2 \\
& & 15 & 2.1 & 630 & 1.4 & \bf 1 & 1.7 & \bf 0.9 & \bf 1 & 2.6 &
\bf 1 \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{D} & 13 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 14 & 26.8 & 21103 & 1.0 & \bf 2 & 1.6 & \bf 0.6 & \bf 2 &
2.4 & \bf 2 \\
& & 15 & 2.1 & 711 & 1.0 & \bf 1 & 1.6 & \bf 0.6 & \bf 1 & 2.1 &
\bf 1 \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{P} & 13 & --- & --- & --- & --- & --- & --- & --- & --- & ---  \\
& & 14 & \bf 14.7 & 21105 & --- & --- & --- & 53.3 & 16275 & 52.1
& \bf 11772 \\
& & 15 & \bf 3.5 & 711 & --- & --- & --- & 9.6 & 593 & 31.5 & \bf 515 \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{HP} & 13 & --- & --- & --- & --- & --- & \bf 0.4 & \bf 0 & --- & ---  \\
& & 14 & \bf 11.8 & 11440 & --- & --- & --- & 25.6 & \bf 4672 & 46.9 &
5773 \\
& & 15 & \bf 3.6 & 630 & --- & --- & --- & 10.6 & \bf 516 & 31.7 & 533 \\
\hhline{||---||--||--||-||--||--||}
\multirow{12}{16pt}{12} & \multirow{3}{16pt}{HD} & 16 & --- & --- & 1.1 & \bf 0 & 1.1 & \bf 0.9 & \bf 0 & --- & --- \\
& & 17 & --- & --- & 6.0 & \bf 9 & 9.3 & \bf 3.0 & \bf 9 & --- & --- \\
& & 18 & 6.1 & 235 & 8.1 & \bf 4 & 9.2 & \bf 3.4 & \bf 4 & 15.7 &
\bf 4 \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{D} & 16 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 17 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 18 & 6.3 & 220 & 6.8 & \bf 4 & 9.0 & \bf 3.0 & \bf 4 & 13.1 &
\bf 4 \\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{P} & 16 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 17 & --- & --- & --- & --- & --- & --- & --- & --- & --- \\
& & 18 & \bf 101.7 & 225 & --- & --- & --- & 153.8 & \bf 224 & --- & ---
\\
\hhline{||~--||--||--||-||--||--||}
& \multirow{3}{16pt}{HP} & 16 & --- & --- & --- & --- & --- & \bf 3.4 & \bf 0 & --- & --- \\
& & 17 & --- & --- & --- & --- & --- & \bf 216.8 & \bf 526 & --- & --- \\
& & 18 & \bf 100.9 & 236 & --- & --- & --- & 166.0 & \bf 224 & --- & ---
\\

\hhline{|b:===:b:==:b:==:b:=:b:==:b:==:b|}
\end{tabular}
\end{center}
\end{table}

Experimental results for the Balanced Academic Curriculum problem are
shown in Table~\ref{tab:bcur-results}. The timing results for the Hybrid
Dual model are comparable to the best results obtained by
\citeA{hnichkw:bacp} and \citeA{choils:propredundancy}. The number of
failures required to solve the problem is orders of magnitude smaller
than the best results presented in either paper, emphasising the value
of domain propagation in this case.  These examples also show a case
where the LU+Lex solver is competitive and clearly the most robust
solver over the different models.

%\pjs{multiset example}


\ignore{
% This is too similar to the Balanced Academic Curriculum problem, and
% so it probably isn't worth including it in the journal article.
\subsection{Progressive Parties}
\label{sect:exp-pparty}

The Progressive Party problem (problem \texttt{prob013} of CSPLib) has
the somewhat frivolous goal of planning a party at a yacht club. At
the yacht club there are $n$ boats, each with a particular crew size
$R_i$ and maximum total capacity $C_i$ ($1 \leq i \leq n$). $h$ of
these are designated as host boats and the remaining $g = n-h$ boats 
as guest boats.
The crews of the host boats remain on board their own boats, while the
crews of the guest boats progressively visit the host boats over $t$
time periods, each guest boat visiting exactly one host boat per time period.
The crew of any given boat must remain
together at all times, and the maximum capacity of any boat must not
be exceeded. In addition, in order to ensure that the participants
meet as many new people as possible, any two boat crews (either hosts
or guests) are permitted to meet no more than once.
The goal of the problem is to
devise a schedule of allocations of guest boats to host boats for each
time period, such that all of these constraints are satisfied.

Smith \ph{cite} investigated the performance of both CP and Integer
Linear Programming (ILP) solvers for this problem. In their model of the
problem the first 13 boats were arbitrarily chosen to be host boats.

We can model this problem as a set problem using a $t \times n$ matrix
of set variables $H_{ij}$ ($1 \leq i \leq t$, $1 \leq j \leq n$),
drawn from the universe $\mathcal{U} = \{1,\dotsc,g\}$. Each set
variable $H_{ij}$ represents the set of guest boats visiting host boat
$j$ at time period $i$. We impose the following constraints:
\begin{itemize}
\item (S1) Capacity constraints:
    $\forall_{i=1}^t \forall_{j=1}^h \wsum(H_{ij}, \n{R_{h+1},
    R_{h+2}, \dotsc, R_n}) \leq C_j - R_j$
\item (S2) Partitioning constraints:
    $\forall_{i=1}^t \partition(H_{i1}, \dotsc, H_{ih})$
\item (S3) Guests may not revisit a host:
    $\forall_{i=1}^{t-1} \forall_{j=i+1}^{t} \forall_{k=1}^h H_{ik}
    \cap H_{jk} = \phi$
\item (S4) Guests may not meet more than once: \\
    \[\forall_{i=1}^{t-1} \forall_{j=i+1}^{t} \forall_{a=1}^h
    \forall_{b \in \{1,\dotsc,h\}  \setminus \{a\}} \, |H_{ia} \cap
    H_{jb}| \leq 1\]
\item (S5) Lexicographic ordering on the first boat:
    $\forall_{i=1}^{t-1} \forall_{j=i+1}^{t} \lexlt(H_{i1}, H_{j1})$
\item (S6) Lexicographic ordering on boats of equal capacity in the
    first time period:\\
    \[\forall_{i=1}^{h-1} \forall_{j=i+1}^{h} ((C_i - R_i) = (C_j -
    R_j)) \ra \lexlt(H_{1i}, H_{1j}) \]
    
\end{itemize}

As before, the ROBDDs representing the $\partition$ and $\wsum$ constraints
can be very large,
so it is impractical to perform domain propagation directly on this
model. We can obtain better performance through the use of
redundant modeling. We define an additional $t \times g$ matrix of 
set variables
$G_{ij}$ ($1 \leq i \leq t$, $1 \leq j \leq g$).

}

\subsection{Comparing the Propagation Performance of the Solvers}

\begin{figure}
\begin{tabular}{ll}
Steiner system $S(2,3,15)$ \\
\includegraphics[height=6.3cm]{st15-domsize-nlabels} &
\includegraphics[height=6.3cm]{st15-bddsize-nlabels}
\\
\\
Social Golfers problem 4-5-4 \\
\includegraphics[height=6.3cm]{gf454-domsize-nlabels} &
\includegraphics[height=6.3cm]{gf454-bddsize-nlabels}
\end{tabular}

\caption{Comparison of total domain and total ROBDD sizes with labeling
step for set bounds, set domain, split set domain, set and
lexicographic bounds, and set and cardinality bounds solvers on the
Steiner System $S(2,3,15)$ and Social Golfers problem 4-5-4. Note that
the $y$ axes of the total domain size graphs have a logarithmic scale
(base 2).}

\label{fig:size-vs-labeling}
\rule{\textwidth}{0.4pt}
\end{figure}

\begin{figure}
\begin{tabular}{cc}
\includegraphics[height=6.3cm]{st15-domsize-time} &
\includegraphics[height=6.3cm]{gf454-domsize-time}
\\
\\
Steiner system $S(2,3,15)$ & Social Golfers problem 4-5-4
\end{tabular}
\caption{Comparison of total domain size with time for set
bounds, set domain, split set domain, set and lexicographic bounds,
and set and
cardinality bounds solvers on the Steiner System $S(2,3,15)$ and
Social Golfers problem 4-5-4. Note
that the $y$ axes have a logarithmic scale (base 2).}
\label{fig:size-vs-time}
\rule{\textwidth}{0.4pt}
\end{figure}

It is instructive to compare the propagation performance 
of the various ROBDD-based solvers graphically.
Here we pick two small examples, namely the Steiner System $S(2,3,15)$ 
and the Social Golfers problem 4-5-4,
using the default labeling for each problem, and graph BDD and domain
sizes against
number of labeling steps and time.


Figure~\ref{fig:size-vs-labeling} depicts how the total domain
sizes and total ROBDD sizes (total number of internal nodes) 
vary with each labeling step for each of
the solvers for the Steiner System $S(2,3,15)$ and for the Social
Golfers problem 4-5-4. In particular, observe
that the domain and split-domain propagators are the
most effective at reducing the domain size and thus restricting the
search space, although
maintaining domain consistency can be costly due to the size of the
ROBDDs required for representing arbitrary domains. This can
be seen in the ROBDD size against labeling step graph, most notably in
the case of the Social Golfers problem 4-5-4.
The lexicographic bounds solver is the next most effective
in domain size reduction, and requires a smaller number of ROBDD nodes
to store its bounds than the domain or cardinality bounds solvers.
The bounds solver is clearly the weakest in terms of domain size
reduction, and as opposed to the other solvers starts from a
zero size domain representation and builds slowly to ROBDDs
representing the answer.

We are not only concerned with the strength of a propagator in the
abstract, since the efficiency of the solver as a whole is dependent on
both the propagation and labeling processes.
Figure~\ref{fig:size-vs-time} depicts graphs of domain size against
time for each of the propagators discussed in this paper. Here we see
that the most effective propagator does not necessarily lead to the
most efficient set solver. Despite comparatively weak propagation, the
set and cardinality bounds solvers lead to the fastest reduction in
domain size per unit time. Nonetheless, the experimental timing
results given above demonstrate that for harder cases the additional
cost of maintaining domain consistency or lexicographic bounds
can be justified through the
consequent reduction in search space.

Interestingly we can see from the graphs in Figure~\ref{fig:size-vs-time}
that the initial domain reduction before
labeling, which is the gap from time 0 to where the line starts, 
can be a significant part of the computation.
In particular, one of the weaknesses of the lexicographic bounds solver is
the large time required to reach an initial fixpoint through lexicographic
bounds propagation. For practical applications it might be preferable to
implement a hybrid solver which utilises one of the other propagators
to generate the initial domains, and uses the lexicographic bounds
propagator during labeling.

\section{Conclusion}
\label{sect:conclusion}

We have demonstrated that ROBDDs form a highly flexible platform for
constructing set constraint solvers. ROBDDs allow a compact
representation of many set domains and set constraints, making them an
effective basis for an efficient set solver. Since ROBDDs can
represent arbitrary Boolean formul\ae{} we can easily conjoin and
existentially quantify them, permitting the removal of intermediate
variables and making the construction of global constraints trivial.
We have demonstrated how to efficiently enforce various levels of consistency,
including set domain, set bounds, cardinality bounds
and lexicographic bounds consistency.
Finally, we have presented experimental results that demonstrate that
the ROBDD-based solver outperforms other common set solvers on a wide
variety of standard set constraint satisfaction problems.

No single set solver is uniformly better than the others.
For many examples simple bounds propagation is the best 
approach, while in other cases, particularly when we ask for all
solutions, domain consistency is preferable.
There are also examples where lexicographic bounds or cardinality bounds
are the best approach.  The split-domain approach is somewhat disappointing,
since it appears that in many cases the overhead 
of the more complicated calculation (Equation~\ref{eq:split}) 
is not rewarded in terms of smaller ROBDD sizes.

In the future we plan to explore a robust general set constraint solver 
that dynamically chooses which level of consistency to maintain
by examining how big the domain ROBDDs are becoming as search
progresses. 
\vskip 0.2in
\bibliography{phlocal}
\bibliographystyle{theapa}
\end{document}

