Fixedpoint：用递归的规则推导可达性

Fixedpoint：用递归的规则推导可达性#

本章的内容#

Z3 的 Fixedpoint 是处理递归规则的引擎，像从 A 能到 B、从 B 能到 C 则从 A 能到 C 这样的推理，可以直接写成规则，默认的引擎是自底向上的 Datalog。

  • 用 rule 写出规则，用 fact 写出事实

  • 用 query() 询问是否成立

  • 到达之后取出实际路径的办法

Datalog 与关系数据库的查询，即 SQL 的递归查询，很相似，从给定的事实出发反复应用规则，不断增加新的事实。这个过程会一直持续到没有新的事实可加为止，因此总能停机。

from z3 import *

例 1：图的连通性#

下面这个无向图中，考察顶点 1 能否到达顶点 6。

        flowchart LR
1 --- 2
1 --- 5
5 --- 2
4 --- 5
2 --- 3
4 --- 3
4 --- 6

用到三个谓词。

  • edge(a, b)：a 与 b 之间有边，作为事实给出

  • connect(a, b)：a 可以直接移动到 b，只要有一条边就双向成立

  • path(a, b)：a 可以到达 b，由递归确定

connect 由 edge 推导出来，两个方向各写一条规则，因此无向图只需要存一份边。

fp = Fixedpoint()
fp.set(engine='datalog')          # 使用 Datalog 引擎

s = BitVecSort(8)                  # 顶点编号用 8 位表示

connect = Function('connect', s, s, BoolSort())
edge = Function('edge', s, s, BoolSort())
path = Function('path', s, s, BoolSort())

a, b, c = Consts('a b c', s)

fp.register_relation(path, edge, connect)
fp.declare_var(a, b, c)

# 规则：有边就能直接移动（无向，所以两个方向都成立）
fp.rule(connect(a, b), edge(a, b))
fp.rule(connect(a, b), edge(b, a))

# 规则：能直接移动，就能到达
fp.rule(path(a, b), connect(a, b))
# 规则：a→b 能移动，且 b→c 能到达，则 a→c 能到达
fp.rule(path(a, c), [connect(a, b), path(b, c)])

# 事实：图的边
edges = [(1, 2), (1, 5), (5, 2), (4, 5), (2, 3), (4, 3), (4, 6)]
for p, q in edges:
    fp.fact(edge(p, q))

代码的说明#

1. 引擎

fp = Fixedpoint()
fp.set(engine='datalog')

engine='datalog' 选择 Datalog 引擎。

2. 顶点的类型

s = BitVecSort(8)

顶点的编号用 8 位的位向量表示。Datalog 只接受有限的类型，IntSort 这类无限的型用不了，因此这里选 BitVecSort：8 位只有 256 个取值，范围有限。

3. 谓词与关系的定义

connect = Function('connect', s, s, BoolSort())

用 Function(名字, 参数的类型..., 返回值的类型) 造出谓词。这是接收两个顶点、返回真假的二元关系，因此返回值取 BoolSort()。三个谓词都是二元关系，类型与元数完全一样。

4. 变量的声明与谓词的注册

fp.register_relation(path, edge, connect)
fp.declare_var(a, b, c)

用 register_relation() 把谓词注册到引擎，用 declare_var() 声明规则中使用的变量。注册的顺序不影响结果，漏注册任何一个都会让规则里的谓词被当作未定义。

5. 规则与事实

fp.rule(path(a, c), [connect(a, b), path(b, c)])

rule(头, 本体) 表示本体全部成立时头部也成立。这条规则读作 a 能到 b、b 能到 c，则 a 能到 c，把 path 自身放进本体就是递归。

fp.fact(edge(1, 2)) 把存在一条边登记为事实。

询问#

print(fp.query(path(1, 6)))     # 能否到达
print(fp.query(path(1, 7)))     # 不存在的顶点

sat
unsat

顶点 1 到 6 可以到达，即 sat，而顶点 7 不可到达，即 unsat。回答这个是否成立的问题，就是 Fixedpoint 的工作。query() 只给出真或假，不给路径也不给证明，与 Solver.check() 的做法相似。

重建路径#

知道能到达之后，还不知道具体的走法。

Z3 的 Datalog 有输出说明的功能，即 datalog.generate_explanations，可以说明某个事实为何成立。不过这个功能在较新的 Z3，即这个环境中的 5.1.0 里并不稳定，会出现 explanations are not supported with undefined predicates 这样的错误，或者答案变成 unsat，甚至 Z3 本身崩溃。

因此本章改为自己构造路径：能否到达问 Fixedpoint，具体路径则对边的列表做 BFS，即广度优先搜索来求。

from collections import deque


def bfs_route(edges, start, goal, undirected=True):
    """从 (from, to) 的列表出发，用 BFS 构造最短路径"""
    neighbors = {}
    for p, q in edges:
        neighbors.setdefault(p, []).append(q)
        if undirected:
            neighbors.setdefault(q, []).append(p)

    prev = {start: None}
    queue = deque([start])
    while queue:
        cur = queue.popleft()
        if cur == goal:
            break
        for nxt in neighbors.get(cur, []):
            if nxt not in prev:
                prev[nxt] = cur
                queue.append(nxt)

    if goal not in prev:
        return None

    route = [goal]
    while prev[route[-1]] is not None:
        route.append(prev[route[-1]])
    return route[::-1]


print(bfs_route(edges, 1, 6))

[1, 5, 4, 6]

路径是 1 → 5 → 4 → 6，由于 edge(4, 5) 是无向边，5 到 4 这个方向也可以使用。bfs_route() 只按边列表找最短路，与 query() 的结果互为印证。

例 2：三个水桶的问题#

有 8 升、5 升、3 升三个水桶，最初只有 8 升的桶是满的，其余都是空的，要让两个桶各装 4 升，应当如何倒水。

这个问题可以看作以状态为顶点、以倒水操作为边的图的连通性，先把可能出现的状态与状态之间的迁移全部用 Python 造出来。

  • 水的总量始终是 8 升

  • 迁移只有把某个桶里的水倒进另一个桶

  • 倒水之后，要么目标桶被倒满，要么源桶被倒空

from itertools import product

sizes = (8, 5, 3)

# 可能的状态（合计 8 升）
states = [(a, b, c)
          for a, b, c in product(*[range(size + 1) for size in sizes])
          if a + b + c == 8]

# 一次倒水能到达的状态对
transitions = []
for s1, s2 in product(states, states):
    changed = [(size, old, new) for size, old, new in zip(sizes, s1, s2) if old != new]
    if len(changed) != 2:
        continue        # 变化的桶正好 2 个（倒出来源和倒进目标）
    if any(new == size or new == 0 for size, old, new in changed):
        transitions.append((s1, s2))    # 目标桶满 or 来源桶空

print(len(states), "种状态，", len(transitions), "种倒法")

24 种状态， 106 种倒法

接着把这种迁移作为 pour 谓词的事实登记，再用 plan 谓词推出从初始状态到目标状态是否存在步骤。状态用三个桶的水量表示，因此 pour 与 plan 都是六元的谓词。

  • pour(a, b, c, a2, b2, c2)：从状态 (a, b, c) 一次就能移到 (a2, b2, c2)

  • plan(a, b, c, a2, b2, c2)：从状态 (a, b, c) 经过若干次可以移到 (a2, b2, c2)

fp = Fixedpoint()
fp.set(engine='datalog')

s = BitVecSort(5)
pour = Function('pour', s, s, s, s, s, s, BoolSort())
plan = Function('plan', s, s, s, s, s, s, BoolSort())
a, b, c, a1, b1, c1, a2, b2, c2 = Consts(['a', 'b', 'c', 'a1', 'b1', 'c1', 'a2', 'b2', 'c2'], s)

fp.register_relation(pour, plan)
fp.declare_var(a, b, c, a1, b1, c1, a2, b2, c2)

# 一步能移，则 plan 也成立
fp.rule(plan(a, b, c, a2, b2, c2), pour(a, b, c, a2, b2, c2))
# 移一步之后还有 plan，则 plan 成立（递归）
fp.rule(plan(a, b, c, a2, b2, c2), [pour(a, b, c, a1, b1, c1), plan(a1, b1, c1, a2, b2, c2)])

for s1, s2 in transitions:
    fp.fact(pour(*s1, *s2))

start, goal = (8, 0, 0), (4, 4, 0)

%time print(fp.query(plan(*start, *goal)))

sat
CPU times: total: 0 ns
Wall time: 1.78 ms

结果是 sat，因此从 8,0,0 到 4,4,0 存在一组步骤。plan 只回答是否存在这样的步骤序列，不回答有几步。

步骤可以对迁移的列表 transitions 做 BFS 求出，例 1 的 bfs_route() 可以直接使用：例 1 的图无向，两个方向都要试，而倒水有先后，倒回去的那一步不算，因此这里要传 undirected=False。

route = bfs_route(transitions, start, goal, undirected=False)

for s1, s2 in zip(route, route[1:]):
    print(s1, "->", s2)

(8, 0, 0) -> (3, 5, 0)
(3, 5, 0) -> (3, 2, 3)
(3, 2, 3) -> (6, 2, 0)
(6, 2, 0) -> (6, 0, 2)
(6, 0, 2) -> (1, 5, 2)
(1, 5, 2) -> (1, 4, 3)
(1, 4, 3) -> (4, 4, 0)

七次倒水之后到达了 4,4,0，人手工解这类问题需要反复试探，而只要把状态与迁移写出来，剩下的部分就可以机械地求解。可能的状态总数是有限的：三个桶的容量确定之后，组合不超过 9×6×4 种。

登记过的规则与事实也可以查看。get_rules() 返回的正是登记时的写法，可以用来核对规则是否有遗漏。

for rule in fp.get_rules()[:5]:
    print(rule)
print("...")

ForAll([A, B, C, D, E, F],
       Implies(pour(A, B, C, D, E, F),
               plan(A, B, C, D, E, F)))
ForAll([A, B, C, D, E, F, G, H, I],
       Implies(And(plan(D, E, F, G, H, I),
                   pour(A, B, C, D, E, F)),
               plan(A, B, C, G, H, I)))
pour(0, 5, 3, 3, 5, 0)
pour(0, 5, 3, 5, 0, 3)
pour(1, 4, 3, 0, 5, 3)
...

例 3：把 pour 写成规则#

例 2 中是用 Python 算出全部迁移，再把 pour 作为事实登记。这一回把 pour 写成规则：把可能出现的状态 state(a, b, c) 作为事实给出，把倒水的条件写成规则。两者的区别在于谁来枚举：前者由 Python 枚举，后者由引擎在推导时枚举。

条件是变化的桶只有两个，并且其中一个是满的或空的。写成规则时，这个条件展开成三条分支，每条对应一个桶作为源桶。

fp = Fixedpoint()
fp.set(engine='datalog')

s = BitVecSort(5)
state = Function('state', s, s, s, BoolSort())
pour = Function('pour', s, s, s, s, s, s, BoolSort())
plan = Function('plan', s, s, s, s, s, s, BoolSort())
a, b, c, a1, b1, c1, a2, b2, c2 = Consts(['a', 'b', 'c', 'a1', 'b1', 'c1', 'a2', 'b2', 'c2'], s)

fp.register_relation(state, pour, plan)
fp.declare_var(a, b, c, a1, b1, c1, a2, b2, c2)

fp.rule(plan(a, b, c, a2, b2, c2), pour(a, b, c, a2, b2, c2))
fp.rule(plan(a, b, c, a2, b2, c2), [pour(a, b, c, a1, b1, c1), plan(a1, b1, c1, a2, b2, c2)])

# 倒水规则：变化的桶有 2 个，其中一个变满或变空
fp.rule(
    pour(a, b, c, a1, b1, c1),
    [
        state(a, b, c),
        state(a1, b1, c1),
        Or(
            And(a == a1, b != b1, Or(b1 == 0, b1 == 5, c1 == 0, c1 == 3)),
            And(b == b1, a != a1, Or(a1 == 0, a1 == 8, c1 == 0, c1 == 3)),
            And(c == c1, a != a1, Or(a1 == 0, a1 == 8, b1 == 0, b1 == 5)),
        ),
    ],
)

# 把可能的状态（合计 8 升）注册为事实
for p in states:
    fp.fact(state(*p))

%time print(fp.query(plan(*start, *goal)))

sat
CPU times: total: 62.5 ms
Wall time: 47.1 ms

把 pour 写成规则、再把状态作为事实给出，结果与例 2 相同，都是 sat，在这个环境里约 0.04 秒，而例 2 约 1 毫秒。

读 pour 的规则时要注意，state(a, b, c) 与 state(a1, b1, c1) 两者都放进了本体，倒水之后的状态也必须是一个可能出现的状态。少了后一个条件，规则就会推出并不存在的状态。

例 4：把 state 也写成规则#

最后把 state 也写成规则，即把每个桶可能的水量作为 bucket_a、bucket_b、bucket_c 的事实给出，再由它们的组合，合计 8 升，推导出 state。

fp = Fixedpoint()
fp.set(engine='datalog')

s = BitVecSort(5)
bucket_a = Function('bucket_a', s, BoolSort())
bucket_b = Function('bucket_b', s, BoolSort())
bucket_c = Function('bucket_c', s, BoolSort())
state = Function('state', s, s, s, BoolSort())
pour = Function('pour', s, s, s, s, s, s, BoolSort())
plan = Function('plan', s, s, s, s, s, s, BoolSort())
a, b, c, a1, b1, c1, a2, b2, c2 = Consts(['a', 'b', 'c', 'a1', 'b1', 'c1', 'a2', 'b2', 'c2'], s)

fp.register_relation(state, pour, plan, bucket_a, bucket_b, bucket_c)
fp.declare_var(a, b, c, a1, b1, c1, a2, b2, c2)

fp.rule(plan(a, b, c, a2, b2, c2), pour(a, b, c, a2, b2, c2))
fp.rule(plan(a, b, c, a2, b2, c2), [pour(a, b, c, a1, b1, c1), plan(a1, b1, c1, a2, b2, c2)])

fp.rule(
    pour(a, b, c, a1, b1, c1),
    [
        state(a, b, c),
        state(a1, b1, c1),
        Or(
            And(a == a1, b != b1, Or(b1 == 0, b1 == 5, c1 == 0, c1 == 3)),
            And(b == b1, a != a1, Or(a1 == 0, a1 == 8, c1 == 0, c1 == 3)),
            And(c == c1, a != a1, Or(a1 == 0, a1 == 8, b1 == 0, b1 == 5)),
        ),
    ],
)

# state 是“各桶水量的组合”，且合计为 8 升
fp.rule(state(a, b, c), [bucket_a(a), bucket_b(b), bucket_c(c), a + b + c == 8])

for i in range(sizes[0] + 1):
    fp.fact(bucket_a(i), str(i))
for i in range(sizes[1] + 1):
    fp.fact(bucket_b(i), str(i))
for i in range(sizes[2] + 1):
    fp.fact(bucket_c(i), str(i))

%time print(fp.query(plan(*start, *goal)))

sat
CPU times: total: 844 ms
Wall time: 832 ms

本例也是 sat，但耗时明显增加，在这个环境里是 0.6 秒，与例 3 的 0.04 秒相比约 15 倍。

把 state 写成规则之后，每次应用 pour 的规则都要重新考虑 state 的组合，探索的量就随之增加。规则越抽象，引擎在推导时要枚举的组合就越多。能写成规则与跑得快是两件事，本例说明已知的关系作为事实给出会更快。

小结#

使用的写法

作用

Fixedpoint()

建立引擎

set(engine='datalog')

选择 Datalog 引擎

Function(名字, 类型..., BoolSort())

造出谓词

register_relation(...)

把谓词注册到引擎

declare_var(...)

声明规则中使用的变量

rule(头, 本体)

追加本体成立则头部成立的规则

fact(...)

追加事实

query(...)

询问是否成立，返回 sat 或 unsat

get_rules()

查看已登记的规则与事实

  • 可达性可以用一次就能移动与递归这两条规则写出

  • Fixedpoint 回答的只是能否到达，具体路径对边的列表做 BFS 即可得到

  • 用说明功能 generate_explanations 取路径在较新的 Z3 中并不稳定，不宜依赖

  • 状态与迁移已知时，作为事实给出更快，全部写成规则虽然灵活却更慢