{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 实验 16：Z 序曲线（Z-order / Morton code）\n",
    "\n",
    "## 为什么做这个实验\n",
    "\n",
    "数据库里最好用、最成熟的索引是 B+ 树，可它只能对**一个**排序键做范围查找。\n",
    "\"找出某个矩形框内的所有点\"是个二维问题，一个 B+ 树区间表达不了。\n",
    "\n",
    "绕过去的办法是找一个映射\n",
    "\n",
    "    f: (x, y) → z（一个整数）\n",
    "\n",
    "让**空间上邻近的点，在一维顺序里也尽量邻近**。这样二维的矩形查询就能近似变成\"若干个 z 值区间\"，\n",
    "直接复用现成的 B+ 树。这类映射叫**空间填充曲线**（space-filling curve，SFC）；\n",
    "Z 序（Z-order，也叫 Morton 码）是其中构造最简单的一种——只需要把两个坐标的二进制位**交织**起来。\n",
    "\n",
    "本实验要弄清楚四件事：\n",
    "\n",
    "1. 位交织怎么算、怎么反算（它是严格的一一映射，可逆）；\n",
    "2. 为什么\"坐标高位相同 ⇒ 排序后仍邻近\"，也就是 Z 序局部性的来源；\n",
    "3. 用普通一维 B+ 树索引二维数据的**两阶段流程**：先取候选、再精确过滤；\n",
    "4. 在仓库自带的真实经纬度数据上，把 Z 序与顺序扫描、与 Hilbert 曲线的**候选率**量出来，\n",
    "   并说清楚这个对比的方法学前提。\n",
    "\n",
    "一句话预告结论：Z 序不是\"最好\"的曲线，但它是**最便宜**的曲线——位交织一次就算完，\n",
    "工业界大量系统（Geohash、S2 的部分环节）仍在用它。\n",
    "\n",
    "## 运行前准备\n",
    "\n",
    "- Python 3.9 及以上。\n",
    "- 必选依赖：`numpy`（整数版 Morton 码与真实数据读取）、`matplotlib`（画图）。\n",
    "- 可选依赖：`ipycanvas`（本 notebook 最后一个 cell 用它画交互式画布）。\n",
    "  装不上也没关系，改用下面给出的 matplotlib 版本即可。\n",
    "- 缺包时安装：\n",
    "\n",
    "  ```bash\n",
    "  python -m pip install numpy matplotlib\n",
    "  ```\n",
    "\n",
    "  （本机已配置环境可用：`python -m pip install numpy matplotlib`）\n",
    "\n",
    "- 数据文件：`notebooks/latitude_longitude.bin`（仓库自带，197 229 个经纬度点）。\n",
    "  读取数据时要 `os.chdir` 到仓库根目录，或用绝对路径。\n",
    "- 打开方式：在仓库根目录执行 `jupyter notebook notebooks/Z-Order-Curve.ipynb`。\n",
    "\n",
    "**关于字符串实现**：下面的教学版用字符串表示 z 码，一眼能看出每一位从哪来。\n",
    "真实系统用整数（位运算），并且要一次算几百万个点——末尾给出了 numpy 向量化版本。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 空间填充曲线家族与位交织规则\n",
    "\n",
    "常见的空间填充曲线：\n",
    "\n",
    "- **Z 序**（Z-order / Morton order）：本实验的主角，构造最简单，只需位交织；\n",
    "- **Gray 码变体**：把 Z 序的码字再套一次格雷码变换，缓解边界跳跃；\n",
    "- **Hilbert 曲线**：递归的\"U\"形拼接，通常被认为局部性更好（**但有争议，见末尾**）；\n",
    "- **Peano 曲线 / Sierpinski 曲线**：其他递归构造。\n",
    "\n",
    "### 位序约定（本 notebook 采用）\n",
    "\n",
    "设递归层级为 L，x、y 各用 L 个二进制位表示，z 码共 2L 位。\n",
    "记 x 的各位为 a_{L−1} … a₁ a₀、y 的各位为 b_{L−1} … b₁ b₀（下标 0 是最低位），则\n",
    "\n",
    "    z = ( b₂ a₂  b₁ a₁  b₀ a₀ )₂        （以 L = 3 为例）\n",
    "\n",
    "即：**y 的一位、x 的一位**交替，从高位排到低位；y 的位落在 z 的奇数位，x 的位落在偶数位。\n",
    "\n",
    "举例（L = 3，x = 2，y = 1）：\n",
    "\n",
    "    x = 2 = 010        y = 1 = 001\n",
    "                    y₂ x₂  y₁ x₁  y₀ x₀\n",
    "    z =               0  0   0  1   1  0   = \"000110\"\n",
    "\n",
    "### 网格规模\n",
    "\n",
    "| 递归层级 L | 网格 | 单元数 | z 码位宽 |\n",
    "| --- | --- | --- | --- |\n",
    "| 1 | 2 × 2 | 4 | 2 |\n",
    "| 2 | 4 × 4 | 16 | 4 |\n",
    "| 3 | 8 × 8 | 64 | 6 |\n",
    "| L | 2^L × 2^L | 2^(2L) | 2L |\n",
    "\n",
    "### 量化：浮点坐标必须先变成整数\n",
    "\n",
    "真实经纬度是浮点数，位交织只认整数，所以要先**量化**到整数网格上。本实验用线性映射：\n",
    "\n",
    "    x = round( (lon + 180) / 360 × (2^L − 1) )\n",
    "    y = round( (lat +  90) / 180 × (2^L − 1) )\n",
    "\n",
    "量化层级越高越精确，但码越长；量化本身也是一次精度损失——这是所有空间索引共同的取舍。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第 1 步：位交织（坐标 → z 码）\n",
    "\n",
    "`scalar_to_z_code(recursion_level, x, y)` 做三件事：\n",
    "\n",
    "1. 把 x、y 转成二进制字符串，并用 `zfill(recursion_level)` 左侧补零到 L 位；\n",
    "2. 断言两者长度都等于 L（防止坐标超出网格范围）；\n",
    "3. 用 `\"\".join(map(add, y_binary, x_binary))` 把 y 与 x **逐位拼接**——\n",
    "   `map(add, y, x)` 把两个字符串的对应字符两两相加，于是第 i 次得到 y[i] + x[i]，\n",
    "   拼起来正好是 `y_{L−1} x_{L−1} … y₀ x₀`。这就是位交织的全部实现。\n",
    "\n",
    "预期看到：`scalar_to_z_code(3, 2, 1)` 输出 `'000110'`，与上一节的手算一致。\n",
    "\n",
    "注意这是**字符串**实现，慢且占内存，只用于教学。真实系统用整数位运算（末尾给出）。"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "from operator import add\n",
    "\n",
    "\n",
    "def scalar_to_z_code(recursion_level: int, x: int, y: int) -> str:\n",
    "    \"\"\"\n",
    "    Converts scalar coordinates into a z-code.\n",
    "\n",
    "    Notice that we use a string-representation of z-codes for educational purposes\n",
    "    a real implementation would represent z-codes in a bit-vector or map them to signed integers\n",
    "    :param recursion_level: the recursion_level of the z-code\n",
    "    :param x: x-coordinate\n",
    "    :param y: y-coordinate\n",
    "    :return: z-code as a string\n",
    "    \"\"\"\n",
    "    if recursion_level == 0:\n",
    "        # return empty string for recursion_level 0\n",
    "        return \"\"\n",
    "\n",
    "    # get the binary representation of x and y as a string and prefix them with zeros:\n",
    "    x_binary = bin(x)[2:].zfill(recursion_level)\n",
    "    y_binary = bin(y)[2:].zfill(recursion_level)\n",
    "\n",
    "    assert len(x_binary) == len(y_binary) == recursion_level\n",
    "\n",
    "    # zip the binary representations of x and y and interleave them#\n",
    "    # this is the actual \"trick\" of z-codes: the bits of the n=2 binary numbers are interleaved\n",
    "    ret = \"\".join(map(add, y_binary, x_binary))\n",
    "\n",
    "    assert len(ret) == 2 * recursion_level\n",
    "\n",
    "    return ret\n",
    "\n",
    "\n",
    "scalar_to_z_code(3, 2, 1)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第 2 步：逆变换（z 码 → 坐标）\n",
    "\n",
    "`z_code_to_scalar(z_code)` 把交织好的码拆回去：\n",
    "\n",
    "- 码长除以 2 就是递归层级 L；\n",
    "- 从**最高位**开始，每 2 位一组：`z_code[i]` 属于 y，`z_code[i+1]` 属于 x；\n",
    "- 每组的权重是 2^exponent，exponent 从 L−1 递减到 0。\n",
    "\n",
    "所以它能还原，是因为位交织是**严格的一一映射**：z 码里从低位起的第 0、2、4… 位就是 x 的各位，\n",
    "第 1、3、5… 位就是 y 的各位，按原权重累加回去即可，信息一点没丢。\n",
    "\n",
    "最后一个 cell 是**往返验证**：`z_code_to_scalar(scalar_to_z_code(2, 3, 1))` 应返回原坐标。\n",
    "如果这里对不上，99% 是位序写反了（把 x 放到了高位）——这是本实验最常见的错误。"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def z_code_to_scalar(z_code: str) -> tuple[int, int, int]:\n",
    "    \"\"\"\n",
    "    inverse of scalar_to_z_code function\n",
    "\n",
    "    converts z-code to scalar-coordinates:\n",
    "    :param z_code: string-representation of a z-code\n",
    "    :return: recursion_level, scalar x, scalar y\n",
    "    \"\"\"\n",
    "    assert len(z_code) % 2 == 0\n",
    "\n",
    "    dimensions: int = 2\n",
    "\n",
    "    recursion_level: int = len(z_code) // dimensions\n",
    "    x: int = 0\n",
    "    y: int = 0\n",
    "    exponent: int = int(recursion_level) - 1\n",
    "\n",
    "    i: int\n",
    "    for i in range(0, len(z_code), 2):\n",
    "        # reconstruct x and y by reversing the interleaving of the z-code\n",
    "        # basically builds sum of 2^i exponents\n",
    "        y += int(z_code[i]) * 2**exponent\n",
    "        x += int(z_code[i + 1]) * 2**exponent\n",
    "        exponent -= 1\n",
    "\n",
    "    return recursion_level, x, y\n",
    "\n",
    "\n",
    "z_code_to_scalar(\"000110\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "z_code_to_scalar(scalar_to_z_code(2, 3, 1))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第 3 步：把整张表的 z 码打印出来\n",
    "\n",
    "辅助函数 `color_code` 只是给码字分段上 ANSI 颜色（终端转义序列），\n",
    "让相邻两位的来源在屏幕上更容易分辨；它不影响码值本身。\n",
    "\n",
    "- `print_z_codes(L)`：按 y 从 0 到 2^L−1、x 从 0 到 2^L−1 的顺序打印整个网格的 z 码；\n",
    "- `compute_z_codes(L)`：返回 `(z_code, x, y)` 三元组列表，供后面画图用。\n",
    "\n",
    "观察要点（用 L = 3 试）：\n",
    "\n",
    "- 沿同一行看，变化的是 z 码的**低位**；\n",
    "- 跨行时变化的是**高位**；\n",
    "- 整体读下来，码值从小到大的路径正是那个\"Z\"字形。\n",
    "\n",
    "这也直观说明了局部性的来源：**码的高位来自坐标的高位**。\n",
    "两个点的高位相同 ⇒ 它们落在同一个粗粒度的象限里 ⇒ 它们的 z 码也共享高位\n",
    "⇒ 按数值排序后必然连成一段连续区间，在一维顺序里仍然邻近。\n",
    "\n",
    "反过来，一旦跨象限，高位一变，一维顺序上就会出现一次**大跳跃**——\n",
    "Z 曲线在象限边界上的那些长跳线，就是这个效应的可视化，也是 Z 序的主要缺点。"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "startEsc: str = \"\\033[1;\"\n",
    "middleEsc: str = \"m\"\n",
    "endEsc: str = \"\\033[1;m\"\n",
    "\n",
    "\n",
    "def color_code(string: str, stepSize: int = 2):\n",
    "    \"\"\"color-codes a string in steps of <stepSize> with different colors\n",
    "\n",
    "    :param string: the string to color-code\n",
    "    :param stepSize: the step-size for coloring\n",
    "    \"\"\"\n",
    "    color: int = 31\n",
    "    i: int\n",
    "    for i in range(0, len(string), stepSize):\n",
    "        prefix: str = string[i : i + stepSize]\n",
    "        print(startEsc + str(color) + middleEsc + prefix + endEsc, end=\"\")\n",
    "        color += 1\n",
    "\n",
    "\n",
    "color_code(\"hello world\")"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def print_z_codes(recursion_level: int):\n",
    "    \"\"\"print all z-codes of a given recursion_level\n",
    "\n",
    "    :param recursion_level: the recursion level\n",
    "    \"\"\"\n",
    "    y: int\n",
    "    for y in range(2**recursion_level):\n",
    "        x: int\n",
    "        for x in range(2**recursion_level):\n",
    "            z_code: str = scalar_to_z_code(recursion_level, x, y)\n",
    "            color_code(z_code)\n",
    "            print(\" \", end=\"\")\n",
    "        print()\n",
    "\n",
    "\n",
    "print_z_codes(3)"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "def compute_z_codes(recursion_level: int):\n",
    "    \"\"\"computes all z-codes for the given recursion_level\n",
    "\n",
    "    @param recursion_level: the recursion level\n",
    "    @return: a list of tuples containing the z-code, scalar x, and scalar y coordinates\n",
    "    \"\"\"\n",
    "    results: list[tuple[str, int, int]] = []\n",
    "    y: int\n",
    "    for y in range(2**recursion_level):\n",
    "        x: int\n",
    "        for x in range(2**recursion_level):\n",
    "            z_code: str = scalar_to_z_code(recursion_level, x, y)\n",
    "            results.append((z_code, x, y))\n",
    "\n",
    "    return results"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第 4 步：把 Z 曲线画出来\n",
    "\n",
    "`visualize_z_codes(L)` 用 `ipycanvas` 画交互式画布：它先算出整个网格的 z 码，\n",
    "**按码值排序**，再依次连线，于是笔尖走出的轨迹就是 Z 曲线。\n",
    "\n",
    "预期看到（L 从 1 到 5 依次画）：每一级递归都在四个象限里重复一个\"Z\"字，\n",
    "层级越高，\"Z\"字越小越密，整体是一条自相似的折线。\n",
    "\n",
    "**如果 `ipycanvas` 装不上**，用下面这个 matplotlib 版本替代（只依赖 numpy + matplotlib）：\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def morton(x, y, levels=16):\n",
    "    \"\"\"向量化的 Morton 码：y 位放奇数位，x 位放偶数位，与教学版位序一致。\"\"\"\n",
    "    x = np.asarray(x, dtype=np.int64)\n",
    "    y = np.asarray(y, dtype=np.int64)\n",
    "    z = np.zeros_like(x)\n",
    "    for i in range(levels):\n",
    "        z |= ((y >> i) & 1) << (2 * i + 1)\n",
    "        z |= ((x >> i) & 1) << (2 * i)\n",
    "    return z\n",
    "\n",
    "\n",
    "def draw_z_curve(L=3):\n",
    "    n = 1 << L\n",
    "    xs, ys = np.meshgrid(np.arange(n), np.arange(n))\n",
    "    xs, ys = xs.ravel(), ys.ravel()\n",
    "    order = np.argsort(morton(xs.astype(np.int64), ys.astype(np.int64), levels=L),\n",
    "                       kind=\"stable\")\n",
    "    plt.rcParams[\"font.sans-serif\"] = [\"Microsoft YaHei\", \"SimHei\", \"DejaVu Sans\"]\n",
    "    fig, ax = plt.subplots(figsize=(6, 6))\n",
    "    ax.plot(xs[order], ys[order], \"-o\", color=\"#1F4E79\", ms=8, lw=2)\n",
    "    ax.set_aspect(\"equal\")\n",
    "    ax.set_xlabel(\"x\")\n",
    "    ax.set_ylabel(\"y\")\n",
    "    ax.set_title(f\"Z 曲线（递归层级 {L}，{n}×{n} 网格）\")\n",
    "    ax.grid(ls=\":\", alpha=0.4)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "draw_z_curve(3)\n",
    "```\n",
    "\n",
    "把 L 改成 1、2、4 各画一次，看递归结构怎么长出来。"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "source": [
    "from ipycanvas import Canvas, hold_canvas\n",
    "from IPython.display import display\n",
    "\n",
    "\n",
    "def visualize_z_codes(recursion_level: int):\n",
    "    \"\"\"visualize the data space in z-code order, i.e. draw the z-curve\n",
    "\n",
    "    :param recursion_level: the recursion level\n",
    "    \"\"\"\n",
    "    results: list[tuple[str, int, int]] = compute_z_codes(recursion_level)\n",
    "    results.sort()\n",
    "    x_width: int = 300\n",
    "    y_width: int = x_width\n",
    "    canvas = Canvas(size=(x_width, y_width), sync_image_data=True)\n",
    "    canvas.line_width = 1\n",
    "    max_x: int = 2**recursion_level\n",
    "    max_y: int = max_x\n",
    "    mult_x: float = x_width / max_x\n",
    "    mult_y: float = y_width / max_y\n",
    "    start: tuple[str, int, int] = results.pop(0)\n",
    "    start_x: int = start[1]\n",
    "    start_y: int = start[2]\n",
    "\n",
    "    with hold_canvas(canvas):\n",
    "        while len(results) > 0:\n",
    "            next: tuple[str, int, int] = results.pop(0)\n",
    "            next_x: int = next[1]\n",
    "            next_y: int = next[2]\n",
    "            canvas.begin_path()\n",
    "            canvas.move_to((start_x + 0.5) * mult_x, (start_y + 0.5) * mult_y)\n",
    "            canvas.line_to((next_x + 0.5) * mult_x, (next_y + 0.5) * mult_y)\n",
    "            canvas.stroke()\n",
    "            start_x, start_y = next_x, next_y\n",
    "\n",
    "    return canvas\n",
    "\n",
    "\n",
    "# visualize the z-codes for recursion levels 1 to 5\n",
    "for _recursion_level in range(1, 6):\n",
    "    display(visualize_z_codes(_recursion_level))"
   ],
   "outputs": [],
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 在真实数据上量一量：197 229 个经纬度点\n",
    "\n",
    "仓库自带 `notebooks/latitude_longitude.bin`：\n",
    "\n",
    "- 197 229 个点；\n",
    "- **float32 小端**存储，纬度、经度**交替**存放：`lat = raw[0::2]`、`lon = raw[1::2]`；\n",
    "- 纬度范围 47.28° – 55.07°，经度范围 −2.01° – 50.51°（欧洲中北部一带）。\n",
    "\n",
    "### 用 z 码做范围查询：两阶段流程\n",
    "\n",
    "    ① 候选阶段（filter）：把查询矩形翻译成一段（或几段）z 区间，用一维索引（B+ 树 / 有序数组 + 二分）取出区间内所有点\n",
    "    ② 精化阶段（refine）：对候选点逐个用真实坐标判断是否在矩形内\n",
    "\n",
    "第 ① 步取到的是矩形的**超集**：一维区间无法精确表达二维矩形，\n",
    "所以第 ② 步**不可省略**。工业实现会用 `BIGMIN` / `LITMAX` 把大区间切成若干真正落在矩形内的小段，\n",
    "跳过大段空洞，从而少取候选。\n",
    "\n",
    "衡量指标：**候选率 = 候选点数 / 真实命中点数**。越接近 1 越好；\n",
    "顺序扫描的候选率 = 总点数 / 命中点数。\n",
    "\n",
    "### 实测结果（同一 2^16 网格，89 个有效查询窗口，取中位数）\n",
    "\n",
    "| 方案 | 候选数 / 真实命中数 |\n",
    "| --- | --- |\n",
    "| 顺序扫描 | 12.64× |\n",
    "| Z 序索引 | 3.19× |\n",
    "| Hilbert 序索引 | **1.77×** |\n",
    "\n",
    "（口径：随机生成小窗口矩形查询；**只统计真实命中数 ≥ 5 的查询**——命中太少的查询不足以体现差异，\n",
    "本次有效窗口 89 个；Z 序与 Hilbert 序建立在**同一个** 2^16 网格上。数字以本案例库\n",
    "《00-实测数字口径表》为准，本机重跑会随随机种子略有差异。）\n",
    "\n",
    "复现要点：量化到 2^16 网格 → 算码 → 按码排序 → 用 `np.searchsorted` 取区间\n",
    "（右端用 `side=\"right\"`）→ 与真实命中数求比值 → 取中位数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 一个必须讲的方法学教训：比较两条曲线要站在同一把尺子上\n",
    "\n",
    "本案例库早期版本做过一次对比：Z 码用 2^16 网格、Hilbert 码用 2^10 网格。\n",
    "结果 Hilbert 看起来领先很多——但那个\"优势\"是假的。\n",
    "\n",
    "原因：网格越粗，同一个查询框在码空间里覆盖的区间越短，候选集天然就更少。\n",
    "也就是说，测出来的差距有相当一部分来自**量化层级不同**，而不是曲线本身的性质。\n",
    "那个数字已经作废，**不要再引用**。\n",
    "\n",
    "现在把两者统一到 2^16 网格后，Hilbert 相对 Z 序的领先幅度从约 2.5 倍回落到约 1.8 倍（3.19 → 1.77）。\n",
    "这才是可复现、可辩护的对比。\n",
    "\n",
    "结论值得单独记住：**任何对比实验，只要两方的\"其他条件\"不一致，读出来的差距就不可信。**\n",
    "比较两条空间填充曲线之前，先确认它们站在同一把尺子上——同网格、同查询集、同过滤阈值、同统计口径。\n",
    "\n",
    "### 顺便说一句：别把\"Hilbert 一定比 Z 序好\"当定论\n",
    "\n",
    "这个说法**有争议**，取决于用什么度量。Dai & Su（DMTCS 2003）在\"平均跨簇距离\"与\n",
    "\"平均总跨簇距离\"两项统计上得到的结论是 **Z 序家族优于 Hilbert 家族**，与很多教材的说法相反。\n",
    "\n",
    "本实验在\"小窗口范围查询的候选率\"这个口径下 Hilbert 更好（1.77× vs 3.19×），\n",
    "但**换一个度量就可能翻转**。工程上的正确做法是：用自己的真实数据和真实查询分布跑一遍，\n",
    "按候选率（或干脆按端到端查询耗时）来选，而不是照搬结论。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Z 序在真实系统里的落地形态\n",
    "\n",
    "| 系统 / 结构 | 做法 | 备注 |\n",
    "| --- | --- | --- |\n",
    "| **Geohash** | 基于 Z 序，把经纬度交替二分后编成 base32 字符串 | 前缀相同表示区域相邻，常用作前缀匹配的粗过滤 |\n",
    "| **Google S2** | 用**六条 Hilbert 曲线**首尾相连覆盖球面（立方体投影的六个面） | 层级 level 0–30，单元 ID 是 64 位整数 |\n",
    "| **Uber H3** | 六边形层次网格，分辨率 0–15 | 总单元数 c(r) = 2 + 120·7^r；**每层恰有 12 个五边形**（六边形无法完美铺满球面，这是拓扑上不可避免的） |\n",
    "| **R 树** | 用最小外接矩形（MBR）组织空间对象，不是线性化路线 | 出处：Guttman，《R-trees: A Dynamic Index Structure for Spatial Searching》，SIGMOD 1984 |\n",
    "| **PostGIS** | GiST / SP-GiST 索引 | 见下方 `ST_DWithin` 的提醒 |\n",
    "| **GeoParquet** | 列式空间数据文件格式，带 bbox 统计便于裁剪 | 配合 DuckDB / Sedona 做列存空间分析 |\n",
    "| **Apache Sedona** | 在 Spark 之上提供空间分区、空间索引与空间连接算子 | 海量空间数据的分布式处理路径 |\n",
    "\n",
    "三者取舍：Geohash / S2 的四边形单元便于位运算与前缀查询；H3 的六边形单元到邻居的距离更均匀\n",
    "（每个单元有 6 个等距邻居），适合邻域聚合，代价是那 12 个五边形。\n",
    "\n",
    "**PostGIS 的一个高频坑**：`ST_DWithin` 是**索引感知**（index-aware）的距离阈值谓词，\n",
    "能用 GiST 索引先做包围盒粗过滤；而 `ST_Distance(...) < d` 只是一个普通表达式，\n",
    "优化器不知道它等价于\"距离阈值\"，只能逐行算距离再比较，**不走索引**，大表上退化成全表扫描。\n",
    "（依据：PostGIS 官方文档 Chapter 5。）\n",
    "\n",
    "```sql\n",
    "-- 推荐：能走 GiST 索引\n",
    "SELECT * FROM restaurants\n",
    "WHERE ST_DWithin(geom, ST_MakePoint(:lon, :lat)::geography, 500);\n",
    "\n",
    "-- 不推荐：全表扫描\n",
    "SELECT * FROM restaurants\n",
    "WHERE ST_Distance(geom, ST_MakePoint(:lon, :lat)::geography) < 500;\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 思考题与详解\n",
    "\n",
    "**1. 手算 `scalar_to_z_code(3, 5, 2)`：写出 x、y 的二进制与交织过程，再用代码验证。**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "    x = 5 = 101        y = 2 = 010\n",
    "    位序   y₂ x₂  y₁ x₁  y₀ x₀\n",
    "           0  1   1  0   0  1\n",
    "    z 码 = \"011001\"\n",
    "\n",
    "即 a₂a₁a₀ = 101（x = 5），b₂b₁b₀ = 010（y = 2），按 y 位在前、x 位在后交织得 `0 1 1 0 0 1`。\n",
    "\n",
    "验证：`print(scalar_to_z_code(3, 5, 2))` 输出 `'011001'`。\n",
    "\n",
    "</details>\n",
    "\n",
    "**2. 求 `z_code_to_scalar(\"011001\")` 的结果，并说明逆变换为什么能还原出 x 与 y。**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "码长 6 ⇒ L = 3。按 2 位一组拆：`01 | 10 | 01`，权重从 2² 递减到 2⁰：\n",
    "\n",
    "    第 1 组 \"01\"（权重 4）：y += 0×4 = 0，x += 1×4 = 4\n",
    "    第 2 组 \"10\"（权重 2）：y += 1×2 = 2，x += 0×2 = 0\n",
    "    第 3 组 \"01\"（权重 1）：y += 0×1 = 0，x += 1×1 = 1\n",
    "    结果：x = 5，y = 2  →  返回 (3, 5, 2)\n",
    "\n",
    "能还原的原因是位交织是**一一映射**：从 z 码最低位起，第 0、2、4… 位恰好是 x 的各位，\n",
    "第 1、3、5… 位是 y 的各位，把它们按原来的权重 2^i 累加回去即可，信息没有丢失。\n",
    "\n",
    "（若你的结果对不上，先检查位序：本 notebook 约定 **y 在奇数位、x 在偶数位**。）\n",
    "\n",
    "</details>\n",
    "\n",
    "**3. 递归层级 L = 16 时，网格边长是多少？单元总数是多少？z 码需要几个二进制位、能否放进一个 64 位整数？**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "- 网格边长 = 2^16 = **65 536**\n",
    "- 单元总数 = 2^16 × 2^16 = 2^32 = **4 294 967 296**（约 43 亿）\n",
    "- z 码位宽 = 2L = **32 位**，**能**放进 64 位整数（也能放进无符号 32 位整数）\n",
    "\n",
    "延伸：若 L 取到 32，z 码就是 64 位，此时必须注意有符号整数溢出，要用 `uint64`；\n",
    "Python 的 `int` 虽然不会溢出，但与 C/NumPy 交互时类型要显式指定。\n",
    "\n",
    "</details>\n",
    "\n",
    "**4. 为什么范围查询必须\"先取候选、再精确过滤\"？举一个具体反例：某点的 z 值落在 `[z(x₀,y₀), z(x₁,y₁)]` 内，但点本身不在查询矩形里。**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "取 L = 2（4×4 网格），查询矩形 x ∈ [1, 2]、y ∈ [0, 1]。\n",
    "\n",
    "    z(1, 0)：x = 01, y = 00  →  y₁ x₁ y₀ x₀ = 0 0 0 1  = 1\n",
    "    z(2, 1)：x = 10, y = 01  →  y₁ x₁ y₀ x₀ = 0 1 1 0  = 6\n",
    "    候选区间 = [1, 6]\n",
    "\n",
    "区间里的 z = 2 对应二进制 `0010`：`y₁ x₁ y₀ x₀ = 0 0 1 0` ⇒ y = `01` = 1，x = `00` = 0，\n",
    "即点 **(0, 1)**。它的 x = 0 不在 [1, 2] 内——**落在 z 区间里却不在矩形里**。\n",
    "\n",
    "原因：**一维区间无法精确表达二维矩形**，z 区间只是矩形的一个超集。所以精化阶段\n",
    "（逐个用真实坐标判断）不可省略。工业实现用 `BIGMIN` 算法把大区间切成若干\"真正落在矩形内\"的小段，\n",
    "跳过大段空洞，从而减少这类无用候选。\n",
    "\n",
    "</details>\n",
    "\n",
    "**5. 关于 Hilbert 与 Z 序的优劣，文献结论一致吗？Dai & Su（DMTCS 2003）得到的是什么结论？工程上应该怎么处理？**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "**不一致。** Dai & Su（DMTCS 2003）在\"平均跨簇距离\"与\"平均总跨簇距离\"两项统计上，\n",
    "得到的结论是 **Z 序家族优于 Hilbert 家族**；而很多教材与工程经验又说 Hilbert 局部性更好。\n",
    "\n",
    "分歧来自**度量口径**不同：用哪种距离、按哪种查询分布来统计，结论会变。\n",
    "本实验在\"小窗口范围查询的候选率\"这个口径下 Hilbert 更好（1.77× vs 3.19×），\n",
    "但**换一个度量就可能翻转**。\n",
    "\n",
    "工程做法：**不要照搬结论，用自己的真实数据 + 真实查询分布实测一遍**，\n",
    "按候选率（或干脆按端到端查询耗时）来选曲线与层级。\n",
    "\n",
    "</details>\n",
    "\n",
    "**6. 有人把 Z 码用 2^16 网格、Hilbert 码用 2^10 网格做对比，得出 Hilbert 大幅领先的结论。这个比较公平吗？请设计一个公平的方案。**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "**不公平。** 网格越粗，同一个查询框在码空间里覆盖的区间越短，候选集天然更少。\n",
    "所以那个\"领先\"里混进了量化层级的差异，不是曲线本身的性质——\n",
    "本案例库早期那个数字已经作废，不要引用。统一到同一 2^16 网格后，\n",
    "Hilbert 的领先幅度从约 2.5 倍回落到约 1.8 倍（3.19 → 1.77）。\n",
    "\n",
    "公平方案：\n",
    "\n",
    "1. **同一量化层级**：两者都用 L = 16（或都用 L = 10），量化公式完全一致；\n",
    "2. **同一查询集**：同一批随机小窗口，同样的\"命中数 ≥ 5 才计入\"过滤；\n",
    "3. **同一统计口径**：同样的候选数计算方式（`searchsorted` 的 `side` 一致）、同样的中位数而非均值；\n",
    "4. 若 Hilbert 的逐点循环太慢，可**抽样**再做，但两侧必须抽**同一批点**。\n",
    "\n",
    "一句话：**比较两条曲线之前，先确认它们站在同一把尺子上。**\n",
    "\n",
    "</details>\n",
    "\n",
    "**7. Z 序在真实系统里有哪些落地形态？请说明 Geohash、Google S2、Uber H3 各自的编码方式，以及 H3 每层有多少个五边形。**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "- **Geohash**：基于 **Z 序**，把经纬度交替二分后编成 base32 字符串；\n",
    "  前缀相同表示区域相邻，常用作前缀匹配的粗过滤。\n",
    "- **Google S2**：把球面（立方体投影后的六个面）用 **六条 Hilbert 曲线首尾相连**覆盖，\n",
    "  层级 level 0–30，单元 ID 是 64 位整数。\n",
    "- **Uber H3**：**六边形**层次网格，分辨率 0–15，总单元数 c(r) = 2 + 120·7^r；\n",
    "  由于六边形无法完美铺满球面，**每层恰有 12 个五边形**（拓扑上不可避免）。\n",
    "\n",
    "取舍：Geohash / S2 的四边形单元便于位运算与前缀查询；H3 的六边形单元到邻居的距离更均匀\n",
    "（每个单元 6 个等距邻居），适合邻域聚合，代价是那 12 个五边形。\n",
    "\n",
    "另外，不走线性化路线的还有 **R 树**（Guttman，SIGMOD 1984，用最小外接矩形组织空间对象），\n",
    "以及 PostGIS 的 GiST / SP-GiST 索引、GeoParquet 列式格式、Apache Sedona（Spark 上的空间计算）。\n",
    "\n",
    "</details>\n",
    "\n",
    "**8. 在 PostGIS 里想查\"距离某点 500 米以内的餐厅\"，为什么应该写 `ST_DWithin` 而不是 `ST_Distance < 500`？**\n",
    "\n",
    "<details><summary>答案要点</summary>\n",
    "\n",
    "`ST_DWithin(geom1, geom2, d)` 是**索引感知**（index-aware）的谓词：\n",
    "PostGIS 会先用 GiST 索引的包围盒做粗过滤，再算精确距离，因此**能走索引**。\n",
    "\n",
    "而 `ST_Distance(geom1, geom2) < d` 只是一个普通表达式：优化器不知道它等价于\"距离阈值\"，\n",
    "只能对每一行都算一次距离再比较，**不使用索引**，在大表上退化成全表扫描。\n",
    "\n",
    "```sql\n",
    "-- 推荐：能走 GiST 索引\n",
    "SELECT * FROM restaurants\n",
    "WHERE ST_DWithin(geom, ST_MakePoint(:lon, :lat)::geography, 500);\n",
    "\n",
    "-- 不推荐：全表扫描\n",
    "SELECT * FROM restaurants\n",
    "WHERE ST_Distance(geom, ST_MakePoint(:lon, :lat)::geography) < 500;\n",
    "```\n",
    "\n",
    "（依据：PostGIS 官方文档 Chapter 5。）\n",
    "\n",
    "这个坑与 Z 序实验是同一个道理：**能不能用上索引，取决于谓词的形式，而不是语义。**\n",
    "想让优化器用索引，就要写成它能识别的谓词形式。\n",
    "\n",
    "</details>"
   ]
  },
  {
   "metadata": {},
   "cell_type": "code",
   "source": "",
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}