iomgaa
60c21ae346
rules 是纯计算，src/rules 共 1761 行：方法.md 第六部分的全部公式、第七部分的余额与扣次、4.2 的入口分配，以及种子与规则登记表。 预期值有三个来源：方法.md 的手算例；tests/fixtures/rules/ 下不看实现写的 Python 脚本；参考项目 calc15、calc14b 里与现行方法一致的部分。审阅做了两轮，每轮都换一个全新上下文的子 Agent，提出的问题都已修复或写明不修的理由。 用户确认的约定（原话"1–8 都同意"）：入池、采纳、消融中发现都用严格不等号；成本用 max(0, ĝ₀)；没测的格附"含它的最好组合"的净值；战绩在分数原值上判。三个数改为严格递增后，删掉了尖峰的全部专门处理（bfb96ba）。vocab 里 Range 的说明同步改为严格递增。 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
785 lines
28 KiB
Python
785 lines
28 KiB
Python
"""rules 预期值的独立参考实现（共享库，本身不生成文件）。
|
||
|
||
用途：被同目录的 gen_forward_info.py、gen_ablation.py 导入，算出预期值写成 JSON，
|
||
供 TypeScript 的 rules 测试读取。只用 Python 标准库。
|
||
|
||
独立性：本文件只照 docs/方法.md 第二部分、4.2、第六、七部分的文字，
|
||
以及用户确认的约定表（scratchpad/conventions.md，下称"约定 N"）写成，
|
||
没有读过 src/rules/ 和参考项目的实现代码。
|
||
|
||
算法上刻意和"照约定逐段积分 1−∏F"的写法分开：
|
||
- 分布全部用有理数（fractions.Fraction，浮点输入先精确转成有理数）表示；
|
||
- 期望一律用 Stieltjes 积分 E[g(M)] = ∫ g(t)·F_M'(t)dt 精确求，
|
||
其中 M 是若干独立变量的最大值，F_M = ∏F（不是 ∫[1−∏F]dt 那条公式）；
|
||
- 最小值通过取负号转成最大值；
|
||
- 部分格信息价值的网格算法是方法规定的确定近似（6.7、约定 27），只能照规定做，用浮点。
|
||
|
||
记号：一条"点列"是 [(x, y), ...]，x 严格递增、y 从 0 到 1 不减。三个数严格递增（方法.md 第五部分），
|
||
折线连续，"低于 x"与"不高于 x"的把握相同。
|
||
"""
|
||
|
||
from fractions import Fraction as Fr
|
||
from itertools import combinations
|
||
import bisect
|
||
import functools
|
||
import math
|
||
import re
|
||
|
||
# ---------------------------------------------------------------- 基础
|
||
|
||
def fr(x):
|
||
"""浮点或整数精确转成有理数。"""
|
||
return x if isinstance(x, Fr) else Fr(x)
|
||
|
||
|
||
def to_float(x):
|
||
return float(x)
|
||
|
||
|
||
def polyline(q1, q5, q9):
|
||
"""6.2：三个数 → 五点折线（有理数）。"""
|
||
q1, q5, q9 = fr(q1), fr(q5), fr(q9)
|
||
if not (q1 < q5 < q9):
|
||
raise ValueError("三个数必须严格递增 q1 < q5 < q9（方法.md 第五部分）")
|
||
return [(q1 - (q5 - q1), Fr(0)), (q1, Fr(1, 10)), (q5, Fr(1, 2)),
|
||
(q9, Fr(9, 10)), (q9 + (q9 - q5), Fr(1))]
|
||
|
||
|
||
def uniform_pts(lo, hi):
|
||
"""[lo, hi] 上均匀分布的点列。"""
|
||
return [(fr(lo), Fr(0)), (fr(hi), Fr(1))]
|
||
|
||
|
||
def reflect(pts, b):
|
||
"""X → b − X 的点列：P(b−X < y) = 1 − P(X ≤ b−y)。点反序，横坐标 b−x，纵坐标 1−y。"""
|
||
b = fr(b)
|
||
return [(b - x, 1 - y) for (x, y) in reversed(pts)]
|
||
|
||
|
||
def shift(pts, b):
|
||
"""X → X − b。"""
|
||
b = fr(b)
|
||
return [(x - b, y) for (x, y) in pts]
|
||
|
||
|
||
def negate(pts):
|
||
"""X → −X。"""
|
||
return reflect(pts, 0)
|
||
|
||
|
||
# 单条点列"低于 x"的把握 F(x)。
|
||
def pts_left(pts, x):
|
||
xs = [p[0] for p in pts]
|
||
i = bisect.bisect_left(xs, x) - 1 # 最后一个 x_i < x
|
||
if i < 0:
|
||
return Fr(0)
|
||
if i == len(pts) - 1:
|
||
return Fr(1)
|
||
(x0, y0), (x1, y1) = pts[i], pts[i + 1]
|
||
return y0 + (y1 - y0) * (x - x0) / (x1 - x0) # x0 < x ≤ x1
|
||
|
||
|
||
def pts_density(pts, x):
|
||
"""约定 4：所在段斜率，段左闭右开（落在折点上取右边那段）；范围外 0。"""
|
||
x = fr(x)
|
||
for (x0, y0), (x1, y1) in zip(pts, pts[1:]):
|
||
if x0 <= x < x1:
|
||
return (y1 - y0) / (x1 - x0)
|
||
return Fr(0)
|
||
|
||
|
||
class Dist:
|
||
"""若干点列的加权混合（权重和为 1）。"""
|
||
|
||
def __init__(self, comps):
|
||
self.comps = [(fr(w), p) for (w, p) in comps if fr(w) != 0]
|
||
|
||
def left(self, x):
|
||
x = fr(x)
|
||
return sum((w * pts_left(p, x) for w, p in self.comps), Fr(0))
|
||
|
||
def breakpoints(self):
|
||
s = set()
|
||
for _, p in self.comps:
|
||
for x, _ in p:
|
||
s.add(x)
|
||
return s
|
||
|
||
def lo(self):
|
||
return min(self.breakpoints())
|
||
|
||
def hi(self):
|
||
return max(self.breakpoints())
|
||
|
||
|
||
def single(pts):
|
||
return Dist([(1, pts)])
|
||
|
||
|
||
# ---------------------------------------------------------------- 多项式（有理系数，低次在前）
|
||
|
||
def pmul(a, b):
|
||
out = [Fr(0)] * (len(a) + len(b) - 1)
|
||
for i, ai in enumerate(a):
|
||
if ai == 0:
|
||
continue
|
||
for j, bj in enumerate(b):
|
||
out[i + j] += ai * bj
|
||
return out
|
||
|
||
|
||
def pderiv(a):
|
||
return [i * a[i] for i in range(1, len(a))] or [Fr(0)]
|
||
|
||
|
||
def pint(a, lo, hi):
|
||
"""∫_lo^hi a(t) dt。"""
|
||
s = Fr(0)
|
||
for i, c in enumerate(a):
|
||
s += c * (hi ** (i + 1) - lo ** (i + 1)) / (i + 1)
|
||
return s
|
||
|
||
|
||
def expect_max(dists, floor0=False):
|
||
"""E[max_i X_i]（各 X_i 独立），floor0=True 时求 E[max(0, max_i X_i)]。
|
||
|
||
Stieltjes 精确积分：F_M = ∏F_i 连续；相邻折点之间各 F_i 都是一次函数，∏ 是多项式，对 t·(∏)' 精确积分。
|
||
floor0 时 0 以下的把握都算在 0 上，贡献 0。
|
||
"""
|
||
bps = sorted(set().union(*[d.breakpoints() for d in dists]))
|
||
total = Fr(0)
|
||
for a, b in zip(bps, bps[1:]):
|
||
if floor0 and b <= 0:
|
||
continue
|
||
# 用平移后的变量 u = t − a 降低有理数的大小
|
||
w = b - a
|
||
poly = [Fr(1)]
|
||
for d in dists:
|
||
fa, fb = d.left(a), d.left(b)
|
||
poly = pmul(poly, [fa, (fb - fa) / w])
|
||
dp = pderiv(poly)
|
||
lo = Fr(0)
|
||
if floor0 and a < 0:
|
||
lo = -a
|
||
# t = u + a
|
||
integrand = pmul([a, Fr(1)], dp)
|
||
total += pint(integrand, lo, w)
|
||
return total
|
||
|
||
|
||
def expect_min(dists):
|
||
"""E[min_i X_i] = −E[max_i (−X_i)]。"""
|
||
neg = [Dist([(w, negate(p)) for w, p in d.comps]) for d in dists]
|
||
return -expect_max(neg)
|
||
|
||
|
||
def mean(d):
|
||
return expect_max([d])
|
||
|
||
|
||
def quantile(d, p):
|
||
"""约定 8：inf{x : F(x) ≥ p}，精确求。F 连续、逐段一次；几条折线之间有空当时 F 走平，取平段左端。"""
|
||
p = fr(p)
|
||
bps = sorted(d.breakpoints())
|
||
for i, a in enumerate(bps):
|
||
if d.left(a) >= p:
|
||
return a
|
||
b = bps[i + 1]
|
||
fa, fb = d.left(a), d.left(b)
|
||
if fb >= p:
|
||
return a + (p - fa) * (b - a) / (fb - fa)
|
||
return bps[-1]
|
||
|
||
|
||
# ---------------------------------------------------------------- 方向与合并
|
||
|
||
def improvement_pts(raw_pts, b, direction):
|
||
"""约定 5：分数原值的点列 → 改进量 Y 的点列。"""
|
||
if direction == "maximize":
|
||
return shift(raw_pts, b)
|
||
if direction == "minimize":
|
||
return reflect(raw_pts, b)
|
||
raise ValueError(direction)
|
||
|
||
|
||
def improve(x, ref, direction):
|
||
"""x 相对 ref 的改进量（按方向）：越大越好时 x − ref。"""
|
||
x, ref = fr(x), fr(ref)
|
||
return x - ref if direction == "maximize" else ref - x
|
||
|
||
|
||
def reserve_bounds(pts_list):
|
||
"""6.3 / 约定 6：全部点列的最低首点 L、最高末点 U，开放预留铺在 [L−W/2, U+W/2]。"""
|
||
L = min(p[0][0] for p in pts_list)
|
||
U = max(p[-1][0] for p in pts_list)
|
||
W = U - L
|
||
return L - W / 2, U + W / 2
|
||
|
||
|
||
def forward_dists(cands, b, direction, lam):
|
||
"""正向合并：cands = {候选编号: [每位研究者的三个数]}，返回 (改进量上的合并分布, 原值上的合并分布, 预留端点)。
|
||
|
||
每人一份、等权；开放预留取本次选批全部候选、全部人的折线（含池外候选）。
|
||
"""
|
||
lam = fr(lam)
|
||
raw = {c: [polyline(*t) for t in triples] for c, triples in cands.items()}
|
||
allraw = [p for ps in raw.values() for p in ps]
|
||
rlo, rhi = reserve_bounds(allraw)
|
||
ylist = [improvement_pts(p, b, direction) for p in allraw]
|
||
ylo, yhi = reserve_bounds(ylist)
|
||
out_y, out_raw = {}, {}
|
||
for c, ps in raw.items():
|
||
n = len(ps)
|
||
out_y[c] = Dist([((1 - lam) / n, improvement_pts(p, b, direction)) for p in ps]
|
||
+ [(lam, uniform_pts(ylo, yhi))])
|
||
out_raw[c] = Dist([((1 - lam) / n, p) for p in ps] + [(lam, uniform_pts(rlo, rhi))])
|
||
return out_y, out_raw, (ylo, yhi), (rlo, rhi)
|
||
|
||
|
||
def explanation_weights(authors, lam):
|
||
"""6.7 / 约定 24：解释合起来 1−λ，研究者等权，同一人几条平分。返回与 authors 同序的权重。"""
|
||
lam = fr(lam)
|
||
n_auth = len(set(authors))
|
||
cnt = {a: authors.count(a) for a in set(authors)}
|
||
return [(1 - lam) / n_auth / cnt[a] for a in authors]
|
||
|
||
|
||
def cell_dist(triples, authors, lam):
|
||
"""消融一格的成分混合（原值）：各解释折线 + 这一格自己的开放预留。"""
|
||
pts = [polyline(*t) for t in triples]
|
||
lo, hi = reserve_bounds(pts)
|
||
ws = explanation_weights(authors, lam)
|
||
return Dist([(w, p) for w, p in zip(ws, pts)] + [(fr(lam), uniform_pts(lo, hi))]), (lo, hi)
|
||
|
||
|
||
# ---------------------------------------------------------------- 正向：入池、V、选批、读数
|
||
|
||
def pool_prob(dy, eps):
|
||
"""约定 9：P(Y < −ε)。"""
|
||
return dy.left(-fr(eps))
|
||
|
||
|
||
def V(dys):
|
||
"""6.4：批内最好期望 E[max(0, max_c Y_c)]。空集为 0。"""
|
||
if not dys:
|
||
return Fr(0)
|
||
return expect_max(dys, floor0=True)
|
||
|
||
|
||
def _str_cmp(a, b):
|
||
"""JS 默认字符串比较（UTF-16 码元）。"""
|
||
x, y = a.encode("utf-16-be"), b.encode("utf-16-be")
|
||
return -1 if x < y else (1 if x > y else 0)
|
||
|
||
|
||
def _seg_cmp(x, y):
|
||
if not (x[0].isdigit() and y[0].isdigit()):
|
||
return _str_cmp(x, y)
|
||
tx, ty = x.lstrip("0") or "0", y.lstrip("0") or "0"
|
||
if len(tx) != len(ty):
|
||
return -1 if len(tx) < len(ty) else 1
|
||
return _str_cmp(tx, ty)
|
||
|
||
|
||
def cmp_ids(a, b):
|
||
"""编号的自然序（方法"编号最小"指数字意义上的大小）：c2 < c10，researcher-2 < researcher-10，p2 < p10。
|
||
|
||
编号切成数字段和非数字段逐段比：都是数字段时先比去掉前导零后的位数、再比字符串；否则按字符串比；
|
||
段先用完的排前面；逐段都相同时按整串字符串比。候选、研究者、入口、因子、解释编号都用它。
|
||
"""
|
||
sa, sb = re.findall(r"\d+|\D+", a), re.findall(r"\d+|\D+", b)
|
||
for x, y in zip(sa, sb):
|
||
c = _seg_cmp(x, y)
|
||
if c:
|
||
return c
|
||
if len(sa) != len(sb):
|
||
return -1 if len(sa) < len(sb) else 1
|
||
return _str_cmp(a, b)
|
||
|
||
|
||
id_key = functools.cmp_to_key(cmp_ids)
|
||
|
||
|
||
def combo_key(combo):
|
||
"""组合的先后：格数少的在前，再按编号逐个比。"""
|
||
return (len(combo), [id_key(c) for c in combo])
|
||
|
||
|
||
def select_batch(dys, pool, k):
|
||
"""约定 11：枚举池中所有 min(k,|池|) 组合，V 最大；V 完全相等时按排序后编号逐个比（自然序），小者优先。"""
|
||
size = min(k, len(pool))
|
||
ids = sorted(pool, key=id_key)
|
||
table = []
|
||
for combo in combinations(ids, size):
|
||
table.append((list(combo), V([dys[c] for c in combo])))
|
||
best = None
|
||
for combo, v in table:
|
||
if best is None or v > best[1]:
|
||
best = (combo, v)
|
||
elif v == best[1] and combo_key(combo) < combo_key(best[0]):
|
||
best = (combo, v)
|
||
return table, best
|
||
|
||
|
||
def adoption_line(eps, B, Kstar, direction):
|
||
"""6.5 / 约定 15。"""
|
||
eps = fr(eps)
|
||
if Kstar is None:
|
||
return eps
|
||
room = improve(Kstar, B, direction) # 越大越好时 K* − B
|
||
line = min(eps, room / 4)
|
||
return max(line, eps / 100)
|
||
|
||
|
||
# ---------------------------------------------------------------- 消融结算
|
||
|
||
def supports(B, full, offs, direction):
|
||
"""6.6 / 约定 17：撑(p) = improvement(关掉 p, 完整)/G；G ≤ 0 抛错。offs = {因子: 分数或 None}。"""
|
||
G = improve(full, B, direction)
|
||
if G <= 0:
|
||
raise ValueError("G ≤ 0")
|
||
return G, {p: (None if x is None else improve(full, x, direction) / G) for p, x in offs.items()}
|
||
|
||
|
||
def label(sup):
|
||
"""6.6 表，从上往下。sup = {因子: 撑 或 None(没有有效分数)}。"""
|
||
vals = [v for v in sup.values() if v is not None]
|
||
if any(v < 0 for v in vals):
|
||
return "negative", None
|
||
if any(v is None for v in sup.values()):
|
||
return "incomplete", None
|
||
i = sum(vals) - 1
|
||
if all(v > Fr(9, 10) for v in vals):
|
||
return "indispensable", i
|
||
if i > Fr(1, 5):
|
||
return "complementary", i
|
||
if i < -Fr(1, 5):
|
||
return "substitute", i
|
||
return "independent", i
|
||
|
||
|
||
def factor_key(p):
|
||
"""因子编号也按自然序。"""
|
||
return id_key(p)
|
||
|
||
|
||
# ---------------------------------------------------------------- 信息价值
|
||
|
||
class AblationModel:
|
||
"""一次消融的成分模型。cells: 因子编号列表（p1…）；expl: [{id, author, cells: {p: (q1,q5,q9)}}]。"""
|
||
|
||
def __init__(self, B, full, Kstar, eps, lam, direction, cells, expl):
|
||
self.B, self.full, self.eps, self.lam, self.dir = fr(B), fr(full), fr(eps), fr(lam), direction
|
||
self.Kstar = None if Kstar is None else fr(Kstar)
|
||
self.cells = sorted(cells, key=factor_key)
|
||
self.G = improve(self.full, self.B, direction)
|
||
if self.G <= 0:
|
||
raise ValueError("G ≤ 0")
|
||
self.h = Fr(1) if self.Kstar is None else max(improve(self.Kstar, self.full, direction), self.eps / 100)
|
||
authors = [e["author"] for e in expl]
|
||
self.weights = explanation_weights(authors, lam) + [self.lam]
|
||
self.names = [e["id"] for e in expl] + ["reserve"]
|
||
# comp_pts[j][p]：成分 j 在格 p 的点列
|
||
self.comp_pts = [{p: polyline(*e["cells"][p]) for p in self.cells} for e in expl]
|
||
self.triples = {p: [tuple(float(q) for q in e["cells"][p]) for e in expl] for p in self.cells}
|
||
res = {}
|
||
self.reserve = {}
|
||
for p in self.cells:
|
||
lo, hi = reserve_bounds([c[p] for c in self.comp_pts])
|
||
self.reserve[p] = (lo, hi)
|
||
res[p] = uniform_pts(lo, hi)
|
||
self.comp_pts.append(res)
|
||
# 成分在每格的期望收益（精确）
|
||
self.Epay = [{p: self.payoff(mean(single(c[p]))) for p in self.cells} for c in self.comp_pts]
|
||
|
||
def payoff(self, x):
|
||
"""约定 25：收益 = h·improvement(x, 完整)/G。"""
|
||
return self.h * improve(self.full, x, self.dir) / self.G
|
||
|
||
def prior_payoffs(self, weights=None):
|
||
w = self.weights if weights is None else weights
|
||
return {p: sum(wj * E[p] for wj, E in zip(w, self.Epay)) for p in self.cells}
|
||
|
||
def g0(self):
|
||
return max(self.prior_payoffs().values())
|
||
|
||
def voi_full_exact(self):
|
||
"""全测 VOI，一维精确：Σ_j w_j E_j[max_p 收益_p] − max_p Σ_j w_j E_j[收益_p]。
|
||
|
||
收益是 x 的同一个仿射函数，max 收益对应 min x（越大越好）或 max x（越小越好）。
|
||
"""
|
||
first = Fr(0)
|
||
for wj, c in zip(self.weights, self.comp_pts):
|
||
ds = [single(c[p]) for p in self.cells]
|
||
ext = expect_min(ds) if self.dir == "maximize" else expect_max(ds)
|
||
first += wj * self.payoff(ext)
|
||
return first - self.g0()
|
||
|
||
# ---- 网格（约定 27），浮点
|
||
@staticmethod
|
||
def m_for(k):
|
||
m = int(round(4_000_000 ** (1 / k))) # 先取浮点近似，再用整数运算校正
|
||
while m ** k > 4_000_000:
|
||
m -= 1
|
||
while (m + 1) ** k <= 4_000_000:
|
||
m += 1
|
||
return m
|
||
|
||
def cell_grid(self, p, m):
|
||
"""一格的小段：返回 [(中点, [各成分质量])] 与小段说明。"""
|
||
fpts = self.float_pts(p)
|
||
bps = sorted(set(x for pts in fpts for x, _ in pts))
|
||
nseg = len(bps) - 1
|
||
s = min(20, m // nseg)
|
||
edges = []
|
||
if s >= 1:
|
||
for a, b in zip(bps, bps[1:]):
|
||
for i in range(s):
|
||
edges.append(a + (b - a) * i / s)
|
||
edges.append(bps[-1])
|
||
else:
|
||
a, b = bps[0], bps[-1]
|
||
for i in range(m):
|
||
edges.append(a + (b - a) * i / m)
|
||
edges.append(b)
|
||
out = []
|
||
n = len(edges) - 1
|
||
for i in range(n):
|
||
lo, hi = edges[i], edges[i + 1]
|
||
mass = []
|
||
for fp in fpts:
|
||
mass.append(f_left(fp, hi) - f_left(fp, lo))
|
||
out.append(((lo + hi) / 2, mass))
|
||
return out, {"breakpoints": len(bps), "segments": nseg, "s": s, "subsegments": n}
|
||
|
||
def float_pts(self, p):
|
||
"""网格用的浮点点列：照方法原式在浮点上算尾巴 q1−(q5−q1)、q9+(q9−q5) 与预留 L−W/2、U+W/2。
|
||
|
||
网格的折点个数决定 s，横坐标差一个 ulp 就可能多出或少掉一段（例如 1.4−(1.75−1.4) 得到
|
||
1.0499999999999998，与另一条折线的 1.05 并存），所以这里必须和实现在同一套浮点式子上取折点。
|
||
"""
|
||
out = []
|
||
for (q1, q5, q9) in self.triples[p]:
|
||
out.append([(q1 - (q5 - q1), 0.0), (q1, 0.1), (q5, 0.5), (q9, 0.9), (q9 + (q9 - q5), 1.0)])
|
||
L = min(o[0][0] for o in out)
|
||
U = max(o[-1][0] for o in out)
|
||
W = U - L
|
||
out.append([(L - W / 2, 0.0), (U + W / 2, 1.0)])
|
||
return out
|
||
|
||
def voi_grid(self, S):
|
||
"""约定 27 的网格 VOI（浮点）。S：测到的格（因子编号列表）。"""
|
||
S = sorted(S, key=factor_key)
|
||
k = len(S)
|
||
m = self.m_for(k)
|
||
grids, info = [], {}
|
||
for p in S:
|
||
g, inf = self.cell_grid(p, m)
|
||
grids.append(g)
|
||
info[p] = inf
|
||
J = len(self.weights)
|
||
w = [float(x) for x in self.weights]
|
||
unmeas = [p for p in self.cells if p not in S]
|
||
Eun = [[float(self.Epay[j][q]) for j in range(J)] for q in unmeas]
|
||
pay = [[float(self.payoff(fr(mid))) for mid, _ in g] for g in grids]
|
||
total = 0.0
|
||
# 逐层展开组合：partial = (各成分联合质量向量, 已测格中点收益的最大值)
|
||
def rec(level, joint, bestpay):
|
||
nonlocal total
|
||
g = grids[level]
|
||
pl = pay[level]
|
||
if level == k - 1:
|
||
for idx, (_, mass) in enumerate(g):
|
||
jv = [joint[j] * mass[j] for j in range(J)]
|
||
marg = sum(jv)
|
||
if marg == 0:
|
||
continue
|
||
best = max(bestpay, pl[idx])
|
||
for E in Eun:
|
||
v = sum(jv[j] * E[j] for j in range(J)) / marg
|
||
if v > best:
|
||
best = v
|
||
total += best * marg
|
||
else:
|
||
for idx, (_, mass) in enumerate(g):
|
||
jv = [joint[j] * mass[j] for j in range(J)]
|
||
if not any(jv):
|
||
continue
|
||
rec(level + 1, jv, max(bestpay, pl[idx]))
|
||
rec(0, w, -math.inf)
|
||
prior = max(float(v) for v in self.prior_payoffs().values())
|
||
return total - prior, {"k": k, "m": m, "cells": info}
|
||
|
||
def voi_partial_exact_k1(self, p):
|
||
"""只测一格时的精确参考值（不是方法规定的算法，只用来估网格误差）。
|
||
|
||
相邻折点之间各成分密度是常数，后验权重不随 x 变，没测格的估计收益是常数 c；
|
||
测到的格收益是 x 的一次函数 A(x)。∫ max(A(x), c)·marg(x) dx 精确求。
|
||
"""
|
||
J = len(self.weights)
|
||
others = [q for q in self.cells if q != p]
|
||
bps = sorted(set(x for c in self.comp_pts for x, _ in c[p]))
|
||
total = Fr(0)
|
||
def best_const(wv):
|
||
s = sum(wv)
|
||
return max(sum(wv[j] * self.Epay[j][q] for j in range(J)) / s for q in others) if others else None
|
||
for a, b in zip(bps, bps[1:]):
|
||
wv = [self.weights[j] * pts_density(self.comp_pts[j][p], a) for j in range(J)]
|
||
dens = sum(wv)
|
||
if dens == 0:
|
||
continue
|
||
c = best_const(wv)
|
||
# A(x) = alpha + beta x
|
||
pa, pb = self.payoff(a), self.payoff(b)
|
||
beta = (pb - pa) / (b - a)
|
||
alpha = pa - beta * a
|
||
def integ(lo, hi, poly):
|
||
return dens * pint(poly, lo, hi)
|
||
if c is None:
|
||
total += integ(a, b, [alpha, beta])
|
||
continue
|
||
if beta == 0:
|
||
total += dens * (b - a) * max(alpha, c)
|
||
continue
|
||
xc = (c - alpha) / beta # A(xc) = c
|
||
cuts = [a] + ([xc] if a < xc < b else []) + [b]
|
||
for lo, hi in zip(cuts, cuts[1:]):
|
||
midA = alpha + beta * (lo + hi) / 2
|
||
if midA >= c:
|
||
total += integ(lo, hi, [alpha, beta])
|
||
else:
|
||
total += dens * (hi - lo) * c
|
||
return total - self.g0()
|
||
|
||
|
||
def f_left(fp, x):
|
||
xs = [q[0] for q in fp]
|
||
i = bisect.bisect_left(xs, x) - 1
|
||
if i < 0:
|
||
return 0.0
|
||
if i == len(fp) - 1:
|
||
return 1.0
|
||
(x0, y0), (x1, y1) = fp[i], fp[i + 1]
|
||
return y0 + (y1 - y0) * (x - x0) / (x1 - x0)
|
||
|
||
|
||
# ---------------------------------------------------------------- 价值、成本、选格
|
||
|
||
def m_prime(balance, size, K, r):
|
||
num = balance - size
|
||
if num <= 0:
|
||
return 0, Fr(0)
|
||
R = -(-num // K)
|
||
r = fr(r)
|
||
return R, sum((r ** j for j in range(R)), Fr(0))
|
||
|
||
|
||
def cost(g0, size, K):
|
||
g = max(Fr(0), fr(g0))
|
||
return 2 * g * (Fr(K + size, K + size + 1) - Fr(K, K + 1))
|
||
|
||
|
||
TOL = 1e-12
|
||
|
||
|
||
def top_combo(rows):
|
||
"""6.7 选格第 2 条：最大净值 M；净值 ≥ M − 1e-12 的视为相同，先取格数少的，再按因子编号逐个比取小的。"""
|
||
M = max(n for _, n in rows)
|
||
near = [r for r in rows if r[1] >= M - TOL]
|
||
return min(near, key=lambda r: combo_key(list(r[0])))
|
||
|
||
|
||
def select_cells(cells, balance, net_of):
|
||
"""约定 29 / 方法 6.7 选格：信息价值减成本选格。net_of(组合元组) → 净值（浮点）。
|
||
|
||
先只留净值 > 1e-12 的组合，一个都没有就不测；在留下的组合里按 top_combo 取。
|
||
“含这一格的净值最大组合”只是记录，不受 > 1e-12 的限制，在含这一格的全部组合里按 top_combo 取。
|
||
返回 (选中组合, 选中净值, 各组合净值表, 没选中格的说明)。
|
||
"""
|
||
cells = sorted(cells, key=factor_key)
|
||
maxn = min(len(cells), balance - 1)
|
||
table = []
|
||
for n in range(1, maxn + 1):
|
||
for combo in combinations(cells, n):
|
||
table.append((combo, net_of(combo)))
|
||
best, bestnet = (), 0.0
|
||
eligible = [r for r in table if r[1] > TOL]
|
||
if eligible:
|
||
best, bestnet = top_combo(eligible)
|
||
unselected = {}
|
||
for c in cells:
|
||
if c in best:
|
||
continue
|
||
reason = "insufficientBudget" if len(best) + 1 > balance - 1 else "notWorthCost"
|
||
cand = [(combo, net) for combo, net in table if c in combo]
|
||
if not cand:
|
||
unselected[c] = {"reason": reason, "bestCombo": None, "bestNet": None}
|
||
else:
|
||
bc, bn = top_combo(cand)
|
||
unselected[c] = {"reason": reason, "bestCombo": list(bc), "bestNet": bn}
|
||
return list(best), bestnet, table, unselected
|
||
|
||
|
||
# ---------------------------------------------------------------- 押的因子
|
||
|
||
def bet(model, outcomes):
|
||
"""6.8 / 约定 30。outcomes = {因子: {"kind": ..., "value": 分数或缺}}。"""
|
||
valid = {p: fr(o["value"]) for p, o in outcomes.items() if o["kind"] == "valid"}
|
||
J = len(model.weights)
|
||
raw = []
|
||
for j in range(J):
|
||
v = model.weights[j]
|
||
for p, x in valid.items():
|
||
v *= pts_density(model.comp_pts[j][p], x)
|
||
raw.append(v)
|
||
s = sum(raw)
|
||
updated = s != 0
|
||
w = [v / s for v in raw] if updated else list(model.weights)
|
||
est = model.prior_payoffs(w)
|
||
pay = {}
|
||
for p in model.cells:
|
||
pay[p] = model.payoff(valid[p]) if p in valid else est[p]
|
||
pick = None
|
||
for p in model.cells:
|
||
if pick is None or pay[p] > pay[pick]:
|
||
pick = p
|
||
neg = []
|
||
for p in model.cells:
|
||
if p in valid and improve(model.full, valid[p], model.dir) < 0:
|
||
neg.append({"factor": p, "offBetterBy": float(improve(model.full, valid[p], model.dir) * -1)})
|
||
if neg:
|
||
status = "noBetNegativeSupport"
|
||
elif pay[pick] <= 0:
|
||
status = "noBetNonPositiveGain"
|
||
else:
|
||
status = "bet"
|
||
return {
|
||
"weightsUpdated": updated,
|
||
"weights": {n: float(x) for n, x in zip(model.names, w)},
|
||
"payoffs": {p: float(v) for p, v in pay.items()},
|
||
"picked": pick,
|
||
"predictedGain": float(pay[pick]),
|
||
"status": status,
|
||
"negativeFactors": neg,
|
||
}
|
||
|
||
|
||
# ---------------------------------------------------------------- 信息账
|
||
|
||
def edge(b, delta, j):
|
||
"""约定 31/32：小格边界一律按浮点 b + j·Δ 算（和实现用同一套边界）。"""
|
||
return b + j * delta
|
||
|
||
|
||
def bin_of(v, b, delta):
|
||
j = math.floor((v - b) / delta)
|
||
while edge(b, delta, j) > v:
|
||
j -= 1
|
||
while edge(b, delta, j + 1) <= v:
|
||
j += 1
|
||
return j
|
||
|
||
|
||
def bin_masses(d, b, delta):
|
||
"""覆盖分布全部折点的小格质量（有理数），返回 [(j, p)]。"""
|
||
lo, hi = float(d.lo()), float(d.hi())
|
||
j0 = bin_of(lo, b, delta)
|
||
j1 = bin_of(hi, b, delta)
|
||
# 浮点的 lo/hi 与有理数端点可能差一点；两侧各多切一格，质量为 0 的格不影响熵
|
||
out = []
|
||
prev = d.left(edge(b, delta, j0 - 1))
|
||
for j in range(j0 - 1, j1 + 2):
|
||
cur = d.left(edge(b, delta, j + 1))
|
||
out.append((j, cur - prev))
|
||
prev = cur
|
||
assert d.left(edge(b, delta, j0 - 1)) == 0 and prev == 1, "切格范围没有覆盖全部质量"
|
||
return out
|
||
|
||
|
||
def world_bits(d, b, delta):
|
||
terms = []
|
||
for _, p in bin_masses(d, b, delta):
|
||
if p > 0:
|
||
pf = float(p)
|
||
terms.append(-pf * math.log2(pf))
|
||
return math.fsum(terms)
|
||
|
||
|
||
def surprise(d, v, b, delta):
|
||
j = bin_of(v, b, delta)
|
||
p = d.left(edge(b, delta, j + 1)) - d.left(edge(b, delta, j))
|
||
if p == 0:
|
||
return {"status": "priorZero", "bin": j, "p": 0.0}
|
||
return {"status": "ok", "bin": j, "p": float(p), "surprise": -math.log2(float(p))}
|
||
|
||
|
||
def disagreement(person_dists, b, delta):
|
||
avg = Dist([(Fr(1, len(person_dists)) * w, p) for d in person_dists for w, p in d.comps])
|
||
hb = world_bits(avg, b, delta)
|
||
each = [world_bits(d, b, delta) for d in person_dists]
|
||
return hb - math.fsum(each) / len(each), hb, each
|
||
|
||
|
||
def landing(q1, q5, q9, v):
|
||
"""约定 34。"""
|
||
if v < q1:
|
||
return "below10"
|
||
if v < q5:
|
||
return "from10to50"
|
||
if v <= q9:
|
||
return "from50to90"
|
||
return "above90"
|
||
|
||
|
||
# ---------------------------------------------------------------- 入口分配、余额、种子
|
||
|
||
def allocate_entries(researchers, entries, prefs):
|
||
"""4.2：researchers 研究者编号；entries = [{id, author}]；prefs = {研究者: [入口编号…]}。
|
||
|
||
研究者和入口都先按编号（自然序）排好，“编号最小”按编号判，与传入次序无关。
|
||
"""
|
||
researchers = sorted(researchers, key=id_key)
|
||
entries = sorted(entries, key=lambda e: id_key(e["id"]))
|
||
author = {e["id"]: e["author"] for e in entries}
|
||
taken = {}
|
||
result = {}
|
||
maxlen = max((len(v) for v in prefs.values()), default=0)
|
||
if len(researchers) == 1:
|
||
r = researchers[0]
|
||
return {r: {"entry": prefs[r][0], "choice": 1}}
|
||
for k in range(maxlen):
|
||
want = {}
|
||
for r in researchers:
|
||
if r in result or k >= len(prefs[r]):
|
||
continue
|
||
e = prefs[r][k]
|
||
if e in taken:
|
||
continue
|
||
want.setdefault(e, []).append(r)
|
||
for e in [x["id"] for x in entries]:
|
||
if e not in want:
|
||
continue
|
||
group = want[e]
|
||
if len(group) == 1:
|
||
win = group[0]
|
||
elif author[e] in group:
|
||
win = author[e]
|
||
else:
|
||
win = min(group, key=id_key)
|
||
result[win] = {"entry": e, "choice": k + 1}
|
||
taken[e] = win
|
||
for r in researchers:
|
||
if r in result:
|
||
continue
|
||
left = [x["id"] for x in entries if x["id"] not in taken]
|
||
e = left[0]
|
||
result[r] = {"entry": e, "choice": None}
|
||
taken[e] = r
|
||
return {r: result[r] for r in researchers}
|
||
|
||
|
||
def balance(total, runs):
|
||
started = sum(1 for s in runs if s == "started")
|
||
reserved = sum(1 for s in runs if s == "reserved")
|
||
return {"charged": started, "occupied": started + reserved, "balance": total - started - reserved}
|