A Sketch of Categorical Relation Algebra Combinators in Z3Py

I’ve discussed relation algebra before. Relations are sets of tuples. There, I implemented the relations naively using lists for sets. This is very simple, and very clean especially with list comprehension syntax. It is however horrifically inefficient, and we could only deal with finitely enumerable domains. The easiest path to fixing these problems is to cash out to an external solver, in this case z3.

There are many beautifully implemented solvers out there and equally beautiful DSL/modeling languages. Examples in mind include sympy, cvxpy, and z3. These modeling languages require you to instantiate variable objects and build expressions out of them and then hand it off to the solver. This is a reasonable interface, but there are advantages to a more categorical/point-free style DSL.

Point-free languages are ones that do not include binding forms that introduce bound/dummy variables. Examples of binding forms like this are \lambda \sum \max \min \int \sup \lim \forall \exists. One problem lies in the fact that the names of bound variables don’t matter, and that they end up accidentally smashing into each other. You may have experienced this in physics or math class as the dummy indices or dummy variable problem causing you to screw up your calculation of some cross product identity or some complicated tensor sum. These are surprisingly subtle problems, very difficult to diagnose and get right. de Bruijn indices is a technique for giving the bound variables canonical names, but it sucks to implement in its own way. When you make a DSL point free, it is a joy to manipulate, optimize, and search. I think this may be the core of why category theory is good language for mathematics and programming.

Point-free style also tends to have significant economy of size, for better or worse. Sometimes it is better to have an expression very dense in information. This is important if you are about the algebraically manipulate an expression with paper and pencil. Every manipulation requires a great deal of mind numbing copying as you proceed line by line, so it can be excruciating if your notation has a lot of unnecessary bulk.

Relations are like functions. The two pieces of the tuple can be roughly thought of as the “input” and the “output”. Relations are only loosely directional though. Part of the point of relations is that the converse (inverse) of a relation is easy to define.

When we are implement relations, we have a choice. Do we want the relation to produce its variables, accept its variable, or accept one and produce the other? There are advantages to each. When relations were [(a,b)], a -> b -> Bool, and a -> [b] converting between these forms was a rather painful enumeration process. The sting of converting between them is taken out by the fact that the conversion is no longer a very computationally expensive process, since we’re working at the modeling layer.

When you’re converting a pointful DSL to pointfree DSL, you have to be careful where you instantiate fresh variables or else you’ll end up with secret relations that you didn’t intend. Every instantiation of id needs to be using fresh variables for example. You don’t want the different id talking to each other. Sometimes achieving this involves a little currying and/or thunking.

There is a pattern that I have notice when I’m using modeling languages. When you have a function or relation on variables, there are constraints produced that you have to keep a record of. The pythonic way is to have a Model or Solver object, and then have that objects mutate an internal record of the set of constraints. I don’t particularly enjoy this style though. It feels like too much boiler plate.

In Haskell, I would use something like a Writer monad to automatically record the constraints that are occurring. Monads are not really all that pleasant even in Haskell, and especially a no go in python without “do” syntax.

However, because we are going point free it is no extra cost at all to include this pipework along for the ride in the composition operation.

Here are implementations of the identity and composition for three different styles. Style 1 is fully receptive, style 2 is mixed (function feeling) and style 3 is fully productive of variables.

Fair warning, I’m being sketchy here. I haven’t really tried this stuff out.

z3 is a simply typed language. You can get away with some polymorphism at the python level (for example the == dispatches correctly accord to the object) but sometimes you need to manually specify the sort of the variables. Given these types, the different styles are interconvertible

We can create the general cadre of relation algebra operators. Here is a somewhat incomplete list

Questions about relation algebra expressions are often phrased in term of relational inclusion. You can construct a relation algebra expression, use the rsub in the appropriate way and ask z3’s prove function if it is true.

Z3 has solvers for

  • Combinatorial Relations
  • Linear Relations
  • Polyhedral Relations
  • Polynomial Relations
  • Interval Relations – A point I was confused on. I thought interval relations were not interesting. But I was interpetting the term incorrectly. I was thinking of relations on AxB that are constrained to take the form of a product of intervals. In this case, the choice of A has no effect on the possible B whatsoever, so it feels non relational. However, there is also I_A x I_B , relations over the intervals of A and B. This is much closer to what is actually being discussed in interval arithmetic.

Applications we can use this for:

  • Graph Problems. The Edges can be thought of as a relation between vertices. Relation composition Using the starn operator is a way to ask for paths through the graph.
  • Linear Relations – To some degree this might supplant my efforts on linear relations. http://www.philipzucker.com/linear-relation-algebra-of-circuits-with-hmatrix/ Z3 is fully capable of understanding linear relations.
  • Safety and liveness of control systems. Again. a transition relation is natural here. It is conceivable that the state space can be heterogenous in time, which is the interesting power of the categorical style. I feel like traditional control systems usually maintain the same state space throughout.
  • Program verification
  • Games? Nash equilibria?

Other Thoughts

  • Maybe we are just building a shitty version of alloy. https://alloytools.org/
  • What about uninterpeted relations? What about higher order relations? What about reflecting into a z3 ADT for a relational language. Then we could do relational program synthesis. This is one style, just hand everything off to smt. https://github.com/nadia-polikarpova/cse291-program-synthesis/tree/master/lectures
  • I should try to comply with python conventions, in particular numpy and pandas conventions. @ should be composition for example, since relation composition has a lot of flavor of matrix composition. I should overload a lot of operators, but then I’d need to wrap in a class ๐Ÿ™
  • Z3 has special support for some relations. How does that play in? https://theory.stanford.edu/~nikolaj/programmingz3.html#sec-special-relations https://z3prover.github.io/api/html/ml/Z3.Relation.html
  • As long as you only use composition, there is a chaining of existentials, which really isn’t so bad.
  • What we’ve done here is basically analogous/identical to what John Wiegley did compiling to the category of z3. Slightly different in that he only allowed for existential composition rather than relational division. http://newartisans.com/2017/04/haskell-and-z3/
  • We can reduced the burden on z3 if we know the constructive proof objects that witness our various operations. Z3 is gonna do better if we can tell it exactly which y witness the composition of operators, or clues to which branch of an Or it should use.
  • It’s a bummer, but when you use quantifiers, you don’t see countermodels? Maybe there is some hook where you can, or in the dump of the proof object.
  • What about recursion schemes? The exact nature of z3 to handle unbounded problems is fuzzy to me. It does have the support to define recursive functions. Also explicit induction predicates can go through sometimes. Maybe look at the Cata I made in fancy relaion algebra post
  • I think most proof assistants have implementations of relation algebra available. I find you can do a surprising amount in z3.

More Stupid Z3Py Tricks: Simple Proofs

Z3 can be used for proofs. The input language isn’t anywhere near as powerful as interactive theorem provers like Coq, Isabelle, or Agda, but you can ask Z3 to prove pretty interesting things. Although the theorems that follow aren’t hard in interactive theorem provers, they would take beyond complete novice level skills to state or prove.

I like to think of the z3 proving process as “failing to find a counterexample”. Z3py has supplies a function prove which is implemented like this.

Basically, it negates the thing you want to prove. It then tries to find a way to instantiate the variables in the expression to make the statement false. If it comes back unsat, then there is no variable assignment that does it. Another way to think about this is rewriting the \forall y. p(y) as \neg \exists y \neg p (y). The first \neg lives at sort of a meta level, where we consider unsat as a success, but the inner \neg is the one appearing in s.add(Not(claim)).

We can prove some simple facts. This is still quite cool, let’s not get too jaded. Manually proving these things in Coq does suck (although is easy if you use the ring, psatz, and lra tactics https://coq.inria.fr/refman/addendum/micromega.html, which you DEFINITELY should. It is a great irony of learning coq that you cut your teeth on theorems that you shouldn’t do by hand).

Ok, here’s our first sort of interesting example. Some properties of even and odd numbers. Even and Odd are natural predicates. What are possible choices to represent predictaes in z3?
We can either choose python functions IntSort -> BoolSort() as predicates or we can make internal z3 functions Function(IntSort(), BoolSort())

All well and good, but try to prove facts about the multiplicative properties of even and odd. Doesn’t go through. ๐Ÿ™

Here’s a simple inductive proof. Z3 can do induction, but you sort of have to do it manually, or with a combinator. Given a predicate f, inductionNat returns

Here’s another cute and stupid trick. Z3 doesn’t have a built in sine or cosine. Perhaps you would want to look into dreal if you think you might be heavily looking into such things. However, sine and cosine are actually defined implicitly via a couple of their formula. So we can instantiate
A slightly counterintuitive thing is that we can’t use this to directly compute sine and cosine values. That would require returning a model, which would include a model of sine and cosine, which z3 cannot express.
However, we can try to assert false facts about sine and cosine and z3 can prove they are in fact unsatisfiable. In this way we can narrow down values by bisection guessing. This is very silly.

A trick that I like to use sometimes is embedding objects in numpy arrays. Numpy slicing is the best thing since sliced bread. A lot, but not all, of numpy operations come for free, like matrix multiply, dot, sum, indexing, slicing, reshaping. Only some are implemented in terms of overloadable operations. here we can prove the Cauchy Schwartz inequality for a particular vector and some axioms of vector spaces.

Defining and proving simple properties of Min and Max functions

Proving the Babylonian method for calculating square roots is getting close to the right answer. I like the to think of the Babylonian method very roughly this way: If your current guess is low for the square root x/guess is high. If your guess is high, x/guess is low. So if you take the average of the two, it seems plausible you’re closer to the real answer. We can also see that if you are precisely at the square root, (x/res + x)/2 stays the same. Part of the the trick here is that z3 can understand square roots directly as a specification. Also note because of python overloading, babylonian with work on regular numbers and symbolic z3 numbers. We can also prove that babylon_iter is a contractive, which is interesting in it’s own right.

A funny thing we can do is define interval arithmetic using z3 variables. Interval arithmetic is very cool. Checkout Moore’s book, it’s good. This might be a nice way of proving facts related to real analysis. Not sure.
This is funny because z3 internally uses interval arithmetic. So what we’re doing is either very idiotically circular or pleasantly self-similar.
We could use a similar arrangement to get complex numbers, which z3 does not natively support

Stupid Z3Py Tricks: Verifying Sorting Networks off of Wikipedia

Sorting networks are a circuit flavored take on sorting. Although you can build circuits for any size input, any particular circuit works for a fixed sized input. They are like an unrolling of the loops or recursion of more familiar sorting algorithms. They come up also in the context of parallel and gpu sorting

Here’s an interesting thing. We can go to Wikipedia and get a little python snippet for the comparison order of a Batcher odd-even mergesort. Kind of a confusing algorithm. Why does it even work? Is it even right? It’s written in some kind of funky, indexful generator style.

Source: Wikipedia

Well we can confirm this relatively straightforwardly using z3 by replacing the implementation of compare_and_swap with its z3 equivalent. We then ask z3 .

This comes back unsat, hence there are no inputs or executions that do not come back sorted. If I delete some elements from pair_to_compare, it comes back sat, showing that it doesn’t always sort.

The trick here is that the circuit is fixed size, so we have no need for induction, one of the main things z3 is rather finicky at.

It’s somewhat interesting to note that the output of odd_even_merge is a sequence of instructions, we can think of this as interpreting a very small 1 instruction programming language.

We can also confirm similarly a simple odd-even bubblesort and other similar algorithms.

Q: What about using uninterpreted sorts rather than integers? Integers is pretty convincing to me.

same_elems is slightly weaker than a permutation predicate. Wasn’t super obvious to me the best way to do a permutation predicate in z3. Would I want to internalize the array?

Edit: Upon further thought, actually the sort IS a nice predicate for permutation. How do we compute if two things are permutations of each other? By sorting them and forcing a zipped equality. Alternatively count the number of each element (a piece of bucket sort). Since this sort is done by composing swaps, it is somewhat intrinsically a permutation

As a bummer though, I think randomized testing on arrays would be equally or perhaps more convincing of the correctness of the algorithm. Oh well.

Programming and Interactive Proving With Z3Py

I’ve been fiddling with z3py, figuring out some functionality and realizing some interesting things you could do with it. I think I’m at a point where it is nice to checkpoint myself with a blog post.

I’m a little surprised z3py doesn’t overload the & and | operators and some kind of implies operator for BoolRef. You can insert them later using this.

from z3 import *
# useful non default operator definitions for z3 bools
BoolRef.__and__ = lambda self, rhs: And(self,rhs)
BoolRef.__or__ = lambda self, rhs: Or(self,rhs)
BoolRef.__xor__ = lambda self, rhs: Xor(self,rhs)
BoolRef.__invert__ = lambda self: Not(self)
BoolRef.__rshift__ = lambda self, rhs: Implies(self,rhs)

Functional Programming

Python is not the best functional programming environment imo. And by functional programming I implicitly mean roughly ML-like FP a la Haskell or OCaml. I don’t venture much into lisp land.

The lack of good algebraic datatypes (the class syntax is so ungainly) and a type system hurts. The lack of pattern matching hurts. The lambda keyword is so long it makes me sad.

But you have full access to z3 from the python bindings. Z3 does have algebraic data types, and a type system. It has built in substitution mechanisms and evaluation. And it has insane search procedures and the ability to prove things. Pretty damn cool!

Unfortunately the type system is rather simplistic, being basically simply typed rather than polymorphic or something else. But using python a a schema/macro system for z3 seems a plausible way forward.

To build templated types, you can have constructor functions in python for the appropriate types.

def Tuple(a,b):
    Type = Datatype('Tuple(f{(a.name(),b.name())})')
    Type.declare('pair', ('fst', a), ('snd', b))
    Type = Type.create()
    return Type
def Either(a,b):
    Type = Datatype('Either(f{(a.name(),b.name())})')
    Type.declare('left', ('getLeft', a))
    Type.declare('right', ('getRight', b))
    Type = Type.create()
    return Type 
def Maybe(a):
    Type = Datatype('Maybe(f{(a.name())})')
    Type.declare('Just', ('fromJust', a))
    Type.declare("Nothing")
    Type = Type.create()
    return Type
def List(a):
    Type = Datatype('List(f{(a.name())})')
    Type.declare('Cons', ('car', a), ('cdr', Type))
    Type.declare("Nil")
    Type = Type.create()
    return Type
''' 
Note in regards to List. Z3 has a built in type Seq that I think it has built in smarts about. You might be better off using that rather than a custom List. YMMV
'''

You can access the constructors from the returned types. Check this out. You get detector functions is_Nothing and is_Just , the extractor function fromJust and constructor functions Nothing and Just. I do a lot of dir exploration with z3py. It’s hard to know what’s available sometimes

# dir(Maybe(IntSort())) returns
[
'Just',
 'Nothing',
... underscore junk ... ,
 'accessor',
 'as_ast',
 'ast',
 'cast',
 'constructor',
 'ctx',
 'ctx_ref',
 'eq',
 'fromJust',
 'get_id',
 'hash',
 'is_Just',
 'is_Nothing',
 'kind',
 'name',
 'num_constructors',
 'recognizer',
 'sexpr',
 'subsort',
 'translate',
 'use_pp']

It’s possible to build a general purpose match combinator on these types since you can introspect the number of constructors of the ADT using num_constructors, constructor, recognizer, and accessor. There might be a match inside z3py somewhere? I think it’s part of the SMTLIB standard now.

def match(x, **kwargs):
    t = x.sort()
    nc = t.num_constructors()
    acc = kwargs["_"] # default argument
    for c in range(nc):
        con = t.constructor(c)
        rec = t.recognizer(c)
        nfields = con.arity()
        if nfields == 0:
            res = kwargs[con.name()]
        else:
            res = kwargs[con.name()](  *[t.accessor(c,a)(x) for a in range(nfields)] )
        acc = If(rec(x), res, acc)
    return acc

Example usage:

match(Const("x", Maybe(IntSort())), Just=lambda y : y + 1, Nothing = IntVal(3), _=IntVal(10))
# returns If(is(Nothing, x), 3, If(is(Just, x), fromJust(x) + 1, 10))

Z3 has a substitution mechanism built in. This is useful for instantiating ForAll and for evaluating Lambda. The substitute_vars function is what you want like so substitute_vars(f.body(), x, y, z)

It is possible to reflect the syntax in a fairly straightforward way back into python via a lambdify function, mimicking the equivalent very useful function from sympy. Lambdify is basically an interp function. Here is a start for such a function. I by no means have implemented interpretation of the entirety of z3. Also I feel like this implementation is very clunky. Some kind of CPS?

def lift1(f,x):
    return lambda *args: f(x(*args))

def lift2(op,l,r):
    return lambda *args: op(l(*args), r(*args))

# interp is useful for transferring expressions into numpy, sympy
# but also for program extraction

from functools import reduce
import operator as op
def interp(a, *args):
    if is_true(a):
        return lambda *args: True
    elif is_false(a):
        return lambda *args: False
    elif is_int_value(a):
        return lambda *args: a.as_long()
    elif is_rational_value(a):
        n = a.numerator_as_long()
        d = a.denominator_as_long()
        return lambda *args: n / d
    #elif is_algebraic_value(a):
    #    pass
    elif is_const(a): # is free variable
        loc = [ind for ind, b in enumerate(args) if eq(a,b)]
        if len(loc) == 0:
            return a
        else:
            ind = loc[0]
            return lambda *args2: args2[ind]   
    b = [interp(c, *args) for c in a.children()]
    if is_and(a):
        return lambda *args: reduce(op.and_, [f(*args) for f in b])
    elif is_or(a):
        return lambda *args: reduce(op.or_, [f(*args) for f in b])
    elif is_app_of(a, Z3_OP_XOR):
        return lambda *args: reduce(op.xor, [f(*args) for f in b])
    elif is_add(a):
        return lambda *args: reduce(op.add, [f(*args) for f in b])
    elif is_mul(a):
        return lambda *args: reduce(op.mul, [f(*args) for f in b])
    elif len(b) == 1:
        n = b[0]    
        if is_not(a):
            return lift1(op.invert , n)# lambda *args: ~n
    elif len(b) == 2:
        l,r = b
        if is_eq(a):
            return lift2(op.eq, l,r) #lambda *args: l == r
        elif is_distinct(a): # This can be multi_argument
            return lift2(op.ne, l,r) # lambda *args: l != r
        elif is_sub(a):
            return lift2(op.sub, l,r) # lambda *args: l - r
        elif is_app_of(a, Z3_OP_POWER):
            return lift2(op.pow, l,r) #lambda *args: l ** r
        elif is_div(a):
            return  lift2(op.truediv, l,r)# lambda *args: l / r
        elif is_idiv(a):
            return lift2(op.floordiv, l,r) # lambda *args: l // r
        elif is_mod(a):
            return lift2(op.mod, l,r) # lambda *args: l % r
        elif is_le(a):
            return lift2(op.le, l,r) # lambda *args: l <= r
        elif is_lt(a):
            return lift2(op.lt, l,r) # lambda *args: l < r
        elif is_ge(a):
            return lift2(op.ge, l,r) #lambda *args: l > r
        elif is_gt(a):
            return lift2(op.gt, l,r) # lambda *args: l >= r
        elif is_implies(a):
            return lambda *args: (~ l(*args) ) & r(*args) 
    print("unrecognized constructor: ", type(a))
    assert(False)
#example usage
a = Bool('a')
interp(Xor(a & a | a, a), a)(True)
x, y = Reals('x y')
interp(x + y + y + y * x, x ,y)(3,2)

There is the ability to define recursive functions in z3. It is also plausible to define them via. In this way you can get a reversible functional programming language, maybe some subset of mercury / curry’s power.

fac = RecFunction('fac', IntSort(), IntSort())
n = Int('n')
RecAddDefinition(fac, n, If(n == 0, 1, n*fac(n-1)))

s = Solver()
s.add(fac(n) == 6)
s.check()
s.model()
#  returns [n = 3, fac = [0 โ†’ 1, else โ†’ fac(-1 + ฮฝ0)ยทฮฝ0]]

Interactive Theorem Proving

Z3 is awesome at thoerem proving. But somethings it just doesn’t handle right and needs human guidance.

Through searching, there are a couple interesting python interactive theorem prover projects. Cody pointed me to a project he worked on a while back, Boole https://github.com/avigad/boole . It has a dependently typed lambda calculus in it with the purpose of gluing together many systems, I think. He implemented a lot of stuff from scratch. I think I want to try to get less and do less. There is also holpy https://arxiv.org/abs/1905.05970 which appears to be being actively developed. It’s roughly a translation of hol to python I think. It’s available from a strange chinese github on the author’s website if you go looking for it.

This suggests an interesting approach. Most interactive theorem provers start unautomated and add it later. Instead we can iteratively build an interface to de-automate z3.

Altogether, this approach is more HOL flavored than Coq/Agda flavored. z3 terms are our logic and python is our manipulation metal language. Ideally, one would want to verify that every.

Python is so unprincipled that I can’t imagine that you could ever build a system up to the trustworthiness of the other theorem provers. But this is freeing in many ways. Since that is off the table, we can just do the best we can.

Using the z3 syntax tree and the z3 proof automation and z3 substitution mechanisms gives us a HUGE step up from implementing them from scratch. Ideally, we’d want to write as little python as possible, and especially as little python as possible that has to be trusted to be implemented correctly.

One big concern is accidental mutation of the proof under our feet by python. Perhaps using hashes and checking them might be a way to at least detect this. I need to have a good think about how to factor out a trusted core from all possible tactics.

I think it helps a little that z3 often will be able to verify the equivalence of small steps in proofs even if it can’t do the entire proof itself.

I think induction principles will need to be injected by hand. Z3 doesn’t really have that built in. There are definitely situations that after you introduce the induction, z3 can slam all the cases no problem. For example, check this out.

Another thing that might be nice is integration/translation to sympy. Sympy has a ton of useful functionality, at the very least differentiation.

Translation and integration with cvxpy for sum of squares proofs would also be quite neat. I already did something with this using sympy. I’m not super sure how you extract exact proofs from the floating point solutions SCS returns? I think there is a thing. I’ve heard the LLL algorithm can be used for this somehow (finding likely algebraic number matches to floating point numbers)?

So here are some sketched out ideas for tactics.

class Proof():
    def __init__(self, goal, name=None): # Taken a name for the theorem?
        self.goals = [([],goal)]
        self.proven = False
        self.name = name
    #def intros(self): #intro_all        
    #    self.goals.append( (ctx, goal.intros())  )
    #    return self
    def equiv(self, goal2):
        ctx, goal1 = self.goals.pop()
        if prove2(Implies(And(*ctx), goal1 == goal2)):
            g = goal2
        else:
            g = goal1
        self.goals.append( (ctx, g))
        return self
    def __eq__(self,rhs):
        return self.equiv(rhs)
    #def assert(): #put new goal in stack with current context. Put into context of 1 below top
    #def assume(): #just put crap in the context.
    def intro_all(self): #name = hint maybe later
        ctx, goal = self.goals.pop()
        assert(goal.is_forall())
        vs = [FreshConst(goal.var_sort(i) , prefix=goal.var_name(i)) for i in range(goal.num_vars())]
        g = instantiate(goal,*vs) 
        self.goals.append( (ctx + vs, g)) # wait. I should keep propositions and variables seperate
        return self
    def intro_imp(self): #intro_impl
        ctx, goal = self.goals.pop()
        if is_implies(goal):
            a, b = goal.children()
            ctx.append(a)
            self.goals.append((ctx,b))
        else:
            self.goals.append((ctx,goal))
        return self
    def split(self): #z3 tactic split-clauses?
        ctx, goal = self.goals.pop()
        if is_and(goal):
            for c in goal.children():
                self.goals.append((ctx,c))
        else:
            self.goals.append((ctx,goal))
        return self
    def z3_tactic(self,t):
        t = Tactic(t)
        ctx, goal = self.goals.pop()
        #g = t(Implies(And(*ctx), goal)).as_expr()
        g = t(goal).as_expr()
        self.goals.append(([],g))
        return self
    def simpl(self):
        return self.z3_tactic("simplify")
    def congruence(self):
        #maybe search for equalities. And put them in the goal
        return self.z3_tactic("solve-eqs")
    def smt(self):
        ctx, goal = self.goals.pop()
        s = Solver()
        #s.set(**keywords)
        claim = Implies(And(*ctx), goal)
        s.add(Not(claim))
        r = s.check()
        if r  == sat:
            print("Countermodel : " + str(s.model()))
        assert(r == unsat)
        return self
    def destruct(self):
        ctx, goal = self.goals.pop()
        if is_bool(goal):
            ctx1 = ctx.copy()
            ctx2 = ctx.copy()
            ctx1.append(goal == True)
            ctx2.append(goal == False)
            self.goals.append((ctx2, BoolVal(False) ))
            self.goals.append((ctx1, BoolVal(True) ))
        else:
            self.goals.append((ctx, goal))    
        return self
    def forget(self,n):
        ctx, goal = self.goals.pop()
        ctx.pop(n)
        self.goals.append((ctx, goal))  
        return self
    def qed(self):
        if len(self.goals) == 0:
            self.proven = True
            # add self to global proof context if self.name is not None
    def get_ctx(self,n):
        return self.goals[-1][0][n]
    def __str__(self):
        if len(self.goals) >= 1:
            ctx, goal = self.goals[-1]
            return "".join([f"[{i}] {str(c)} : {str(c.sort())}\n" for i, c in enumerate(ctx)]) + "----------------\n" + f"{str(goal)} : {str(goal.sort())}"
        else:
            return "No Goals Left"
    def __repr__(self):
        return str(self)
x = Real("x")
Proof(x**2 - 1 == 0).equiv((x+1)*(x-1) == 0).equiv((x == 1) | (x == -1))
a, b = Bools('a b')
p = Proof((a & b) > b)
p.intro_imp().destruct() 
   .smt() \
   .smt() \
.qed()

Another question is how to implement an apply tactic gracefully. Fully deconstructing syntax trees and unifying ourselves is not utilizing z3 well. If you have a good idea how to get unification out of z3, I’d be interested to hear from you here. https://stackoverflow.com/questions/59398955/getting-z3-instantiations-of-quantified-variables/59400838#59400838

Here’s an idea though. In the cold light of day, I am still not sure this reasoning makes much sense. Suppose we’re trying to apply forall x. a(x) -> b(x) to a c(y). If forall x. b(x) -> c(y) we’re good and by assumption that is obvious for some reason, like the syntactic instantiation of b gives c. We can ask z3 to prove that and it will hopefully easy. If we can prove forall x. a(x) in the current context, that would be sufficient, but not true typically. It is an overly difficult request. We really only need to prove a(x) for values pertinent to the proof of c(y). Here’s a suspicious strategem. Any a -> b can be weakened to (q -> a) -> (q -> b). In particular we can choose to weaken forall x. a(x) -> b(x) to forall x. ((c(y) -> b(x)) -> a(x)) -> ((c(y) -> b(x)) -> b(x)). Then we can replace the goal with forall x. ((c(y) -> b(x)) -> a(x)) after we prove that (forall x. (c(y) -> b(x)) -> b(x)) -> c(y). Maybe c(y) -> b(x) is sufficient to restrict the values of x? Not sure.

Another rough sketch of induction on Nat. Not right yet.

def inductionNat(self):
    assert(self.num_vars() == 1 and self.var_sort(0) == IntSort() and self.is_forall())
    n = FreshInt()
    return instantiate(self, IntVal(0)) & ForAll([n],instantiate(self, n) & (n > 0) > instantiate(self, n+1))

We could also make a simple induction for ADTs based on the similar introspection we used for match above. It’s ugly but I think it works.

def induction(self):
    assert(is_quantifier(self) and self.is_forall() and self.num_vars() == 1) #we can eventually relax vars = 1
    t = self.var_sort(0)
    nc = t.num_constructors()
    th = []
    for i in range(nc):
        con = t.constructor(i)
        nfields = con.arity()
        if nfields == 0:
            th += [substitute_vars(self.body(), con())]
        else:
            hyp = []
            args = []
            for d in range(nfields):
                td = con.domain(d)
                x = FreshConst(td)
                print(x)
                if td == t:
                    hyp += [substitute_vars(self.body(), x)]
                args += [x]
            th += [ForAll(args, Implies(And(*hyp), substitute_vars(self.body(), con(*args))))]
        print(th)
    return And(*th)

I haven’t really though much about tacticals yet.

# describe_tactics() gives a list of all z3 tactics
ackermannize_bv : A tactic for performing full Ackermannization on bv instances.
subpaving : tactic for testing subpaving module.
horn : apply tactic for horn clauses.
horn-simplify : simplify horn clauses.
nlsat : (try to) solve goal using a nonlinear arithmetic solver.
qfnra-nlsat : builtin strategy for solving QF_NRA problems using only nlsat.
nlqsat : apply a NL-QSAT solver.
qe-light : apply light-weight quantifier elimination.
qe-sat : check satisfiability of quantified formulas using quantifier elimination.
qe : apply quantifier elimination.
qsat : apply a QSAT solver.
qe2 : apply a QSAT based quantifier elimination.
qe_rec : apply a QSAT based quantifier elimination recursively.
psat : (try to) solve goal using a parallel SAT solver.
sat : (try to) solve goal using a SAT solver.
sat-preprocess : Apply SAT solver preprocessing procedures (bounded resolution, Boolean constant propagation, 2-SAT, subsumption, subsumption resolution).
ctx-solver-simplify : apply solver-based contextual simplification rules.
smt : apply a SAT based SMT solver.
psmt : builtin strategy for SMT tactic in parallel.
unit-subsume-simplify : unit subsumption simplification.
aig : simplify Boolean structure using AIGs.
add-bounds : add bounds to unbounded variables (under approximation).
card2bv : convert pseudo-boolean constraints to bit-vectors.
degree-shift : try to reduce degree of polynomials (remark: :mul2power simplification is automatically applied).
diff-neq : specialized solver for integer arithmetic problems that contain only atoms of the form (<= k x) (<= x k) and (not (= (- x y) k)), where x and y are constants and k is a numeral, and all constants are bounded.
eq2bv : convert integer variables used as finite domain elements to bit-vectors.
factor : polynomial factorization.
fix-dl-var : if goal is in the difference logic fragment, then fix the variable with the most number of occurrences at 0.
fm : eliminate variables using fourier-motzkin elimination.
lia2card : introduce cardinality constraints from 0-1 integer.
lia2pb : convert bounded integer variables into a sequence of 0-1 variables.
nla2bv : convert a nonlinear arithmetic problem into a bit-vector problem, in most cases the resultant goal is an under approximation and is useul for finding models.
normalize-bounds : replace a variable x with lower bound k <= x with x' = x - k.
pb2bv : convert pseudo-boolean constraints to bit-vectors.
propagate-ineqs : propagate ineqs/bounds, remove subsumed inequalities.
purify-arith : eliminate unnecessary operators: -, /, div, mod, rem, is-int, to-int, ^, root-objects.
recover-01 : recover 0-1 variables hidden as Boolean variables.
bit-blast : reduce bit-vector expressions into SAT.
bv1-blast : reduce bit-vector expressions into bit-vectors of size 1 (notes: only equality, extract and concat are supported).
bv_bound_chk : attempts to detect inconsistencies of bounds on bv expressions.
propagate-bv-bounds : propagate bit-vector bounds by simplifying implied or contradictory bounds.
propagate-bv-bounds-new : propagate bit-vector bounds by simplifying implied or contradictory bounds.
reduce-bv-size : try to reduce bit-vector sizes using inequalities.
bvarray2uf : Rewrite bit-vector arrays into bit-vector (uninterpreted) functions.
dt2bv : eliminate finite domain data-types. Replace by bit-vectors.
elim-small-bv : eliminate small, quantified bit-vectors by expansion.
max-bv-sharing : use heuristics to maximize the sharing of bit-vector expressions such as adders and multipliers.
blast-term-ite : blast term if-then-else by hoisting them.
cofactor-term-ite : eliminate term if-the-else using cofactors.
collect-statistics : Collects various statistics.
ctx-simplify : apply contextual simplification rules.
der : destructive equality resolution.
distribute-forall : distribute forall over conjunctions.
dom-simplify : apply dominator simplification rules.
elim-term-ite : eliminate term if-then-else by adding fresh auxiliary declarations.
elim-uncnstr : eliminate application containing unconstrained variables.
injectivity : Identifies and applies injectivity axioms.
snf : put goal in skolem normal form.
nnf : put goal in negation normal form.
occf : put goal in one constraint per clause normal form (notes: fails if proof generation is enabled; only clauses are considered).
pb-preprocess : pre-process pseudo-Boolean constraints a la Davis Putnam.
propagate-values : propagate constants.
reduce-args : reduce the number of arguments of function applications, when for all occurrences of a function f the i-th is a value.
reduce-invertible : reduce invertible variable occurrences.
simplify : apply simplification rules.
elim-and : convert (and a b) into (not (or (not a) (not b))).
solve-eqs : eliminate variables by solving equations.
special-relations : detect and replace by special relations.
split-clause : split a clause in many subgoals.
symmetry-reduce : apply symmetry reduction.
tseitin-cnf : convert goal into CNF using tseitin-like encoding (note: quantifiers are ignored).
tseitin-cnf-core : convert goal into CNF using tseitin-like encoding (note: quantifiers are ignored). This tactic does not apply required simplifications to the input goal like the tseitin-cnf tactic.
qffd : builtin strategy for solving QF_FD problems.
pqffd : builtin strategy for solving QF_FD problems in parallel.
smtfd : builtin strategy for solving SMT problems by reduction to FD.
fpa2bv : convert floating point numbers to bit-vectors.
qffp : (try to) solve goal using the tactic for QF_FP.
qffpbv : (try to) solve goal using the tactic for QF_FPBV (floats+bit-vectors).
qffplra : (try to) solve goal using the tactic for QF_FPLRA.
default : default strategy used when no logic is specified.
sine-filter : eliminate premises using Sine Qua Non
qfbv-sls : (try to) solve using stochastic local search for QF_BV.
nra : builtin strategy for solving NRA problems.
qfaufbv : builtin strategy for solving QF_AUFBV problems.
qfauflia : builtin strategy for solving QF_AUFLIA problems.
qfbv : builtin strategy for solving QF_BV problems.
qfidl : builtin strategy for solving QF_IDL problems.
qflia : builtin strategy for solving QF_LIA problems.
qflra : builtin strategy for solving QF_LRA problems.
qfnia : builtin strategy for solving QF_NIA problems.
qfnra : builtin strategy for solving QF_NRA problems.
qfuf : builtin strategy for solving QF_UF problems.
qfufbv : builtin strategy for solving QF_UFBV problems.
qfufbv_ackr : A tactic for solving QF_UFBV based on Ackermannization.
ufnia : builtin strategy for solving UFNIA problems.
uflra : builtin strategy for solving UFLRA problems.
auflia : builtin strategy for solving AUFLIA problems.
auflira : builtin strategy for solving AUFLIRA problems.
aufnira : builtin strategy for solving AUFNIRA problems.
lra : builtin strategy for solving LRA problems.
lia : builtin strategy for solving LIA problems.
lira : builtin strategy for solving LIRA problems.
skip : do nothing tactic.
fail : always fail tactic.
fail-if-undecided : fail if goal is undecided.
macro-finder : Identifies and applies macros.
quasi-macros : Identifies and applies quasi-macros.
ufbv-rewriter : Applies UFBV-specific rewriting rules, mainly demodulation.
bv : builtin strategy for solving BV problems (with quantifiers).
ufbv : builtin strategy for solving UFBV problems (with quantifiers).

Proving some Inductive Facts about Lists using Z3 python

Z3 is an SMT solver which has a good rep. Here are some excellent tutorials.

https://rise4fun.com/z3/tutorial

https://theory.stanford.edu/~nikolaj/programmingz3.html

http://homepage.cs.uiowa.edu/~ajreynol/VTSA2017/

SMT stands for satisfiability modulo theories. The exact nature of power of these kinds of solvers has been and is still hazy to me. I have known for a long time that they can slam sudoku or picross or other puzzles, but what about more infinite or logic looking things? I think I may always be hazy, as one can learn and discover more and more encoding tricks to get problems and proofs that you previously thought weren’t solvable into the system. It’s very similar to learning how to encode to linear programming solvers in that respect.

SMT solvers combine a general smart search procedure with a ton of specialized solvers for particular domains, like linear algebra, polynomials, linear inequalities and more.

The search procedure goes by the name DPLL(T). It is an adjustment of the procedure of SAT solvers, which are very powerful and fast. SAT solvers find an assignment of boolean variables in order to make a complicated boolean expression true, or to eventually find that it can never be made true. SAT solvers work on the principle of guessing and deducing. If a OR b needs to be true and we already know a is false, we can deduce b must be true. When the deduction process stalls, we just guess and then backtrack if it doesn’t work out. This is the same process you use manually in solving Sudoku.

The modern era of SAT solvers came when a couple new tricks vastly extended their power. In particular Conflict Driven Clause Learning (CDCL), where when the solver finds itself in a dead end, it figures out the choices it made that put it in the dead end and adds a clause to its formula so that it will never make those choices again.

https://sahandsaba.com/understanding-sat-by-implementing-a-simple-sat-solver-in-python.html

SMT works by now having the boolean variables of the SAT solver contain inner structure, like the boolean p actually represents the fact x + y < 5. During the search process it can take the pile of booleans that have been set to true and ask a solver (maybe a linear programming solver in this case) whether those facts can all be made true in the underlying theory. This is an extra spice on top of the SAT search.

Something that wasn’t apparent to me at first is how important the theory of uninterpreted formulas is to SMT solvers. It really is their bread and butter. This theory is basically solved by unification, which is the fairly intuitive process of finding assignments to variables to make a set of equations true. If I ask how to make fred(x,4) = fred(7,y), obviously the answer is y=4, x=7. That is unification. Unification is a completely syntax driven way to deduce facts. It starts to give you something quite similar to first order logic.

https://eli.thegreenplace.net/2018/unification/

https://cs.mtsu.edu/~rbutler/courses/sdd/automated_deduction/unification.pdf

I was also under the impression that quantifiers \forall, \exists were available but heavily frowned upon. I don’t think this is quite true. I think they are sort of a part of the entire point of the SMT solver now, although yes, they are rather flaky. There are a number of methods to deal with the quantifier, but one common one is to look for a pattern or parts of the subformula, and instantiate a new set of free variables for all of the quantified ones and add the theorem every time the patterns match. This is called E-matching.

Here are a couple tutorials on proving inductive facts in SMT solvers. You kind of have to hold their hand a bit.

http://lara.epfl.ch/~reynolds/vmcai15.pdf

http://homepage.divms.uiowa.edu/~ajreynol/pres-ind15.pdf

SMT solvers queries usually have the flavor of finding something, in this case a counterexample. The idea is that you try to ask for the first counterexample where induction failed. Assuming that proposition P was true for (n-1), find n such that P is not true. If you can’t find it, then the induction has gone through.

And here is some code where I believe I’m showing that some simple list functions like reverse, append, and length have some simple properties like \forall t. rev (rev(t)) == t .

from z3 import *

# Very useful reference
# https://theory.stanford.edu/~nikolaj/programmingz3.html

f = Function('f', IntSort(), IntSort())
s = Solver()
#s.add(f(True) == False, f(False) == True)
x = Int('x')
s.add(ForAll([x], f(x) >= x)) #> and < do not seem to be returning
print(s.sexpr())
s.check()
print(s.model())

# Rolling my own list data type
# Z3 has a built in which will probably be better?
s = Solver()

L = Datatype("MyList")
L.declare("Nil")
L.declare("Cons", ("hd", IntSort()), ("tail", L))
L = L.create()

t = Const('t', L)
u = Const('u', L)

y = Int('y')


rev = Function('reverse', L, L)
app = Function('append', L, L, L)
leng = Function('length', L, IntSort())

#defining my functions. Micro Haskell, BABY
#length
s.add( leng(L.Nil) == 0 )
s.add(  ForAll([u,y],  leng(L.Cons(y,u)) == 1 + leng(u))) #  patterns = [leng(L.Cons(y,u))] 

#append
s.add(  ForAll([u], app(L.Nil, u) == u)) 
s.add(  ForAll([t, u, y] , app(L.Cons(y,t), u)  ==  L.Cons(y, app(t, u ))))

#reverse
s.add( rev(L.Nil) == L.Nil)
s.add(  ForAll([y,t],  rev(L.Cons(y,t)) == app(rev(t), L.Cons(y, L.Nil))))


s.push()
print("proving leng(t) >= 0")
#s.add( Or(And(t == L.Cons(y,u),  leng(u) >= 0 ), t == L.Nil))
s.add(  Not(leng(t) >= 0 ))
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),  leng(L.tail(t)) >= 0 ))
print(s.check())

s.pop()
s.push()
#s.add( leng(app(L.Nil, u)) == leng(u) )

print("prove length is preserved under app.")
s.add( leng(app(t,u)) != leng(t) + leng(u))
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   leng(app(L.tail(t),u)) == leng(L.tail(t)) + leng(u)  ))
print(s.check())

s.pop()
s.push()

print("reverse preserves length")
#Lemma Could place in clause with the above proof of this theorem
s.add( ForAll([t,u], leng(app(t,u)) == leng(t) + leng(u)   )) #subgoal needed
s.add( leng(rev(t)) != leng(t))
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   leng(rev(L.tail(t))) == leng(L.tail(t)) ))


s.pop()
s.push()

print("reverse reverse = id")
s.add( ForAll( [t,u], rev(app(t,u)) ==  app(rev(u), rev(t)) ) ) #subgoal needed
s.add( rev(rev(t)) != t )
s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   rev(rev(L.tail(t))) == L.tail(t) ))
print(s.check())


#junk

#s.add(t != L.Nil ) 
#s.add( ForAll([t], rev(L.Cons(y,t)) ==   ) , rev(L.Nil) == L.Nil)
#s.add( leng(L.Cons()) == 0 )
#s.add(  ForAll(y,  leng(L.Cons(y,L.Nil)) == 1 + leng(L.Nil)))
#s.add(  Not( leng(app(t,u))  == leng(t) + leng(u) )) 

# prove length of app + nil is same
#s.add( leng(app(t,L.Nil)) != leng(t))
#s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   leng(app(L.tail(t),L.Nil)) == leng(L.tail(t))))

#s.add( app(L.Nil,L.Nil) == L.Nil)
#s.add( app(t,L.Nil) != t)
#s.add( Implies(t == L.Cons(L.hd(t), L.tail(t)),   app(L.tail(t),L.Nil) == L.tail(t)))


#s.add(Or(t == L.Nil , And( t == L.Cons(y, u),   app(u,L.Nil) == u)  ))
#s.add( Implies(t == L.Nil, leng(app(L.Nil, u)) == leng(u) ))
# all of these don't work
#s.add( Implies(u == L.tail(t),  leng(u) >= 0 ))
#s.add( ForAll(u, Implies(t == L.Cons(y,u),  leng(u) >= 0 )))
#print(s.sexpr())
#print(s.check())
#print(s.model())

#print(tactics())
'''
def induction(freevar, f, construtors?):
    x = Int('x')
    return And(Not(f(x)), 
'''

'''
# controller synthesis
pos = Array(Real, 10) 
for t in range(10):
    s.add(pos[t] == pos[t] +  v * dt)
    s.add(v[t] == v[t] + f(pos[t]) * dt)
    s.add(f(pos[t] <= 1)) # constraints on force
s.add(ForAll([init], Implies(And(init <= 1, pos[0] == init), pos[9] <= 0.1) ))
 
s.add(Forall[x], f(x) == a * x) # parametrizing. 
'''
#s.set("produce-proofs", True
#s.set("mbqi", False)
#s.set(mbqi=False)
#print(s.proof())

Gettin’ that Robot some Tasty Apples: Solving a simple geometrical puzzle in Z3 python

At work there is a monthly puzzler.

“Design a robot that can pick up all 9 apples arranged on a 3 by 3 rectangular grid, and spaced 1m apart. The robot parts they have are faulty. The robot can only turn three times”

I think the intent of the puzzle is that the robot can drive in forward and reverse, but only actually turn 3 times. It’s not very hard to do by hand. I decided to take a crack at this one using Z3 for funzies. Z3 is an SMT solver. It is capable of solving a interesting wide variety of problems.

I interpret this as “find 4 lines that touch all points in the grid, such that each subsequent line intersects.”

It is fairly easy to directly translate this into a Z3 model.

from z3 import *
import matplotlib.pyplot as plt
import numpy as np
s = Solver()

linenum = 4

def linepoint(l,p): #line point product. Is zero if point is on line
   return l[0]*p[0] + l[1]*p[1] + l[2]

lines = [[Real(s + str(i)) for s in ('a','b','c')] for i in range(linenum)]

applepoints = [(x,y) for x in range(5) for y in range(3)]
turnpoints = [[Real(s + str(i)) for s in ('x','y')] for i in range(linenum-1)]

for pt in applepoints:  #Every point needs to be touched by one line
	s.add(Or([linepoint(l,pt)==0 for l in lines]))

for l in lines: #Exclude invalid lines (all zeros)
		s.add(Or([a != 0 for a in l]))

for i in range(linenum-1): #Consecutive lines intersect (aren't parallel). There are other ways of doing this
	s.add(linepoint( lines[i], turnpoints[i]) == 0)
	s.add(linepoint( lines[i+1], turnpoints[i]) == 0)









def convert(x):
	if is_int_value(x):
		return x.as_long()
	elif is_rational_value(x):
		return x.numerator_as_long()/x.denominator_as_long()
	elif is_algebraic_value(x):
		return convert(x.approx(7))


print(s.check())

m = s.model()
#print("x = %s" % m[x])

#print "traversing model..."
for d in m.decls():
    print("%s = %s" % (d.name(), m[d]))

#nplines = np.array([[m[a].numerator_as_long()/m[a].denominator_as_long() for a in l] for l in lines])
nplines = np.array([[convert(m[a]) for a in l] for l in lines])

#xs = np.random.randn(2,3)
#xs = np.




print(nplines)
for l in nplines:
	pts = []
	for j in [-1,1]:
		# We can use Z3 to draw a nice set of lines too
		opt = Optimize()
		x, y = Reals('x y')
		opt.add([0 <=  x, x <= 2, 0 <= y, y <= 2])
		opt.add(l[0] * x + l[1]*y + l[2] == 0)
		opt.minimize(j*0.453 * x + j*0.3234 * y) # random objective direction hopefully not perpendicular to a line
		print(opt.check())
		m = opt.model()
		pts.append([convert(m[x]), convert(m[y])])
		print(l)
	pts = np.array(pts).T
	plt.plot(pts[0,:], pts[1,:])
plt.show()

A couple comments:

If we ask z3 to use only 3 lines, it returns unsat. Try to prove that by hand.

However, If the robot is on the projective plane, it is possible with 3 lines. It only needs to drive to infinity and turn twice. All lines intersect exactly once on the projective plane. How convenient.

The problem only seems somewhat difficult to computerize because of the seemingly infinite nature of geometry. If we only consider the lines that touch at least two points, all possible robot paths becomes extremely enumerable. Is there a proof that we only need these lines?

Another interesting approach might be to note that the points are described by the set of equations x*(x-1)*(x-2)=0 and y*(y-1)*(y-2)=0. I think we could then possibly use methods of nonlinear algebra (Groebner bases) to find the lines. Roughly an ideal containment question? Don’t have this one fully thought out yet. I think z3 might be doing something like this behind the scenes.