# 实验 04：位图索引与 Bloom 过滤器 —— 把 WHERE 子句变成位运算

> 用一张 50 万行的地址表，测出"逐行判断"和"位图 AND"的差距，再压缩位图、再把空间压到极限变成概率结构。

## 一、实验目标

完成本实验后，你应当能够：

1. 手工构造等值编码（equality encoding）位图索引，并用位运算回答组合条件的查询。
2. 解释同一份位图用"布尔列表"和"Python 大整数"两种表示时，为什么查询性能差一个量级。
3. 说清楚范围编码（range encoding）与等值编码的差别，知道什么时候该用哪一个。
4. 实现 WAH（Word-Aligned Hybrid）字对齐游程压缩，量化"数据越稀疏、越聚簇压缩越划算"。
5. 实现 Bloom 过滤器，实测假阳性率并与理论公式对比，理解它多出两个可控旋钮：位数 m 与哈希个数 k。
6. 理解三者的取舍：精确位图快但费空间，WAH 省空间，Bloom 最省但有假阳性、需要回盘复核。

## 二、预备知识

### 2.1 位序列

位序列（bit sequence）是长度为 N 的 0/1 串，第 j 位表示"第 j 行是否满足某个条件"。数据库里一张 N 行的表，一个谓词（predicate）就能生成一条位序列。

```
行号 j      0  1  2  3  4  5  6  7
city='NY'   1  0  0  1  0  1  0  1
```

这就是"第 0、3、5、7 行在纽约"。

位序列的关键性质：**逻辑运算可以按位并行做**。

| SQL 逻辑 | 位运算 | Python 写法 |
| --- | --- | --- |
| `p AND q` | 按位与 AND | `a & b` |
| `p OR q` | 按位或 OR | `a \| b` |
| `NOT p` | 按位取反 NOT | `~a` |
| `p AND NOT q` | AND NOT | `a & ~b` |

### 2.2 位序列的两种物理表示

| 表示方式 | 每个 0/1 占多少空间 | AND 一次要做多少事 | 备注 |
| --- | --- | --- | --- |
| 布尔列表 `list[bool]` | Python 里每个元素是对象引用，实际远大于 1 位（约 8 字节指针 + 单例对象） | N 次 Python 层循环 + N 次布尔判断 | 直观但极慢 |
| Python 大整数 `int` | 1 位（紧凑二进制位，连续的机器字） | 一次 CPU 位运算指令，覆盖全部 N 位 | `bin(x).count("1")` 统计置位数 |

第二种是数据库里的做法：把 N 行压缩成 N 个连续的比特，CPU 的一条 `AND` 指令一次处理 64 位，Python 的大整数运算在幕后把这件事扩展到任意长度。

### 2.3 等值编码与范围编码

设属性 A 有取值 v₁…v_d，表有 N 行。

**等值编码（equality encoding）**：每个取值一条位序列。

```
E_i[j] = 1  ⟺  A[j] = v_i
```

查 `A = v_i` 只需取第 i 条；查 `A ∈ {v_i, v_k}` 是 `E_i | E_k`；查 `A ≤ v_i` 需要把 ≤ i 的所有位序列都 OR 起来。

**范围编码（range encoding）**：每个取值一条位序列，存的是"小于等于"。

```
R_i[j] = 1  ⟺  A[j] ≤ v_i
```

查 `A ≤ v_i` 只需一次取 `R_i`；查 `A = v_i` 变成 `R_i & ~R_{i-1}`。

两种方案谁更好取决于查询混合比，系统的比较可以参考 Chan 与 Ioannidis 1998 年发表的位图索引比较工作。

### 2.4 WAH 压缩

WAH（Word-Aligned Hybrid）是一种对位图做游程压缩（run-length encoding）的方案，原始论文是 Wu、Otoo、Shoshani 发表于 ACM TODS 2006 的 *Optimizing Bitmap Indices with Efficient Compression*。

做法是把位序列按机器字长（本实验 32 位）切成若干字（word）：

```
字类型        含义                              存储开销
──────────────────────────────────────────────────────
0-fill       连续若干个字全是 0                1 个字（头里记录重复次数）
1-fill       连续若干个字全是 1                1 个字
literal      字内 0/1 混杂，无法压缩            1 个头 + 该字原样存
```

```
原始:  00000000...0000 | 0010...1001 | 1111...1111 | 0000...0000
        └─ 30 个字全 0 ─┘ └ 混杂，原样 ┘ └ 20 个全 1 ┘ └ 剩余全 0 ┘
WAH:   [0-fill, 30]       [literal ...]  [1-fill, 20]  [0-fill, ...]
       1 个字              2 个字         1 个字        1 个字
```

图 1：WAH 的编码方式。字对齐的好处是某些运算可以直接在压缩态上做，不必解压。

另一种主流实现是 Roaring Bitmap：把 32 位整数按高 16 位分到 65536 个桶，每个桶内根据稀疏程度选用 array（16 位整数数组）、bitset（1024 个 64 位字）或 run（游程）三种容器之一。它和 WAH 是同一问题的两条技术路线。

### 2.5 Bloom 过滤器

Bloom 过滤器（Bloom filter）是 Bloom 于 1970 年提出的近似集合成员判定结构。它的目标不是"回答准确"，而是"用极小空间挡掉绝大多数不可能命中的查询"。

结构：长度为 m 的位数组 + k 个哈希函数。

- 插入键 x：算出 k 个位置 h₁(x)…h_k(x)，把这些位置置 1。
- 查询键 y：算出 k 个位置，若其中有任何一个为 0 → **一定不在**（无假阴性）；若全为 1 → **可能在**，需要回原始数据复核。

```
插入 x        h1(x)    h2(x)    h3(x)
              ↓        ↓        ↓
位数组  ... 0  1  0  0  1  0  0  1  0 ...
                                    ↑ 这些位置可能被别的键碰出来

查询 y （y 没插过，但碰巧落在三个已置 1 的位置上）→ 判定"可能存在"
                                        ↑ 这就是假阳性（false positive）
```

假阳性率：

```
p ≈ (1 − e^{−kn/m})^k
```

符号含义：

- `n`：插入的元素个数（个）
- `m`：Bloom 过滤器的位数（bit）
- `k`：哈希函数个数（个）
- `e^{−kn/m}`：插入全部 n 个元素后，某个特定位仍为 0 的概率
- `p`：一次"不存在键"查询被误判为存在的概率

最优哈希函数个数：

```
k* = (m/n) · ln 2
```

给定目标假阳性率 p 时所需位数：

```
m = − n · ln p / (ln 2)²
```

应用：Google Bigtable（OSDI'06）用它避免读不存在的 SSTable；Apache Cassandra 用它减少磁盘读取；PostgreSQL 有 bloom 扩展索引，属于有损索引，命中后必须回表复核。

![位图运算示意](../figures/fig_bitmap_ops.png)

图 2：一个简单的例子，把三条谓词位图用 OR、AND 组合成结果。

![位图索引性能](../figures/fig_bitmap_bench.png)

图 3：50 万行表上的真实测量。左图是查询耗时对比，右图把"建索引"也放进来，说明索引的代价要靠多次查询摊薄。

## 三、实验原理

### 3.1 为什么位图查询能快一个量级

朴素逐行扫描要做 N 次 Python 解释层循环，每次做若干次比较和布尔判断。位图方案把这些全部下推到 CPU 的位运算上：`city[0] & (street[1] | street[2])` 这一句里，Python 只是发出一条指令，真正干活的是 C 层的大整数运算代码，一次覆盖 500,000 位（62.5 KiB）。统计命中数用 `bin(x).count("1")`（Python 3.8 之后可用 `int.bit_count()`）。

代价有两个：建索引要花钱（把 N 行扫一遍位操作）；结果是有位点还要还原成行号（此处只统计数量，跳过了这一步，真实系统里要额外做一次位扫描）。

### 3.2 索引要多少次查询才回本

设：

- `C_build`：构建索引耗时
- `C_naive`：一次朴素查询耗时
- `C_bm`：一次位图查询耗时

则查询次数 Q 达到下面的值之后，索引才"赚回来"：

```
Q ≥ C_build / (C_naive − C_bm)
```

代入本实验的参考实测值，后面会算出一个具体的数字。这是数据库决定是否建索引的核心算式：写多读少的表不该建索引。

### 3.3 空间的三档取舍

| 方案 | 每行的空间开销 | 是否精确 | 运算能否在压缩态上做 |
| --- | --- | --- | --- |
| 精确位图（大整数） | 1 bit × 取值个数 | 精确 | 天然支持 |
| WAH 压缩 | 取决于数据（稀疏时远小于 1 bit） | 精确 | 部分支持（字对齐设计的目的） |
| Bloom 过滤器 | m/n bit（通常 4–16） | 概率的（有假阳性） | 只支持成员判定 |

三者不是替代关系：位图索引用在列上（取值个数少的低基数列），Bloom 用作某个数据块或文件的摘要，回答"这个块里到底有没有这个键"。

## 四、实验环境与准备

### 4.1 软件要求

- Python 3.10 及以上
- 必装：`numpy`、`matplotlib`
- 选装：`pandas`（只用于把输出表格化，不装也能跑）

```bash
pip install numpy matplotlib pandas
```

### 4.2 目录约定

本指导书位于 `…/教学产出/03-实验指导书/实验04-位图与Bloom过滤器.md`，图件在 `…/教学产出/figures/`，用相对路径 `../figures/…` 引用。

### 4.3 两种学法

- **A 类（必做，自包含）**：只用 `numpy` + 标准库。
- **B 类（选做）**：调用 ExplainDB 教学库的 `system.bit_sequences`、`system.indexes.bitmap_indexes`、`system.indexes.bloom_filters`。

所有 A 类代码都可以直接 `python xxx.py` 运行，无需仓库文件。

## 五、实验步骤

### 步骤 1：造数据，先做一次最朴素的扫描

```python
import time
import numpy as np

rng = np.random.default_rng(11)
n = 500_000
city = rng.integers(0, 5, n)        # 5 个不同的城市编号
street = rng.integers(0, 10, n)     # 10 条不同的街道编号


def naive_scan(rep=5):
    """逐行判断：city=0 AND (street=1 OR street=2)"""
    ts, cnt = [], 0
    for _ in range(rep):
        t0 = time.perf_counter()
        cnt = 0
        for i in range(n):
            if city[i] == 0 and (street[i] == 1 or street[i] == 2):
                cnt += 1
        ts.append((time.perf_counter() - t0) * 1000)
    return float(np.median(ts)), cnt


def one_pass():
    """只跑一次的干净版本，用于核对结果"""
    cnt = 0
    for i in range(n):
        if city[i] == 0 and (street[i] == 1 or street[i] == 2):
            cnt += 1
    return cnt


print("命中行数 =", one_pass())
print("期望行数 =", n * (1 / 5) * (2 / 10))
naive_ms, hits = naive_scan()
print(f"朴素逐行扫描：{naive_ms:.3f} ms（5 次中位数），命中 {hits} 行")
```

输出的"期望行数"是 20,000（= 500,000 × 1/5 × 2/10），实测命中会在它附近浮动，这是随机数据的正常波动。

参考实测：命中 19,845 行，朴素扫描 153.588 ms（5 次中位数）。耗时和你的机器有关，但命中数是确定的——随机种子固定，所以每次应该都是 19,845。

### 步骤 2：构建位图索引，用一句话回答同一个查询

每个取值一条位序列，用 Python 大整数表示（第 i 行对应第 i 个比特）。

```python
import time
import numpy as np

rng = np.random.default_rng(11)
n = 500_000
city = rng.integers(0, 5, n)
street = rng.integers(0, 10, n)

# ---- 建索引：扫一遍数据，把每一行的位置置到对应取值的位序列上 ----
t0 = time.perf_counter()
city_bm = [0] * 5
street_bm = [0] * 10
for i in range(n):
    city_bm[city[i]] |= (1 << i)
    street_bm[street[i]] |= (1 << i)
build_ms = (time.perf_counter() - t0) * 1000


def bitmap_query(rep=5):
    ts, cnt = [], 0
    for _ in range(rep):
        t0 = time.perf_counter()
        mask = city_bm[0] & (street_bm[1] | street_bm[2])
        cnt = mask.bit_count()
        ts.append((time.perf_counter() - t0) * 1000)
    return float(np.median(ts)), cnt


bm_ms, bm_hits = bitmap_query()
print(f"构建索引：{build_ms:.1f} ms")
print(f"位图查询（int.bit_count）：{bm_ms:.4f} ms（5 次中位数），命中 {bm_hits} 行")

# ---- 同一个问题：换用字符串方式的 popcount ----
mask = city_bm[0] & (street_bm[1] | street_bm[2])
ts = []
for _ in range(5):
    t0 = time.perf_counter()
    cnt_bin = bin(mask).count("1")          # 先转成二进制字符串再数字符
    ts.append((time.perf_counter() - t0) * 1000)
print(f"仅统计数量：bin().count() {np.median(ts):.4f} ms，"
      f"int.bit_count() 更快；两者结果一致：{cnt_bin == bm_hits}")

# ---- 核对：位图结果必须和朴素扫描完全一致 ----
naive_cnt = sum(1 for i in range(n)
                if city[i] == 0 and (street[i] in (1, 2)))
assert naive_cnt == bm_hits, (naive_cnt, bm_hits)
print("结果与朴素扫描一致 ✓")
```

参考实测：

| 项目 | 耗时 | 备注 |
| --- | --- | --- |
| 构建索引 | 7,660 ms | 一次性投入，扫 50 万行做 50 万次位操作 |
| 朴素逐行扫描 | 153.588 ms | 每次查询都要全表走一遍 |
| 位图索引查询 | 1.5184 ms | AND + OR + 用 `bin(x).count("1")` 统计置位数 |
| 位图查询（换 `int.bit_count()`） | 比上一行快一到两个量级 | 不构造字符串，直接在整数上做 popcount |
| 加速比 | 约 101× | 153.588 ÷ 1.5184 |
| 命中行数 | 19,845 | 与朴素扫描一致 |

注意最后两行：查询耗时的大头往往不是 AND/OR 本身，而是"怎么数 1"。`bin(mask).count("1")` 要先生成一个 50 万字符的字符串，这就是那 1.5184 ms 的来源；Python 3.8 之后提供的 `int.bit_count()` 直接在二进制表示上做 popcount，快得多。参考实测表格采用的是前者（比较保守的口径），换成后者加速比会明显上升——请把你测到的两个数字都记录下来。

同一份代码在另一台机器（Python 3.12）上得到 7,247.9 / 136.793 / 1.4858 ms，加速比约 92×，命中数同样是 19,845。命中数不随机，是因为 numpy 的种子固定；耗时随机器波动很正常。

### 步骤 3：回本次数与"该不该建索引"

把步骤 1、2 的数字代入 3.2 节的算式：

```python
build_ms = 7660.0      # 换成你自己测出来的值
naive_ms = 153.588
bm_ms = 1.5184

breakeven = build_ms / (naive_ms - bm_ms)
print(f"需要查询 {breakeven:.1f} 次以上，索引才回本")
print(f"单次节省 {(naive_ms - bm_ms):.3f} ms")
print(f"加速比 {naive_ms / bm_ms:.1f}×")
print(f"索引读取压缩：{(500_000 / 8 / 1024):.1f} KiB / 每条位序列")
```

按参考实测，答案是约 50 次。也就是说这张表只要被按同样的模式查询 50 次以上，建位图索引就是划算的；如果这张表每天只查两三次，这笔投资要一年才回本。

### 步骤 4：布尔列表 vs 大整数，同一个结构的两种物理表示

```python
import time
import numpy as np

rng = np.random.default_rng(11)
n = 200_000                       # 列表实现太慢，这里用小一号的 N
col = rng.integers(0, 5, n)

# ---------- 表示 1：布尔列表 ----------
t0 = time.perf_counter()
bm_list = [[False] * n for _ in range(5)]
for i in range(n):
    bm_list[col[i]][i] = True
build_list = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
r_list = [a and b for a, b in zip(bm_list[0], bm_list[1])]
q_list = (time.perf_counter() - t0) * 1000

# ---------- 表示 2：Python 大整数 ----------
t0 = time.perf_counter()
bm_int = [0] * 5
for i in range(n):
    bm_int[col[i]] |= (1 << i)
build_int = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
mask_int = bm_int[0] & bm_int[1]
q_int = (time.perf_counter() - t0) * 1000

print(f"{'':12s}{'构建 (ms)':>12s}{'一次 AND (ms)':>16s}{'存储空间（粗估）':>20s}")
print(f"{'布尔列表':<12s}{build_list:>12.1f}{q_list:>16.4f}{'≈ 每个 bool 1 字节以上':>20s}")
print(f"{'大整数':<14s}{build_int:>10.1f}{q_int:>16.5f}{'≈ N/8 字节':>20s}")
print(f"AND 加速比：{q_list / q_int:.0f}×")
```

两种表示的信息内容完全一样，差别只在物理布局：

- 布尔列表的 AND 要在 Python 层做 N 次循环，每次解释器都要耗费若干十几个字节码；
- 大整数的 AND 是一条指令，处理 N 位只需 N/64 个机器字的操作。

这也解释了为什么 Python 教学库里的 `BitListBitSequence` 和 `IntegerBitSequence` 是同一个抽象的两个实现：上层逻辑不动，换底层实现就能拿到数量级的性能差异。

### 步骤 5：等值编码 vs 范围编码

同一个范围查询，两种编码的工作量完全不同。

```python
import time
import numpy as np

rng = np.random.default_rng(11)
n = 200_000
house_number = rng.integers(1, 1001, n)     # 1000 个不同取值

# ---------- 等值编码：E_i[j] = 1 ⟺ A[j] = v_i ----------
t0 = time.perf_counter()
eq_bm = [0] * 1001
for i in range(n):
    eq_bm[house_number[i]] |= (1 << i)
build_eq = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
acc = 0
for v in range(1, 300):          # house_number < 300，要 OR 掉 299 条
    acc |= eq_bm[v]
eq_cnt = acc.bit_count()
q_eq = (time.perf_counter() - t0) * 1000

# ---------- 范围编码：R_i[j] = 1 ⟺ A[j] ≤ v_i ----------
# 先按取值把行号分组，再从最小取值开始向前缀累加。
# 注意：分组时用 .tolist() 转成 Python int，避免 numpy 整数位移溢出。
t0 = time.perf_counter()
house_number_py = house_number.tolist()
groups = {v: [] for v in range(1, 1001)}
for row, v in enumerate(house_number_py):
    groups[v].append(row)

range_bm = [0] * 1001
acc2 = 0
for v in range(1, 1001):
    for row in groups[v]:
        acc2 |= (1 << row)
    range_bm[v] = acc2
build_range = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
res = range_bm[299]                     # 直接取出 "A ≤ 299"，即 house_number < 300
range_cnt = res.bit_count()
q_range = (time.perf_counter() - t0) * 1000

print(f"等值编码：构建 {build_eq:8.1f} ms，范围查询（OR 299 条）{q_eq:9.3f} ms，命中 {eq_cnt}")
print(f"范围编码：构建 {build_range:8.1f} ms，范围查询（取 1 条）{q_range:9.5f} ms，命中 {range_cnt}")
print("理论命中数（200,000 × 299/1000）=", 200_000 * 299 / 1000)
print("结果一致：", eq_cnt == range_cnt)
```

关键观察：

- **等值编码**做范围查询，需要把区间内所有位序列 OR 起来（本例 299 条），查询时间与区间宽度成正比；等值查询则只要 1 条。
- **范围编码**做范围查询只要 1 条位序列，代价是构建更贵（需要按序累加），且更新一个值可能影响后面所有位序列。
- 范围编码要求属性值可排序。这也是为什么教学库里 `RangeEncodedBitmapIndex` 只允许整数属性。

这两种方案的系统比较见 Chan & Ioannidis 1998 的工作。

### 步骤 6：WAH 压缩 —— 稀疏散布 vs 聚簇

```python
import math
import numpy as np

WORD = 32          # 机器字长（位）
FILL_0, FILL_1, LITERAL = 0, 1, 2


def wah_encode(bits: np.ndarray, word: int = WORD) -> int:
    """返回 WAH 编码后占用的字（word）数。bits 为 0/1 数组。"""
    n_words = math.ceil(len(bits) / word)
    out, i = 0, 0
    while i < n_words:
        w = bits[i * word: (i + 1) * word]
        if len(w) < word:
            w = np.concatenate([w, np.zeros(word - len(w), dtype=np.uint8)])
        if np.all(w == 0):
            kind = FILL_0
        elif np.all(w == 1):
            kind = FILL_1
        else:
            kind = LITERAL
        if kind == LITERAL:
            out += 1            # 1 个头 + 1 个字面字（本实现各计 1）
            i += 1
        else:
            run = 0
            while i < n_words:
                w = bits[i * word: (i + 1) * word]
                if len(w) < word:
                    w = np.concatenate([w, np.zeros(word - len(w), dtype=np.uint8)])
                if (kind == FILL_0 and np.all(w == 0)) or (kind == FILL_1 and np.all(w == 1)):
                    run += 1
                    i += 1
                else:
                    break
            out += 1            # 一个填充字表示 run 个同值字
    return out


n = 100_000
rng = np.random.default_rng(21)
print(f"{'密度':>8s}{'未压缩字数':>12s}{'WAH 随机分布':>14s}{'占比':>9s}"
      f"{'WAH 聚簇分布':>14s}{'占比':>9s}")
for d in (0.001, 0.005, 0.01, 0.05, 0.1, 0.25, 0.5):
    rand_bits = (rng.random(n) < d).astype(np.uint8)
    clust_bits = np.zeros(n, dtype=np.uint8)
    n_set = int(n * d)
    start = rng.integers(0, n - n_set)
    clust_bits[start:start + n_set] = 1

    unc = math.ceil(n / WORD)
    r = wah_encode(rand_bits)
    c = wah_encode(clust_bits)
    print(f"{d * 100:>7.1f}%{unc:>12d}{r:>14d}{r / unc * 100:>8.1f}%"
          f"{c:>14d}{c / unc * 100:>8.2f}%")
```

参考实测（N = 100,000 位，机器字 32 位，未压缩 = 100,000/32 = 3,125 字）：

| 置位密度 | 未压缩字数 | 随机分布 WAH | 占比 | 聚簇分布 WAH | 占比 |
| --- | --- | --- | --- | --- | --- |
| 0.1% | 3,125 | 204 | 6.5% | 5 | 0.16% |
| 0.5% | 3,125 | 878 | 28.1% | 5 | 0.16% |
| 1% | 3,125 | 1,532 | 49.0% | 5 | 0.16% |
| 5% | 3,125 | 3,016 | 96.5% | 5 | 0.16% |
| 10% | 3,125 | 3,123 | 99.9% | 5 | 0.16% |
| 25% | 3,125 | 3,125 | 100.0% | 5 | 0.16% |
| 50% | 3,125 | 3,125 | 100.0% | 5 | 0.16% |

两条结论很明显：

1. **随机分布**下，密度超过约 5% 之后压缩基本失效，因为几乎每个字里都既有 0 又有 1，无法凑出长游程。
2. **聚簇分布**下压缩率与密度无关，恒定在 4~5 个字左右：一个连续块只需要"前导 0 游程 + 字面字 + 尾随 0 游程"几个字就能描述。真实数据（按某个键排序后再建的位图）往往聚簇，这正是位图索引在实践中可行的基础。

![WAH 压缩效果](../figures/fig_wah_compress.png)

图 4：WAH 随密度与分布变化的实测结果。

### 步骤 7：Bloom 过滤器 —— 实测假阳性率

用 md5 做 double hashing（一个哈希值拆出两个搬运，再线性组合出 k 个位置），这是 Kirsch–Mitzenmacher 提出的常见做法：不必真的实现 k 个独立哈希函数。

```python
import math
import hashlib


class BloomFilter:
    """教学用 Bloom 过滤器：double hashing，位数组用 Python 大整数。"""

    def __init__(self, m: int, k: int):
        self.m = m
        self.k = k
        self.bits = 0

    def _indices(self, key):
        d = hashlib.md5(str(key).encode()).digest()
        a = int.from_bytes(d[0:4], "big")
        b = int.from_bytes(d[4:8], "big") | 1
        return [(a + i * b) % self.m for i in range(self.k)]

    def add(self, key):
        for i in self._indices(key):
            self.bits |= (1 << i)

    def may_contain(self, key):
        for i in self._indices(key):
            if not (self.bits >> i) & 1:
                return False
        return True


def theory_fpr(m, n, k):
    return (1.0 - math.exp(-k * n / m)) ** k


def empirical_fpr(m, k, n_insert, n_probe, seed=3):
    bf = BloomFilter(m, k)
    for x in range(n_insert):
        bf.add(f"key-{seed}-{x}")
    # 检查抽样的一部分已插入键，确认没有假阴性
    for x in range(0, n_insert, max(1, n_insert // 200)):
        assert bf.may_contain(f"key-{seed}-{x}"), "出现了假阴性，实现有误"
    fp = 0
    for x in range(n_probe):
        if bf.may_contain(f"absent-{seed}-{x}"):
            fp += 1
    return fp / n_probe


n_insert, n_probe = 5_000, 5_000
print(f"{'m/n':>6s}{'k*':>5s}{'实测 FPR':>12s}{'理论 FPR':>12s}")
for bpe in (4, 6, 8, 10, 14, 18, 24):
    m = int(bpe * n_insert)
    k_opt = max(1, round((m / n_insert) * math.log(2)))
    emp = empirical_fpr(m, k_opt, n_insert, n_probe)
    th = theory_fpr(m, n_insert, k_opt)
    print(f"{bpe:>6d}{k_opt:>5d}{emp:>12.4f}{th:>12.4f}")
```

参考实测（n = 5,000 插入，5,000 次"不存在键"探测，k 取理论最优的取整值）：

| m/n（每元素位数） | k | 实测假阳性率 | 理论值 (1−e^{−kn/m})^k |
| --- | --- | --- | --- |
| 4 | 3 | 0.1476 | 0.1469 |
| 6 | 4 | 0.0566 | 0.0561 |
| 8 | 6 | 0.0224 | 0.0216 |
| 10 | 7 | 0.0098 | 0.0082 |
| 14 | 10 | 0.0016 | 0.0012 |
| 18 | 12 | 0.0000 | 0.0002 |
| 24 | 17 | 0.0000 | 0.0000 |

实测值略高于理论值，是正常的：理论公式假设哈希完全随机且忽略了 k 个位置之间的相关性，double hashing 的位置之间存在线性关系，加上 5,000 次探测本身有采样噪声（m/n=14 时理论期望只有 6 个假阳性，多一个就差 0.0002）。

![Bloom 假阳性率](../figures/fig_bloom_fpr.png)

图 5：每元素位数越多，假阳性率指数级下降。

### 步骤 8：k 不是越大越好

固定一个空间预算，扫描不同的 k，看看假阳性率怎么变。

```python
import math

# 下面沿用步骤 7 里定义的 BloomFilter / theory_fpr / empirical_fpr
# （把步骤 7 的代码保存成 bloom_lab.py 时改为：from bloom_lab import ...）
n_insert, n_probe, bpe = 5_000, 5_000, 12
m = bpe * n_insert
k_star = (m / n_insert) * math.log(2)

print(f"固定 m/n = {bpe}，理论最优 k* = {k_star:.2f}")
print(f"{'k':>4s}{'实测 FPR':>12s}{'理论 FPR':>12s}")
emp_all = []
for k in range(1, 15):
    emp = empirical_fpr(m, k, n_insert, n_probe)
    th = theory_fpr(m, n_insert, k)
    emp_all.append(emp)
    print(f"{k:>4d}{emp:>12.4f}{th:>12.4f}")
print("实测最优 k =", emp_all.index(min(emp_all)) + 1)
```

参考实测：理论最优 k* = 8.32；在 k ∈ [1,14] 的网格里，实测最小值落在 k = 10。

这看起来像是"公式失效了"，其实是采样噪声。k = 8 与 k = 10 的理论假阳性率只相差约 0.0002，换算成 5,000 次探测只差 1 个样本的期望计数（15.7 vs 16.7），而这次实验的标准差约 4 个样本。差异完全被噪声淹没。加大探测次数到 50 万，实测最优点会重新靠近 8。

![k 的选择](../figures/fig_bloom_kopt.png)

图 6：固定 m/n = 12 时扫描 k 的结果。曲线是"先降后升"，最小值附近很平。

### 步骤 9（选做）：调用教学库的实现

课程自带的实现把这些抽象包装成了类，行为与本实验的 A 类代码一致，多了"位序列抽象"、"范围/等值编码的公共基类"、"用布尔列表或大整数两种后端"这些工程化内容。

```python
import sys, pathlib
import time

for _p in (pathlib.Path.cwd(), *pathlib.Path.cwd().parents):
    if (_p / "system").is_dir():
        if str(_p) not in sys.path:
            sys.path.insert(0, str(_p))
        break

import itertools as it
from system.data_classes import Address
from system.bit_sequences import IntegerBitSequence, BitListBitSequence, WAHBitSequence
from system.indexes.bitmap_indexes import (
    SortedEqualityEncodedBitmapIndex,
    RangeEncodedBitmapIndex,
)
from system.indexes.bloom_filters import BloomFilter

cities = ["New York", "Cuppertino", "Berlin", "Singapore"]
streets = ["Unistreet", "Macstreet", "Longstreet", "Msstreet", "BigDataStreet"]
house_numbers = list(range(1, 201))

addresses = []
for idx, comb in enumerate(it.product(cities, streets, house_numbers)):
    addresses.append(Address(idx, comb[0], comb[1], comb[2]))

t0 = time.perf_counter()
idx_city = SortedEqualityEncodedBitmapIndex(bit_sequence_type=IntegerBitSequence)
idx_city.bulkload(((a.city, a.id) for a in addresses))
idx_street = SortedEqualityEncodedBitmapIndex(bit_sequence_type=IntegerBitSequence)
idx_street.bulkload(((a.street, a.id) for a in addresses))
build = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
result = idx_city.get_equal("New York") & (
    idx_street.get_equal("Unistreet") | idx_street.get_equal("BigDataStreet")
)
q = (time.perf_counter() - t0) * 1000
print(f"构建 {build:.1f} ms，查询 {q:.4f} ms，命中 {len(result.as_set())} 行")

# 范围编码：范围查询只需一次访问
idx_house = RangeEncodedBitmapIndex(bit_sequence_type=IntegerBitSequence)
idx_house.bulkload(((a.house_number, a.id) for a in addresses))
print("house_number < 100 命中：", len(idx_house.get_smaller(100).as_set()))

# Bloom 过滤器：注意 contains_key 返回真时要回原始数据复核
bf = BloomFilter(
    data_objects=addresses[:400],
    key_attribute_name="id",
    number_of_available_bits=4000,
    bit_sequence_type=BitListBitSequence,
)
print("BloomFilter id=0 可能存在：", bf.contains_key(0))
print("BloomFilter id=3999 可能存在：", bf.contains_key(3999))
```

注意 `BloomFilter` 的构造函数里，哈希函数个数是按 `⌈ln2 · m/n⌉` 自动算的（见 `system/indexes/bloom_filters.py` 的 `__init__`），不需要手动传 k。

## 六、观测记录表

**表 6-1：位图索引 vs 朴素扫描（建议 N = 500,000，city 5 值，street 10 值）**

| 项目 | 你的实测 | 参考实测 |
| --- | --- | --- |
| 构建索引耗时 (ms) | | 7,660 |
| 朴素逐行扫描耗时 (ms) | | 153.588 |
| 位图查询耗时 (ms) | | 1.5184 |
| 加速比 | | 约 101× |
| 命中行数 | | 19,845 |
| 期望命中行数 N×(1/5)×(2/10) | | 20,000 |
| 相对期望的偏差 | | −0.78% |
| 回本所需查询次数 | | 约 50 |

**表 6-2：位序列表示对比（建议 N = 200,000）**

| 表示方式 | 构建耗时 (ms) | 一次 AND (ms) | 粗估存储（每 N 行） |
| --- | --- | --- | --- |
| 布尔列表 | | | ≥ N 字节 |
| Python 大整数 | | | N/8 字节 |

**表 6-3：WAH 压缩（N = 100,000 位，字长 32 位）**

| 密度 | 未压缩字数 | 随机分布字数 | 占比 | 聚簇分布字数 | 占比 |
| --- | --- | --- | --- | --- | --- |
| 0.1% | 3,125 | | | | |
| 1% | 3,125 | | | | |
| 10% | 3,125 | | | | |
| 50% | 3,125 | | | | |

**表 6-4：Bloom 过滤器（n = 5,000，探测 5,000 次）**

| m/n | k | 实测 FPR | 理论 FPR |
| --- | --- | --- | --- |
| 4 | 3 | | 0.1469 |
| 6 | 4 | | 0.0561 |
| 8 | 6 | | 0.0216 |
| 10 | 7 | | 0.0082 |
| 14 | 10 | | 0.0012 |

## 七、思考题

1. 为什么 Bloom 过滤器不会有假阴性（false negative），只会有假阳性？如果把"删除"也支持进去，会出现什么问题？
2. 参考实测里 m/n = 14、k = 10 时理论值是 0.0012，实测是 0.0016。这个差异合理吗？请用统计涨落的角度算一算，5,000 次探测、p = 0.0012 时，期望的假阳性个数和标准差分别是多少。
3. 固定 m/n = 12 时理论最优 k* = 8.32，但实测最优点落在 k = 10。请判断这是公式错了吗，并给出论证。
4. 表 6-3 里随机分布在密度 0.1% 时压缩到 204 个字。请估算这个数字的来源：一个 32 位的字里至少含一个 1 的概率是多少？由此推算"字面字"的期望个数。
5. 为什么聚簇分布下的 WAH 字数几乎与密度无关（一直是 5 个字左右）？请说明一条连续块被 WAH 编码成哪几个字。
6. 范围编码 schemata 构造出的 `R_i` 之间满足什么单调关系？这个关系对"更新一行数据"意味着什么代价？
7. 本实验中的位图查询快是因为`bin(x).count("1")`（或 `int.bit_count()`）。如果查询要求"返回命中的所有行号"而不是"返回数量"，代价会变成什么？请写出提取行号的代码并口述其复杂度。
8. 位图索引适合低基数（取值个数少）的列。如果给一列 100% 唯一的键（例如主键）建`equality位图`，空间会是多大？请用 N = 500,000 算一算，并说明为什么这种情况要用 B+ 树或哈希索引。

## 八、思考题答案与解析

**1. Bloom 没有假阴性的原因与删除问题**

插入键 x 时会把 h₁(x)…h_k(x) 全部置 1，之后再也不会被清 0。下次查询同一个键时，这些位置必然全是 1，所以 `may_contain(x)` 一定返回真。**没有哪个操作会把已置的位清掉，就不存在假阴性。**

加入删除就会破坏这个性质：把某个位置清 0 时，不知道还有多少个其它键的哈希落在这个位置上，清多了会让别的键产生假阴性。解决路线有两条：

- 计数 Bloom（counting Bloom filter）：每个位置用一个小计数器代替 1 位，删除时减 1。代价是空间放大（每个计数器通常 4 位）。
- 不允许删除：定期重建整个过滤器（很多系统就是这么做的）。

**2. m/n = 14 的差异是噪声吗**

一次"不存在键"判定为假阳性是伯努利试验，概率 p = 0.0012，试验次数 n_probe = 5,000：

```
期望个数  E = n_probe · p = 5000 × 0.0012 = 6.0
标准差    σ = √(n_probe · p · (1−p)) = √(5000 × 0.0012 × 0.9988) = √5.993 = 2.45
实测假阳性率 0.0016  ⇒  实测个数 = 5000 × 0.0016 = 8.0 个
偏离        (8 − 6) / 2.45 = 0.82 σ
```

0.82 个标准差是很常见的涨落。反过来说，p 的测量值本身的标准差是 σ/n_probe = 2.45/5000 = 0.00049，所以实测 0.0016 ± 0.0005 完全覆盖理论值 0.0012。差异合理，不需要修正公式。

要想把这条曲线测准，把 n_probe 提到 500,000（增加 100 倍），相对误差就降到原来的 1/10。

**3. k* = 8.32 而实测最优落在 10**

公式没错。原因是这个区域的函数太平坦，被噪声捂住了。用同一个公式算两个点（m/n = 12，即 kn/m = k/12）：

```
k = 8:  p = (1 − e^{−8/12})^8  = (1 − 0.51342)^8  = 0.48658^8  = 0.003142
k = 10: p = (1 − e^{−10/12})^10 = (1 − 0.43460)^10 = 0.56540^10 = 0.003339
```

两者相差 0.0002。换算成 5,000 次探测的期望个数：

```
k = 8  →  5000 × 0.003142 = 15.7 个
k = 10 →  5000 × 0.003339 = 16.7 个
```

而这一次实验的涨落 σ = √(5000 × 0.0031 × 0.9969) ≈ 3.9 个。差 1 个归类噪声。

正确的验证办法是增大 n_probe，或者直接看理论曲线——理论曲线的最小值确实在 8.32 附近。图 6 里理论曲线的谷底形态很平滑，实测曲线的锯齿完全是噪声，两者的最低点是一致的。

顺带注意一个实用点：k 必须是整数，8.32 附近只能用 8 或 9，两者的理论差只有 0.00002，工程上取哪个都行。

**4. 为什么 0.1% 密度下是 204 个字**

一个字的 32 位全为 0 的概率是 (1 − p)^32，p = 0.001：

```
P(字内含至少一个 1) = 1 − (1 − 0.001)^32 = 1 − 0.999^32
0.999^32 = e^{32·ln 0.999} = e^{−0.032016} = 0.96849
⇒ P = 1 − 0.96849 = 0.03151
```

总共 3,125 个字，期望的"字面字"个数：

```
3,125 × 0.03151 = 98.5 个
```

WAH 遇到字面字要额外存一个头字，两个字面字之间还要穿插 0-fill 字。粗略地数：`98.5 个字面字 + 98.5 个头/间隔字 + 首尾两个边界` ≈ 200 上下，与实测 204 吻合。

密度升到 1% 时同一个算式给出：

```
P(字内含至少一个 1) = 1 − 0.99^32 = 1 − 0.72490 = 0.27510
⇒ 字面字期望 = 3,125 × 0.27510 = 860 个
```

按"每个字面字配一个间隔字"估算会是 1,720，但实测只有 1,532。原因是密度上升后相邻的字面字变多，它们会连成一片，共享一个头字，不需要每个都配间隔。这个估算在稀疏时准确，密集时偏高——正好说明 WAH 的字数由"字面字的游程结构"决定，而不只是字面字的个数。

**5. 聚簇分布为什么恒定在 5 个字**

聚簇位序列长这样（一整段连续的 1，其余全是 0）：

```
0 … 0 | 1 1 1 … 1 1 1 | 0 … 0
 └ 前导全 0 字 ┘ └ 若干全 1 字 + 最多两个混杂字 ┘ └ 尾随全 0 字 ┘
```

编码结果是：

```
1 个 0-fill 字（覆盖前导的全部 0 字）
1~2 个字面字（块的头尾例外）：若块恰好字对齐则没有字面字
1 个 1-fill 字（覆盖块内部的全 1 字）
1 个 0-fill 字（覆盖尾随的 0 字）
```

字数是常数，与块的长度无关——这就是游程编码的本质。实测拿到 5 个字，说明这批数据的块边界没有字对齐，头尾各多出一个字面字。

这一条对实践很重要：把表按某个键排序（cluster）之后建位图，压缩率会有数量级的改善。多维聚簇、Z 序曲线之所以有用，根源也在这里。

**6. 范围编码的单调性与更新代价**

范围编码的定义直接给出包含关系：

```
R_1 ⊆ R_2 ⊆ … ⊆ R_d   （因为 A ≤ v_i 蕴含 A ≤ v_{i+1}）
```

也就是说这 d 条位序列构成一个单调链。好处是任意区间查询都能用两三次运算凑出来：

```
A ≤ v_i        →  R_i
A = v_i        →  R_i & ~R_{i-1}
v_i < A ≤ v_j  →  R_j & ~R_i
```

坏处在更新：改一行 A[j] 从 v_a 变成 v_b（设 v_a < v_b），那么 R_a, R_{a+1}, …, R_{b−1} 这 (b−a) 条位序列的第 j 位都要翻转。更新代价与取值距离成正比。对比等值编码，同样这次更新只碰两条位序列。

这是范围编码的定位：适合只读为主的分析型负载（数据仓库、OLAP），不适合高频点更新（OLTP）。Chan & Ioannidis 1998 的工作系统地比较了这些编码在不同查询/更新混合下的表现。

**7. 返回行号要付多少**

统计数量只需要一次 popcount。要拿到行号就得把每一个置位的位置找出来：

```python
def set_positions(x: int):
    """返回整数 x 中所有为 1 的位下标（从低位向高位）。"""
    pos = []
    while x:
        lsb = x & -x              # 取出最低的一个置位
        pos.append(lsb.bit_length() - 1)
        x ^= lsb                  # 清掉这位
    return pos


# 造一个小例子：16 行里第 2、5、9 行命中
demo = (1 << 2) | (1 << 5) | (1 << 9)
print("命中行的行号：", set_positions(demo))

# 若已经算出了 city_bm / street_bm，则：
# mask = city_bm[0] & (street_bm[1] | street_bm[2])
# rows = set_positions(mask)
```

复杂度是 O(h)，h 为命中行数（每条输出一次循环），而不是 O(N)。注意别写成 `[(mask >> i) & 1 for i in range(n)]`——那样是 O(N)，会把位图查询的优势全丢掉。

真实系统里还有一层：拿到行号后要去堆表回表取完整记录，那才是真正的大头。

**8. 给唯一值列建位图有多大**

一列有 V 个不同取值就 V 条位序列，每条 N 位。位图索引的总位数为 V · N。

V = N = 500,000 时：

```
总位数 = 500,000 × 500,000 = 2.5 × 10^11 bit
       = 2.5 × 10^11 / 8 字节 = 3.125 × 10^10 字节
       ≈ 31.25 GB（按 10^9 换算）/ 29.1 GiB（按 2^30 换算）
```

而原始数据区区几十 MB。空间放大了几千倍。

所以低基数列（性别、城市、状态码）适合位图（V 小 ⇒ V·N 小，且 V 条位序列足够表达）；高基数列（主键、身份证号、时间戳）要用 B+ 树或哈希索引，它们的空间是 O(N) 而不是 O(N·V)。经验阈值视系统而定，通常 V 到几千以上就要重新考虑。

## 九、扩展任务（选做）

1. **把 WAH 的 AND 也在压缩态上实现**：不要在解压后做 AND，试着直接遍历两个 WAH 字序列并合并同类型字。比较"解压 → AND → 再压"这条路径的耗时。
2. **加入"脏数据"实验**：把 `street` 的分布从均匀改成 Zipf（幂律），重做 6-1。观察命中行数是否还接近 20,000，讨论倾斜数据下位图选择率的计算为什么更难。
3. **实现 Counting Bloom**：把位数组换成 4 位计数器的数组（可用 `array` 模块），支持 `remove()`，测空间膨胀的倍数。
4. **测 wordsize 的影响**：把 `WAHBitSequence.WORD_LENGTH` 从 4 一路调到 256，记录每种配置下的存储位数与一次 AND 的耗时，画出"字长 vs 空间"、"字长 vs 时间"两条曲线，讨论为什么字长存在拐点。
5. **对比 Roaring**：安装 `pyroaring`（若网络允许）或自己实现一个简化版（按高 16 位分桶 + 容器内用 sorted array），与 WAH 在同一批数据上比较空间与 AND 耗时。
6. **做一次端到端的流水线**：把位图查询结果还原成行号，用这些行号去"回表"取出完整记录，再用 Bloom 过滤器挡住不可能命中的 id。统计整条链路各阶段耗时，画出占比图，看看时间到底花在哪一步。

## 十、常见错误与排查

| 现象 | 原因 | 处理 |
| --- | --- | --- |
| 位图查询与朴素扫描结果不一致 | 位对应关系错位：习惯性的第 i 行写成 `1 << i` 却在某个循环里用了 `1 << (i+1)` | 保持"行号 = 比特下标"的唯一约定；在两种实现里都用一个小型的对拍脚本验证 |
| `1 << i` 这一步跑得特别慢 | 用位循环逐位构建 N = 500,000 的大整数会反复重新分配大对象 | 这是本实验里最贵的部分（构建 7,660 ms 主要花在这）；减少 N 或用 `numpy.packbits` + `int.from_bytes` 一次性构造 |
| `int.bit_count()` 报 `AttributeError` | Python 版本低于 3.8（3.10 起 `int.bit_count` 稳定可用） | 换成 `bin(x).count("1")`；或升级到 Python 3.10+ |
| `np.floor` 类型错误：`ufunc 'bitwise_and' not supported` | 把 numpy 的 bool 数组和 Python 大整数混在一起做 `&` | 二者选一：要么全程纯 Python 大整数（本实验推荐），要么全程 numpy + `np.packbits` |
| 加了 `assert` 却依然不生效 | numpy 标量的比较结果不是 Python 的 bool，而是 numpy 布尔类型 | 断言前先转换：`bool(cond)` 或用 `int(...)` 比较 |
| Bloom 实测出现了"假阴性"（已插入的键查不到） | 哈希函数实现里 k 个位置重复，或 `add` 与 `may_contain` 用了不同的哈希逻辑 | `_indices` 必须只定义一次，`add` 与 `may_contain` 都调用它；确认 b 为奇数，防止退化成同一个位置 |
| Bloom 实测假阳性率明显高于理论值（比如 2 倍以上） | m 与 k 不匹配：m/n = 10 却用了 k = 2，或 n 的计数算错（把成功了插入的重算为不同值） | 用 `k_opt = round((m/n) · ln2)` 核对；检查 `n_insert` 实际插入了多少个不同的键 |
| Bloom 的 hashed key 写法变了以后，实测数字对不上表格 | 键串生成器（"key-seed-x"）改变会使哈希分布改变 | 保持与本文一致的 `f"key-{seed}-{x}"` 写法；这属于实验配置，修改后数字当然会变 |
| WAH 压缩率总是 100% | 传入的位序数组密度太高（比如随机分布、密度 50%），或 word 遍历时 `np.all(w == 0)` 判断块错位 | 确认 `bits` 是 0/1 的 `uint8` 数组；用人工构造的 `[0]*100 + [1]*5 + [0]*100` 样例先验证 `wah_encode` 是否返回预期的小数字 |
| 调用 `system.indexes.bitmap_indexes` 时 `bulkload` 卡住 | 数据量设得太大：house_number 到 1000 会生成约 4000 行，范围编码要建 1000 条位序列 | 先用 `range(1, 201)` 的小规模跑通，再调大 |
| 使用 `system` 包时 `ModuleNotFoundError` | 未把仓库根加入 `sys.path` | 按 4.3 节/步骤 9 的代码片段先插入路径 |
| `RangeEncodedBitmapIndex` 传入字符串属性报类型错误 | 范围编码要求可排序，教学库限定整数属性 | 给 house_number 这类有序属性建范围编码；字符串属性用 `SortedEqualityEncodedBitmapIndex` |