SMT-for-IN3070/committees.org

#+TITLE: SMT for IN3070
#+AUTHOR: Lars Tveito
#+HTML_HEAD: <script type="text/javascript" src="js/script.js"></script>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="Rethink/rethink.css" />
#+OPTIONS: toc:nil num:nil html-style:nil

At the Department of Informatics (University of Oslo), all exams are corrected
by a committee consisting of two examiners. For large courses, there are often
many examiners where some want to correct more than others. The administration
is responsible for forming these committees. Sometimes there are additional
constraints on which examiners can and cannot form a committee, for example,
due to different levels of experience.

Before digitizing exams at the department, the administration would have
physical copies of the exam to distribute. This would actually make it easier
to form the committees because the constraints could be handled on the fly.
When digitized, the problem would essentially turn into a math problem which in
the general case is not particularly easy to solve.

This is an actual email (in Norwegian) forwarded to me from someone in the
administration:

#+BEGIN_QUOTE
Mine mattekunnskaper er tydeligvis fraværende. Jeg klarer ikke finne en
fornuftig løsning på dette:

| A | 160 |
| B | 150 |
| C | 110 |
| D |  60 |
| E |  60 |
| F |  30 |

Det er snakk om sensur i inf1300 med 283 besvarelser. D og E kan ikke rette
mot hverandre. De bør helst rette mot B eller C.

Har du et bra forslag til meg? Jeg blir GAL. Det var bedre før, da jeg hadde
besvarelsene fysisk å kunne telle ark.

Har du mulighet til å hjelpe en stakkar?

Takk
#+END_QUOTE

We want to answer this cry for help with a general solution for this problem
using SMT-solving.

* Satisfiability modulo theories (SMT)

  SMT-solvers are tools for solving satisfiability problems, i.e. given a
  first- order logical formula $\phi$, decide whether or not there exists a
  model $\mathcal{M}$ such that $\mathcal{M} \models \phi$. In general, this is
  an undecidable problem. However, there are theories within first-order logic
  that are decidable. SMT solvers can produce models that satisfy a set of
  formulas for many useful theories, some of which are decidable. It is natural
  to think of SMT as a generalization of SAT, which is satisfiability for
  propositional logic.

  The solver we will be using is [[https://github.com/Z3Prover/z3][Z3]].

** Theories

   Examples of theories can be the theory of booleans (or propositional logic),
   integers or real numbers with equations and inequations, or other common
   programming concepts like arrays or bitvectors. Z3 supports solving
   constraint problems in such theories. More formally, we define theories as
   follows:

   #+BEGIN_definition
   A theory is a set of first-order logic formulas, closed under logical
   consequence.
   #+END_definition

   We can imagine how this might work. The natural numbers can, for instance,
   be axiomatized with the Peano axioms.

   1. $\forall x \in \mathbb{N} \ (0 \neq  S ( x ))$
   2. $\forall x, y \in \mathbb{N} \ (S( x ) =  S( y ) \Rightarrow x = y)$
   3. $\forall x \in \mathbb{N} \ (x  + 0 = x )$
   4. $\forall x, y \in \mathbb{N} \ (x + S( y ) =  S( x + y ))$
   5. $\forall x \in \mathbb{N} \ (x \cdot 0 = 0)$
   6. $\forall x, y \in \mathbb{N} \ (x \cdot  S ( y ) = x \cdot y + x )$

   In addition, one axiom is added to formalize induction. Because a theory is
   closed under logical consequence, the theory consists of all true
   first-order sentences that follow from these axioms, which correspond to the
   true sentences about natural numbers.

   However, in Z3, we don't see such axioms, but axiomatizations provide the
   formal justification for implementing special solvers for commonly used
   theories. In Z3, we could write something like this:

   #+BEGIN_SRC z3
   (declare-const a Int)
   (declare-const b Int)
   (declare-const c Int)

   (assert (< 0 a b c))

   (assert (= (+ (* a a) (* b b)) (* c c)))

   (check-sat)
   (get-model)
   #+END_SRC

   This encodes two constraints
   - $0 < a < b < c$
   - $a^2 + b^2 = c^2$
   where $a,b,c$ are whole numbers. Then we ask Z3 to produce a model
   $\mathcal{M}$ such that $\mathcal{M} \models (0 < a < b < c) \land (a^2 +
   b^2 = c^2)$, which outputs:

   #+BEGIN_EXAMPLE
   sat
   (model
     (define-fun c () Int
       5)
     (define-fun b () Int
       4)
     (define-fun a () Int
       3)
   )
   #+END_EXAMPLE

   The first line ~sat~ indicates that the formula is satisfiable, and produce
   a model where $a^\mathcal{M}=3$, $b^\mathcal{M}=4$ and $c^\mathcal{M}=5$.

   Note that we would get a different answer if we declared the constant
   symbols as real numbers because Z3 would use the theory for reals to satisfy
   the constraints.

** Many-sorted first-order logic

   Z3 implements [[http://smtlib.cs.uiowa.edu/papers/smt-lib-reference-v2.6-r2017-07-18.pdf][SMT-LIB]], a standardized syntax and semantics for SMT solvers.
   Its underlying logic is a /Many-sorted first-order logic/, where values must
   have an associated sort (which is a basic form of type). Think of it as
   partitioning the domain, where each sort corresponds to a part. A signature
   in a many-sorted first logic is defined as follows.

   #+BEGIN_definition
   A signature $\Sigma = (S, F, P)$ consists of a countable set of
   - Sorts $S$.
   - Function symbols $F$, where each member is a function symbol $f$ with an
     associated type $s_1 \times \dots \times s_n \to s$, where $s \in S$ and
     $s_1, \dots, s_n \in S$. Constants are simply zero-arity function symbols.
   - Predicate symbols $P$, where each predicate has an associated type $s_1
     \times \dots \times s_n$. We assume an equality $=_s$ predicate with type
     $s \times s$ for all sorts in $S$.
   #+END_definition

   The equality relation will be denoted $=$, letting the sort remain implicit.

   For example, the signature for the integers can be formalized as
   $\Sigma_{int} = (S_{Int}, F_{Int}, P_{Int})$ where
   - $S_{Int} = \{Int\}$
   - $F_{Int} = \{0, 1, +, -, *\}$ where the constant symbols $0, 1$ has a type
     $Int$ and the function symbols $+,-,*$ has a type signature $Int \times
     Int \to Int$.
   - $P_{Int} = \{<, =\}$ where the predicate symbols $<, =$ has type signature
     $Int \times Int$.

   In Z3, the type signature of function- and predicate symbols informs Z3 of
   what theory it should apply.

* Back to the problem

  We have 283 exams. Every exam must be corrected by a committee consisting of
  two examiners. Each examiner has an associated capacity of exams they want to
  correct. Examiners D and E can't be in the same committee, and should rather
  be in committee with B or C. We prefer a smaller number of committees.

  We use the [[https://ericpony.github.io/z3py-tutorial/guide-examples.htm][Python API for Z3]]. Install with:

  #+BEGIN_SRC sh
   pip3 install z3-solver
  #+END_SRC

  Create a Python file and populate it with:

  #+BEGIN_SRC python :tangle committees.py
  from z3 import *
  #+END_SRC

  This allows us to generate instances with Python that Z3 can solve.

** Instances

   Let's formulate an instance as a four-tuple $(N, C, S, A)$ where
   - $N$ is the number of exams to correct
   - $C$ is a list of capacities, where each examiner is identified by
     their position of the list
   - $S$ is a mapping from a single examiner to a set of examiners they /should/ form a committee with
   - $A$ is a symmetric relation, relating examiners that we should /avoid/
     placing in the same committee

   We define a committee as a set of exactly two examiners (identified by their
   index in the list of capacities).

   The code below suggests a Python representation of a problem instance. It
   is, as you must have noticed, blurred (until you click it). This is to
   encourage the reader to solve the problem on their own and emphasize that
   what will be presented is a mere suggestion on how to attack the problem.

   #+BEGIN_SRC python :tangle committees.py
   def example_instance():
       N = 283
       #    A    B    C    D   E   F
       C = [160, 150, 110, 60, 60, 30]
       S = {3 : {1, 2}, 4 : {1, 2}}
       A = {frozenset([3, 4])}
       return (N, C, S, A)
   #+END_SRC

** Constraint modeling

   We need to capture our intention with first-order logic formulas, and
   preferably quantifier-free. In the context of SMT-solving, quantifier-free
   means that we only try to solve a set of constraints where no variable is
   bound by a quantifier; these are usually much easier to solve. Rather, we
   use a finite set of constant symbols, with some associated sort, and try to
   find an interpretation for them.

   The end result needs to be a set of committees, where each committee
   consists of two examiners with a number of exams to correct. An important
   part of finding a reasonable encoding is to balance what part of the problem
   should be solved with Python and what should be solved by the SMT-solver. My
   experience is that a good rule of thumb is to move as much structural
   complexity to Python and encode the Z3 instance with simple structures.

** Modeling committees

   A natural encoding could be modeling a committee as an integer constant,
   where the value assigned to a committee corresponds to the number of exams
   they correct. If the committee is not assigned any exams, we discard it
   completely. It is quite easy to compute all possible committees and make one
   integer constant for each of them.

   Let's write a function that takes a list of capacities, and return a
   dictionary, associating committees to their corresponding integer constant.
   Remember that we represent a committee as a set of exactly two examiners.

   #+BEGIN_SRC python :tangle committees.py
   def committees(C):
       cs = {frozenset([i,j])
             for i in range(len(C))
             for j in range(i+1, len(C))}
       return {c : Int(str(c)) for c in cs}
   #+END_SRC

** Capacity constraints

   Now we must ensure that no examiner receives more exams than their capacity.
   Given an examiner $i$, where $0 <= i < |C|$, we let $c_i$ denote the set of
   all committees $i$ participates in. Then $\sum{c_i} <= C[i]$, i.e. the sum
   of the exams corrected by committees in $c_i$ does not exceed the capacity
   of the examiner $i$. We write a function that encodes these constraints.

   #+BEGIN_SRC python :tangle committees.py
   def capacity_constraint(comms, C):
       return [sum(comms[c] for c in comms if i in c) <= C[i]
               for i in range(len(C))]
   #+END_SRC

   Because we are modeling committees as integers, we have to be careful not to
   allow committees correcting a negative number of exams.

   #+BEGIN_SRC python :tangle committees.py
   def non_negative_constraint(comms):
       return [0 <= comms[c] for c in comms]
   #+END_SRC

** Committee constraints

   The $S$ relation is sort of odd. That one examiner /should/ form a committee
   with someone they relate to by $S$. This is not an absolute requirement,
   which is not ideal for a satisfiability problem, so we will ignore this
   constraint for now. The $A$ relation is similar but clearer. For any pair
   $(i,j) \in A$, we don't form a committee consisting of those examiners.

   #+BEGIN_SRC python :tangle committees.py
   def avoid_correct_with_constraint(comms, A):
       return [comms[frozenset([i, j])] == 0 for i, j in A]
   #+END_SRC

** All exams are corrected constraint

   Each committee corrects their exams two times (once by each examiner), so if
   the sum of all the committees is $N$, then all exams have been corrected
   twice (presumably by two different examiners). Let's encode that as a
   constraint.

   #+BEGIN_SRC python :tangle committees.py
   def all_corrected_constraint(comms, N):
       return sum(comms.values()) == N
   #+END_SRC

** Invoking Z3

   Now that we have functions that model our problem, we can invoke Z3.

   #+BEGIN_SRC python :tangle committees.py
   def check_instance(instance):
       N, C, S, A = instance
       comms = committees(C)

       s = Solver()

       s.add(capacity_constraint(comms, C))
       s.add(non_negative_constraint(comms))
       s.add(all_corrected_constraint(comms, N))
       s.add(avoid_correct_with_constraint(comms, A))

       s.check()
       return s.model()
   #+END_SRC

   Calling ~check_instance(example_instance())~ returns a model:

   #+BEGIN_EXAMPLE
   [frozenset({2, 4}) = 0,
    frozenset({0, 2}) = 0,
    frozenset({2, 3}) = 0,
    frozenset({1, 3}) = 0,
    frozenset({2, 5}) = 0,
    frozenset({3, 5}) = 0,
    frozenset({0, 5}) = 13,
    frozenset({1, 2}) = 110,
    frozenset({4, 5}) = 0,
    frozenset({1, 5}) = 17,
    frozenset({0, 3}) = 60,
    frozenset({0, 4}) = 60,
    frozenset({0, 1}) = 23,
    frozenset({3, 4}) = 0,
    frozenset({1, 4}) = 0]
   #+END_EXAMPLE

   This is not especially readable, so let's write a quick (and completely
   unreadable) prettyprinter.

   #+BEGIN_SRC python :tangle committees.py
   def prettyprint(instance, m):
       N, C, S, A = instance
       comms = committees(C)
       exams = [sum(m[comms[c]].as_long() for c in comms if i in c)
                for i in range(len(C))]
       examiners = '\n'.join(['%s: %d/%d' % (chr(ord('A') + i), exams[i], C[i])
                              for i in range(len(C))])
       cs = [(c, m[comms[c]].as_long()) for c in sorted(comms, key=sorted)]
       csstr = '\n'.join([', '.join(map(lambda i: chr(ord('A') + i),
                                        sorted(c))) + ': ' + str(cv)
                          for c, cv in cs if cv > 0])
       print(examiners + '\n\n' + csstr)
   #+END_SRC

   This outputs the something like:

   #+BEGIN_EXAMPLE
   A: 156/160
   B: 150/150
   C: 110/110
   D: 60/60
   E: 60/60
   F: 30/30

   A, B: 23
   A, D: 60
   A, E: 60
   A, F: 13
   B, C: 110
   B, F: 17
   #+END_EXAMPLE

   Note the /something like/. There are multiple ways to satisfy this set of
   constraints, and Z3 only provide /some/ solution (if one exists).

* Optimization

  So far, we have found a way to model the problem and satisfy the constraints.
  However, it is preferable to have fewer committees, because all committees
  have to discuss the exams, causing administrative overhead. Z3 also provides
  optimization, meaning that we can find a smallest or largest solution for
  numeric theories. The underlying theory for optimization is MaxSMT.

** Minimize committees

   In our case, we want to minimize the number of committees. First, we write a
   function to find the number of committees which we will soon minimize.

   #+BEGIN_SRC python :tangle committees.py
   def number_of_committees(comms):
       return sum(If(0 < comms[c], 1, 0) for c in comms)
   #+END_SRC

   Now we can invoke Z3, using an ~Optimize~ instance and adding our
   minimization constraint.

   #+BEGIN_SRC python :tangle committees.py
   def optimize_instance(instance):
       N, C, S, A = instance
       comms = committees(C)

       o = Optimize()

       o.add(capacity_constraint(comms, C))
       o.add(non_negative_constraint(comms))
       o.add(all_corrected_constraint(comms, N))
       o.add(avoid_correct_with_constraint(comms, A))

       o.minimize(number_of_committees(comms))

       o.check()
       return o.model()
   #+END_SRC

   There is still more than one way to satisfy this model, but we are
   guaranteed to get a minimal number of committees (which is 6 in our
   example).

   #+BEGIN_EXAMPLE
   A: 160/160
   B: 150/150
   C: 110/110
   D: 56/60
   E: 60/60
   F: 30/30

   A, B: 57
   A, D: 43
   A, E: 60
   B, C: 93
   C, F: 17
   D, F: 13
   #+END_EXAMPLE

** Dealing with /should/

   Remember $S$, which maps examiners to other examiners they /should/ form a
   committee with. With optimization, we now have a way of expressing that some
   solution is more preferable than another. One way to model this is
   maximizing the number of exams given to committees that consists of an
   examiner $i$ that should be in a committee with examiner $j$. We want this
   for all such pairs $i,j$, and can achieve this by summing all such
   committees.

   #+BEGIN_SRC python :tangle committees.py
   def should_correct_with_weight(comms, S, C):
       return sum(comms[frozenset([i, j])] for i in S for j in S[i])
   #+END_SRC

   When adding multiple optimization objectives (or goals), Z3 defaults to
   order the objectives lexicographically, i.e. in the order they appear. If we
   place the minimization of committees before the ~should_correct_with_weight~, then we are still guaranteed to get a minimal
   number of committees.

   #+BEGIN_SRC python :tangle committees.py
   def optimize_instance(instance):
       N, C, S, A = instance
       comms = committees(C)

       o = Optimize()

       o.add(capacity_constraint(comms, C))
       o.add(non_negative_constraint(comms))
       o.add(all_corrected_constraint(comms, N))
       o.add(avoid_correct_with_constraint(comms, A))

       o.minimize(number_of_committees(comms))
       o.maximize(should_correct_with_weight(comms, S, C))

       o.check()
       return o.model()
   #+END_SRC

   #+BEGIN_EXAMPLE
   A: 156/160
   B: 150/150
   C: 110/110
   D: 60/60
   E: 60/60
   F: 30/30

   A, B: 90
   A, C: 43
   A, F: 23
   B, E: 60
   C, D: 60
   C, F: 7
   #+END_EXAMPLE

** Optimize for capacities

   Maybe we can try to satisfy (🙃) all the examiners by trying to close the
   gap between their capacity and the number of exams they end up correcting.
   Usually, there is quite a lot of flex in these capacities; if you are
   willing to correct $50$ exams, then you will most likely be okay with
   correcting $40$ and /actually/ willing to correct $52$. Therefore, we can
   try to add some slack to the capacity.

   In reality, the numbers from the original email were

   | A | 158 |
   | B | 150 |
   | C | 108 |
   | D |  60 |
   | E |  60 |
   | F |  15 |

   But when we add them up, it turns out that they only have capacity to
   correct $551$ exams (and we need $2*N = 566$).

   We create a new instance with the original values.

   #+BEGIN_SRC python :tangle committees.py
   def original_instance():
       N = 283
       #    A    B    C    D   E   F
       C = [158, 150, 108, 60, 60, 15]
       S = {3 : {1, 2}, 4 : {1, 2}}
       A = {frozenset([3, 4])}
       return (N, C, S, A)
   #+END_SRC

   Now we can compute a "badness"-score (or weight) for the examiners
   capacities, rather than just stating we cannot surpass their capacity.

   #+BEGIN_SRC python :tangle committees.py
   def capacity_weight(comms, i, C):
       a = sum(comms[c] for c in comms if i in c)
       return If(a > C[i], a - C[i], C[i] - a)
   #+END_SRC

   For the total weight of the capacities, we try to just sum the weights for
   each examiner.

   #+BEGIN_SRC python :tangle committees.py
   def capacity_weight(comms, C):
       return sum(capacity_weight(comms, i, C) for i in range(len(C)))
   #+END_SRC

   We can now add all of the optimization objectives, stating that it most
   important to respect the capacities of the examiners, then prefer a small
   number of committees, and lastly the /should/ requirement from the previous
   section.

   #+BEGIN_SRC python :tangle committees.py
   def optimize_instance(instance):
       N, C, S, A = instance
       comms = committees(C)

       o = Optimize()

       o.add(non_negative_constraint(comms))
       o.add(all_corrected_constraint(comms, N))
       o.add(avoid_correct_with_constraint(comms, A))

       o.minimize(capacity_weight(comms, C))
       o.minimize(number_of_committees(comms))
       o.maximize(should_correct_with_weight(comms, S, C))

       o.check()
       return o.model()
   #+END_SRC

   We now get something like:

   #+BEGIN_EXAMPLE
   A: 158/160
   B: 158/150
   C: 110/110
   D: 65/60
   E: 60/60
   F: 15/30

   A, B: 158
   C, D: 65
   C, E: 45
   E, F: 15
   #+END_EXAMPLE

   If we were to prioritize the /should/ requirement over minimizing the number
   of committees, then we would get something like:

   #+BEGIN_EXAMPLE
   A: 158/160
   B: 158/150
   C: 109/110
   D: 65/60
   E: 60/60
   F: 16/30

   A, B: 98
   A, C: 44
   A, F: 16
   B, E: 60
   C, D: 65
   #+END_EXAMPLE

   At this point, I hope you have realized that we now have a tool that we can
   use to derive a very flexible and general solution to this sort of problem.

* Wrapping up

  The goal of this example was to show that when presented a problem where
  there is no obvious algorithm that suits it, then a tool like Z3 allows you
  to describe a solution declaratively and provide a satisfying answer.

** When not to use SMT

   SAT is an NP-complete problem, and solving for richer theories does not
   reduce this complexity. So in general, SMT solving is NP-hard and not even
   decidable in all cases. If you are presented with a problem that has a known
   polynomial algorithm, then don't use an SMT solver.

   It is important to try to compartmentalize your SMT-instances; solving many
   small SMT-instances is likely to be more efficient than solving one large.
   Look for ways to divide your problem into sub-problems, and try to exclude
   the "obvious" part of a problem from the SMT-instance and move it to the
   host language (in our case Python).

   An example of where we did not follow this advice, were with the requirement
   that examiners $(i,j) \in A$ can not form a committee. Rather than encoding
   that those committees correct zero exams, we could simply remove those
   integer constants. Note that this is not a dramatic example, as the
   constraint is very simple and most likely trivial for Z3 to handle.

** When to use SMT

   If your problem is known to be NP-complete and has an elegant formulation in
   SMT, then using a tool like Z3 could be a very good idea.

   Another situation is when you currently don't know how hard the problem is.
   Specifying your problem in terms of constraints helps you understand the
   problem. Often, you will be able to solve small instances of the problem,
   which can give you insights into how you might solve the problem more
   efficiently with a more fine-tuned algorithm.

   A similar situation is when you don't exactly know what your problem is.
   This might sound like a weird situation, but my guess is that it happens
   quite frequently. Using an SMT solver as a part of a prototype gives a lot
   of flexibility because of its declarative nature. Changing your problem only
   slightly often leads to a major rewrite of your algorithm; with SMT solving,
   this is usually not the case, because it is just a matter of adding or
   removing some constraints. Once you have a well-functioning prototype, you
   can start looking for a more efficient solution if necessary.

** Exercises for the curious

   If you found this interesting, then consider solving some problems with SMT
   solving.

*** The exam committee problem

    Try to walk through the problem we have discussed here. Feel free to sneak
    a peek at the code whenever you get stuck. You might find a more efficient
    encoding or a more elegant one. Maybe you want to make it accessible
    through a web page so that this example actually ends up helping the
    administration with this problem. Play around, and let me know if you do
    something cool with it!

    Another exercise, which is by no means an easy one, is to show that this
    problem is in P or is NP-complete. Currently, we have not been able to
    prove it either way. Note that this is far from the interest area of
    IN3070, but I find it interesting, and think maybe you do too.

*** Puzzles

   Many puzzle games are NP-complete and have a nice encoding in SMT.

   Perhaps the most common example used when presenting SMT is [[https://en.wikipedia.org/wiki/Sudoku][Sudoku]]. Write
   one yourself, and if you get stuck there are many nice, and easily
   googleable, resources.

   Another example is [[https://en.wikipedia.org/wiki/Mastermind_(board_game)][Mastermind]]; if it's too hard, make the rules simpler.
   [[https://projecteuler.net/problem=185][This problem from Project Euler]] presents a simplified version of Mastermind
   and can be solved quite elegantly with Z3.

   Do you have a favorite puzzle game? See if you can model it as an SMT
   problem, and write a solver for it.

* COMMENT Local variables
  # Local Variables:
  # eval: (add-hook 'after-save-hook 'org-html-export-to-html nil t)
  # End: