Programming Z3

Satisfiability Modulo Theories and Z3

Z3

• A state-of-art Satisfiability Modulo Theories (SMT) solver
  
  
  
  
• On GitHub: https://github.com/Z3prover/z3.git
  • MIT open source license
  • Developments: Scalable simplex solver, Strings, Horn clauses, floating points, SAT extensions, word-level bit-vector reasoning, propagators

Z3 for software ++

SAT/SMT examples

• Book draft by Dennis Yurichev

• Online Programming Z3

• Online in your browser and in your other browser window.

Basic SAT

  Tie, Shirt = Bools('Tie Shirt')
  s = Solver()
  s.add(Or(Tie, Shirt), Or(Not(Tie), Shirt), Or(Not(Tie), Not(Shirt)))
  print(s.check())
  print(s.model())

Basic SMT

I = IntSort()
f = Function('f', I, I)
x, y, z = Ints('x y z')
A = Array('A',I,I)

fml = Implies(x + 2 == y, f(Select(Store(A, x, 3), y - 2)) == f(y - x + 1))

s = Solver()
s.add(Not(fml))
print(s.check())

SMT in a nutshell

Is formula satisfiable under theory ?

SMT solvers use specialized algorithms for .

Complexity

SAT/SMT

Logics and Theories

EUF

  S = DeclareSort('S')
  f = Function('f', S, S)
  x = Const('x', S)
  solve(f(f(x)) == x, f(f(f(x))) == x)
  solve(f(f(x)) == x, f(f(f(x))) == x, f(x) != x)

Linear Arithmetic

  x, y = Reals('x y')
  solve([x >= 0, Or(x + y <= 2, x + 2*y >= 6), 
                 Or(x + y >= 2, x + 2*y > 4)])

Z3 introduces auxiliary variables and represents the formula as

Only bounds (e.g., ) are asserted during search.

An Overview of Arithmetic Theories

Arrays

The declaration

  A = Array('A', IntSort(), IntSort())

  solve(A[x] == x, Store(A, x, y) == A)

which produces a solution where x necessarily equals y.

Arrays as Combinators [8]

  A = Array(Index, Elem) # array sort 
  
  a[i]             # index array 'a' at index 'i'
                   # Select(a, i)
  
  Store(a, i, v)   # update array 'a' with value 'v' at index 'i'
                   # = lambda j: If(i == j, v, a[j])
    
  Const(v, A)      # constant array
                   # = lambda j: v
  
  Map(f, a)        # map function 'f' on values of 'a'
                   # = lambda j: f(a[j])

  Ext(a, b)        # Extensionality
                   # Implies(a[Ext(a, b)] == b[Ext(a, b)], a == b)

Axiomatizing Arrays

   def array_axioms(s, index, range):
       A = ArraySort(index, range)
       a, b = Const('a b', A)
       i, j = Const('i j', index)
       v = Const('v', range)
       s.add(ForAll([a, i, j, v], Store(a, i, v)[j] == If(i == j, v, a[j])))       
       s.add(ForAll([a, b], Implies(a[Ext(a,b)] == b[Ext(a,b)], a == b)))

Bit-Vectors

def is_power_of_two(x):
    return And(x != 0, 0 == (x & (x - 1)))
x = BitVec('x', 4)
prove(is_power_of_two(x) == Or([x == 2**i for i in range(4)]))

The absolute value of a variable can be obtained using addition and xor with a sign bit.

v = BitVec('v',32)
mask = v >> 31
prove(If(v > 0, v, -v) == (v + mask) ^ mask)

Floating points

x = FP('x', FPSort(3, 4))
print(10 + x)

Sequences

s, t, u = Strings('s t u')
prove(Implies(And(PrefixOf(s, t), SuffixOf(u, t), 
                  Length(t) == Length(s) + Length(u)), 
             t == Concat(s, u)))

One can concatenate single elements to a sequence as units:

s, t = Consts('s t', SeqSort(IntSort()))
solve(Concat(s, Unit(IntVal(2))) == Concat(Unit(IntVal(1)), t))
prove(Concat(s, Unit(IntVal(2))) != Concat(Unit(IntVal(1)), s))

Algebraic Data-types

Tree = Datatype('Tree')
Tree.declare('Empty')
Tree.declare('Node', ('left', Tree), ('data', Z), ('right', Tree))
Tree = Tree.create()
t = Const('t', Tree)
solve(t != Tree.Empty)

It may produce the solution

[t = Node(Empty, 0, Empty)]

Similarly, one can prove that a tree cannot be a part of itself.

prove(t != Tree.Node(t, 0, t))

Extended SAT queries

Slide is by Joao Marques-Silva

Logical Queries - SAT+

sat unsat

model (clausal) proof

correction set core

local min correction set local min core

min correction set min core

Logical Queries

Programming Z3

Main interfaces

s = Solver()


s.add()
s.assert_and_track()

s.push()
s.pop()

s.check()
s.check()

s.model()
s.proof()
s.unsat_core()   

s.statistics()     

Assertions

s.add()
s.assert_and_track()

Scopes

s.push()

s.pop()

Check Satisfiability

s.check()

s.check()

Certificates

s.model()

s.proof()

s.unsat_core()   

Models

When s.check() returns sat Z3 can provide a model.

f = Function('f', Z, Z)
x, y = Ints('x y')
s.add(f(x) > y, f(f(y)) == y)
print(s.check())
print(s.model())

A possible model for s is:

[y = 0, x = 2, f = [0 -> 3, 3 -> 0, else -> 1]]

All SMT

def block_model(s):
    m = s.model()
    s.add(Or([ f() != m[f] for f in m.decls() if f.arity() == 0]))           

It is naturally also possible to block models based on the evaluation of only a selected set of terms, and not all constants mentioned in the model. The corresponding block_model is then

def block_model(s, terms):
    m = s.model()
    s.add(Or([t != m.eval(t, model_completion=True) for t in terms]))

A loop that cycles through multiple solutions can then be formulated:

def all_smt(s, terms):
    while sat == s.check():
       print(s.model())
       block_model(s, terms)

All SMT example

x, y = Ints('x y')
s = Solver()
s.add(1 <= x, x <= y, y <= 3)
s.push()
print("all ", x)
all_smt(s, [x])
s.pop()
print("all ", x, y)
all_smt(s, [x, y])

all  x
[x = 1, y = 1]
[x = 2, y = 2]
[x = 3, y = 3]
all  x y
[x = 3, y = 3]
[x = 2, y = 2]
[x = 1, y = 1]
[x = 1, y = 2]
[x = 1, y = 3]
[x = 2, y = 3]

All SMT - avoiding lemmas by using scopes

def all_smt(s, initial_terms):
    def block_term(s, m, t):
        s.add(t != m.eval(t, model_completion=True))
    def fix_term(s, m, t):
        s.add(t == m.eval(t, model_completion=True))
    def all_smt_rec(terms):
        if sat == s.check():
           m = s.model()
           yield m
           for i in range(len(terms)):
               s.push()
               block_term(s, m, terms[i])
               for j in range(i):
                   fix_term(s, m, terms[j])
               yield from all_smt_rec(terms[i:])
               s.pop()   
    yield from all_smt_rec(list(initial_terms))

All SMT - user propagator

class BlockTracked(UserPropagateBase):
    def __init__(self, s):
        UserPropagateBase.__init__(self, s)
        self.trail = []
        self.lim   = []
        self.add_fixed(lambda x, v : self._fixed(x, v))
        self.add_final(lambda : self._final())

    def push(self):
        self.lim += [len(self.trail)]

    def pop(self, n):
        self.trail = self.trail[0:self.lim[len(self.lim) - n]]
        self.lim = self.lim[0:len(self.lim)-n]

    def _fixed(self, x, v):    
        self.trail += [(x, v)]

    def _final(self):
        print(self.trail)
        self.conflict([x for x, v in self.trail])

All SMT - user propagator instance

s = SimpleSolver()
b = BlockTracked(s)

x, y, z, u = Bools('x y z u')
b.add(x)
b.add(y)
b.add(z)

s.add(Or(x, Not(y)), Or(z, u), Or(Not(z), x))

print(s.check())

python propagator.py 
[(z, True), (x, True), (y, False)]
[(z, True), (x, True), (y, True)]
[(z, False), (y, False), (x, False)]
[(z, False), (y, False), (x, True)]
[(z, False), (y, True), (x, True)]
unsat

Other Propagation Facilities

 add_eq
 propagate

Cores, Correction Sets, Satisfying Assignments

• (M)US (Minimal) unsatisfiable subset

  • (minimal) subset of assertions that are unsatisfiable. Also known as a core

• (M)SS (Maximal) satisfiable subset

  • (maximal) subset of assertions that are satisfiable.

• (M)CS (Minimal) correction set

  • complement of an (M)SS.

• (Prime) implicant

  • iff is a core for .

Maximizing Satisfying Assignments [7] - background

• Find maximal subset such that is satisfiable.

• Idea: speed up checks by assuming soft constraints that cannot be added to .

• Initialize , .

• For each

  • If is SAT, set
  • Otherwise, add

Maximizing Satisfying Assignments [7] - code

def tt(s, f): 
    return is_true(s.model().eval(f))

def get_mss(s, ps):
    if sat != s.check():
       return []
    mss = { q for q in ps if tt(s, q) }
    return get_mss(s, mss, ps)

def get_mss(s, mss, ps):
    ps = ps - mss
    backbones = set([])
    while len(ps) > 0:
       p = ps.pop()
       if sat == s.check(mss | backbones | { p }):
          mss = mss | { p } | { q for q in ps if tt(s, q) }
          ps  = ps - mss
       else:
          backbones = backbones | { Not(p) }
    return mss          

Unsatisfiable Cores

def get_mus(s, seed):
    return s.unsat_core()

Minimizing Cores [1, 2, 4, 6]

Use built-in core minimization:

s.set("sat.core.minimize","true")  # For Bit-vector theories
s.set("smt.core.minimize","true")  # For general SMT 

Or roll your own:

def quick_explain(test, sub):
    return qx(test, set([]), set([]), sub)

def qx(test, B, D, C):
    if {} != D:
       if test(B):
          return set([])
    if len(C) == 1:
       return C
    C1, C2 = split(C)
    D1 = qx(test, B | C1, C1, C2)
    D2 = qx(test, B | D1, D1, C1)
    return D1 | D2
    
def test(s):
    return lambda S: s.check([f for f in S]) == unsat


s = Solver()
a, b, c, d, e = Bools('a b c d e')
s.add(Or(a, b))
s.add(Or(Not(a), Not(b)))
s.add(Or(b, c))
s.add(Or(Not(c), Not(a)))

print s.check([a,b,c,d,e])
print s.unsat_core()

mus = quick_explain(test(s), {a,b,c,d})

A Duality: MUS MCS [10]

All Cores and Correction Sets [5]

def ff(s, p): 
    return is_false(s.model().eval(p))

def marco(s, ps):
    map = Solver()
    while map.check() == sat:
        seed = {p for p in ps if not ff(map, p)}
        if s.check(seed) == sat:
           mss = get_mss(s, seed, ps)
           map.add(Or(ps - mss))
           yield "MSS", mss
        else:
           mus = get_mus(s, seed)
           map.add(Not(And(mus)))
           yield "MUS", mus

Bounded Model Checking

Bounded Model Checking example

Let us check whether there is some , such that when numbers are represented using 4 bits. The corresponding transition system uses a state variable x0 which is named x1 in the next state. Initially x0 == 0 and in each step the variable is incremented by 3. The goal state is x0 == 10.

x0, x1 = Consts('x0 x1', BitVecSort(4))
bmc(x0 == 0, x1 == x0 + 3, x0 == 10, [], [x0], [x1])

Pebble games

x, y, z, u, a, b, c, d, e, f = Bools('x y z u a b c d e f')
x0, y0, z0, u0, a0, b0, c0, d0, e0, f0 = Bools('x0 y0 z0 u0 a0 b0 c0 d0 e0 f0')
x1, y1, z1, u1, a1, b1, c1, d1, e1, f1 = Bools('x1 y1 z1 u1 a1 b1 c1 d1 e1 f1')
init = And([x0, y0, z0, u0, Not(a0), Not(b0), Not(c0), Not(d0), Not(e0), Not(f0)])
final = And([x0, y0, z0, u0, Not(a0), Not(b0), Not(c0), Not(d0), e0, f0])
edges = [(x0,f0),(y0,a0),(z0,a0),(z0,c0),(z0,d0),(z0,b0),(u0,b0),(a0,f0),
         (a0,c0),(b0,d0),(d0,e0),(c0,e0)]
next = { x0 : x1, y0 : y1, z0 : z1, u0 : u1, a0 : a1, b0 : b1, c0 : c1,
         d0 : d1, e0 : e1, f0 : f1 }

Pebble Game Transitions

x, y, z, u, a, b, c, d, e, f = Bools('x y z u a b c d e f')
x0, y0, z0, u0, a0, b0, c0, d0, e0, f0 = Bools('x0 y0 z0 u0 a0 b0 c0 d0 e0 f0')
x1, y1, z1, u1, a1, b1, c1, d1, e1, f1 = Bools('x1 y1 z1 u1 a1 b1 c1 d1 e1 f1')
init = And([x0, y0, z0, u0, Not(a0), Not(b0), Not(c0), Not(d0), Not(e0), Not(f0)])
final = And([x0, y0, z0, u0, Not(a0), Not(b0), Not(c0), Not(d0), e0, f0])
edges = [(x0,f0),(y0,a0),(z0,a0),(z0,c0),(z0,d0),(z0,b0),(u0,b0),(a0,f0),
         (a0,c0),(b0,d0),(d0,e0),(c0,e0)]
next = { x0 : x1, y0 : y1, z0 : z1, u0 : u1, a0 : a1, b0 : b1, c0 : c1,
         d0 : d1, e0 : e1, f0 : f1 }

If the pebbling state changes on a node, then both children must be pebbled.

add_pebble = [Implies(Xor(parent, next[parent]), And(child, next[child]))
                 for (child, parent) in edges]

def transition(bound):
    max_pebbles = AtMost(x0,y0,z0,u0,a0,b0,c0,d0,e0,f0, bound)
    return And(add_pebble + [max_pebbles])

Encoding Bounded Model Checking

index = 0
def fresh(s):
    global index
    index += 1
    return Const("!f%d" % index, s)

def zipp(xs, ys):
    return [p for p in zip(xs, ys)]

def bmc(init, trans, goal, fvs, xs, xns):
    s = Solver()
    s.add(init)
    count = 0
    while True:
        print("iteration ", count)
        count += 1
        p = fresh(BoolSort())
        s.add(Implies(p, goal))
        if sat == s.check(p):
            print (s.model())
            return
        s.add(trans)
        ys = [fresh(x.sort()) for x in xs]
        nfvs = [fresh(x.sort()) for x in fvs]
        trans = substitute(trans, 
                           zipp(xns + xs + fvs, ys + xns + nfvs))
        goal = substitute(goal, zipp(xs, xns))
        xs, xns, fvs = xns, ys, nfvs

Propositional Interpolation

I am given solvers A and B

• A.add(B.assertions()) is unsat

• A, B share propositional atoms xs

• Compute reverse interpolant I

  • Implies(B, I)

  • Implies(A, Not(I))

Propositional Interpolation Example

A = SolverFor("QF_FD")
B = SolverFor("QF_FD")
a1, a2, b1, b2, x1, x2 = Bools('a1 a2 b1 b2 x1 x2')
A.add(a1 == x1, a2 != a1, a2 != x2)
B.add(b1 == x1, b2 != b1, b2 == x2)
print(list(interpolate(A, B, [x1, x2])))

    [Not(And(Not(x2), Not(x1))), Not(And(x2, x1))]

Propositional Interpolation - code

def mk_lit(m, x):
    if is_true(m.eval(x)):
       return x
    else:
       return Not(x)

def pogo(A, B, xs):   
    while sat == A.check():
       m = A.model()
       L = [mk_lit(m, x) for x in xs]
       if unsat == B.check(L):
          notL = Not(And(B.unsat_core()))
          yield notL
          A.add(notL)
       else:
          print("expecting unsat")
          break

interpolate = pogo

MaxSAT

MaxSAT example

a, b, c = Bools('a b c')
opt = Optimize()
opt.add_soft(a, 1)
opt.add_soft(b, 2)
opt.add_soft(c, 3)
opt.add(a == c)
opt.add(Not(And(a, b)))
print(opt.check())
print(opt.model())

MaxSAT flattened

a, b, c = Bools('a b c')
opt = Optimize()
opt.add_soft(a, 1)
opt.add_soft(b, 1); opt.add_soft(b, 1)
opt.add_soft(c, 1); opt.add_soft(c, 1); opt.add_soft(c, 1)
opt.add(a == c)
opt.add(Not(And(a, b)))
print(opt.check())
print(opt.model())

• NB. Implementations of MaxSAT flatten weights on demand.

Core Guided MaxSAT

The initial idea is to replace a core

by soft constraints such that cost of new system is decreased.

MaxSAT with Cores - example

MaxSAT with Cores - example (2)

MaxSAT with Cores [9]

A:

A':

Lemma: for every model of ,

Proof: min:

Core Guided MaxSAT - maxres

def add_def(s, fml):
    name = Bool("%s" % fml)
    s.add(name == fml)
    return name

def relax_core(s, core, Fs):
    prefix = BoolVal(True)
    Fs -= { f for f in core }
    for i in range(len(core)-1):
        prefix = add_def(s, And(core[i], prefix))
        Fs |= { add_def(s, Or(prefix, core[i+1])) }

def maxsat(s, Fs):
    cost = 0
    Fs0 = Fs.copy()
    while unsat == s.check(Fs):
        cost += 1
        relax_core(s, s.unsat_core(), Fs)    
    return cost, { f for f in Fs0 if tt(s, f) }

Other MaxSAT algorithms

• PM - original core guided MaxSAT algorithm by Fu and Malik
• RC2 - state-of-art replaces core by a single relaxation bound (fewer soft constraints)
• MaxHS - does not relax soft constraints, instead enumerates cores and finds minimal hitting sets

A different twist on MaxRes

def relax_core(s, core, Fs):
    Fs -= { f for f in core }
    while len(core) >= 2:
        v = add_def(s, And(core[0], core[1]))
        u = add_def(s, Or(core[0], core[1]))
        Fs |= { u }
        core = core[2:] + [ v ]

Loop invariant:

• The cost of core | Fs is constant
• core has cost at least 1

On exit, len(core)-1 elements are re-added to Fs

Hypothesis: allow different clusters of soft constraints.

Quantifiers

Quantifier Solving - algorithm

• Maintain solver state

  • Corresponding to instances of with instantiated.

• Let be a value satisfying .

• If is unsat, is sat.

• Otherwise, let

  • .

• If becomes unsat, then is unsat.

Quantifier Solving - code

Modified from an example in pysmt.

from z3.z3util import get_vars


def efsmt(ys, phi, maxloops = None):
    """Solving exists xs. forall ys. phi(x, y)"""
    xs = [x for x in get_vars(phi) if x not in ys]
    E = Solver()
    F = Solver()
    E.add(BoolVal(True))
    loops = 0
    while maxloops is None or loops <= maxloops:
        loops += 1
        eres = E.check()
        if eres == sat:
            emodel = E.model()
            sub_phi = substitute(phi, [(x, emodel.eval(x, True)) for x in xs])
            F.push()
            F.add(Not(sub_phi))
            fres = F.check()
            if fres == sat:
                fmodel = F.model()
                sub_phi = substitute(phi, [(y, fmodel.eval(y, True)) for y in ys])
                E.add(sub_phi)
            else:
                return fres, [(x, emodel.eval(x, True)) for x in xs]
            F.pop()
        else:
            return eres
    return unknown

Quantifier Solving - example

x, y, z = Reals("x y z")
print(efsmt([y], Implies(And(y > 0, y < 10), y - 2 * x < 7)))
print(efsmt([y], And(y > 3, x == 1)))

Output:

(sat, [(x, 3/2)])
unsat

Monadic Decomposition

• How can I convert a formula into monadic form ?

Monadic Decomposition - idea

• Keep sampling values from solver state .

• Update solver state to exclude COI of sample:

• Decompose by two cases:

• Sample within the cases.

Monadic Decomposition - code 1

from z3 import *
def nu_ab(R, x, y, a, b):
    x_ = [ Const("x_%d" %i, x[i].sort()) for i in range(len(x))]
    y_ = [ Const("y_%d" %i, y[i].sort()) for i in range(len(y))]
    return Or(Exists(y_, R(x+y_) != R(a+y_)), Exists(x_, R(x_+y) != R(x_+b)))

def isUnsat(fml):
    s = Solver(); s.add(fml); return unsat == s.check()

def lastSat(s, m, fmls):
    if len(fmls) == 0: return m
    s.push(); s.add(fmls[0])
    if s.check() == sat: m = lastSat(s, s.model(), fmls[1:])       
    s.pop(); return m

Monadic Decomposition - code 2

def mondec(R, variables):
    print(variables)
    phi = R(variables);    
    if len(variables)==1: return phi
    l = int(len(variables)/2)
    x, y = variables[0:l], variables[l:]
    def dec(nu, pi):
        if isUnsat(And(pi, phi)): 
           return BoolVal(False)
        if isUnsat(And(pi, Not(phi))): 
           return BoolVal(True)
        fmls = [BoolVal(True), phi, pi] 
        #try to extend nu
        m = lastSat(nu, None, fmls)              
        #nu must be consistent 
        assert(m != None)                         
        a = [ m.evaluate(z, True) for z in x ]
        b = [ m.evaluate(z, True) for z in y ]
        psi_ab = And(R(a+y), R(x+b))
        phi_a = mondec(lambda z: R(a+z), y)
        phi_b = mondec(lambda z: R(z+b), x)
        nu.push()
        #exclude: x~a and y~b
        nu.add(nu_ab(R, x, y, a, b))              
        t = dec(nu, And(pi, psi_ab)) 
        f = dec(nu, And(pi, Not(psi_ab)))
        nu.pop()
        return If(And(phi_a, phi_b), t, f)
    #nu is initially true
    return dec(Solver(), BoolVal(True))             

Cube and Conquer

Using Cubes

s = SolverFor("QF_FD")
s.add()
s.set("sat.restart.max", 100)
def cube_and_conquer(s):
    for cube in s.cube():
       if len(cube) == 0:
          return unknown
       if is_true(cube[0]):
          return sat     
       is_sat = s.check(cube):
       if is_sat == unknown:
          s1 = s.translate(s.ctx)
          s1.add(cube)
          is_sat = cube_and_conquer(s1)
       if is_sat != unsat:
          return is_sat
    return unsat

Application: Product Configuration

Query

• Find all unit literals - with explanations
• Useful for finding fixed parameters [3]

a, b, c, d = Bools('a b c d')

s = Solver()
s.add(Implies(a, b), Implies(c, d))   # background formula
print s.consequences([a, c],          # assumptions
                     [b, c, d])       # what is implied?

(sat, [Implies(c, c), Implies(a, b), Implies(c, d)])

Algorithm

1. Assert hard constraints
2. Add assumptions at
3. Let model of
4. Case splits:
5. Unit propagation
   • if conflict, remove fixed variable with explanation from
6. Run CDCL until next restart or sat
   • if sat: use resulting model to remove non-fixed from
   • if , exit, else goto 4

Modeling with SMT

• Most examples we have surveyed are essentially propositional SAT.

• But SMT > SAT

• Advantage of SMT is at the modeling level.

Modeling with EUF

• Assign items to urns such that each urn has at most two items.

• With indicator variables:

  • if item is assigned to urn .
  • Size of encoding: .

• With EUF

Some other fragments

• EPR

  • Decidable and succinct
  • Can encode many domains

• Algebraic Datatypes and EUF

  • ADTs have built-in occurs check to enforce acyclicity
  • Example: Read-only fields in OO classes.