import json
import math
from pathlib import Path

from PIL import Image, ImageDraw, ImageFont

ROOT = Path(__file__).parent
r = json.loads((ROOT / 'content_statistics_results.json').read_text())
data = json.loads((ROOT / 'analysis_data.json').read_text())
FONT = '/System/Library/Fonts/Supplemental/Arial Unicode.ttf'
im = Image.new('RGB', (1660, 1000), '#F8FAFC')
dr = ImageDraw.Draw(im)


def txt(x, y, value, size=24, color='#334155', anchor=None):
    dr.text((x, y), value, font=ImageFont.truetype(FONT, size), fill=color, anchor=anchor)


txt(65, 38, '视频与图文的观看差异：统计检验与稳健性', 40, '#0F172A')
txt(65, 105, '24篇有效笔记，7篇视频、17篇图文；观察性比较，不代表形式的因果效应', 24, '#64748B')
txt(95, 185, '观看量分布', 29, '#0F172A')
left, right, top, bottom = 150, 690, 265, 755


def yview(v):
    return bottom - math.log10(v) / 5.2 * (bottom - top)


for exponent in range(6):
    yy = yview(10**exponent)
    dr.line((left, yy, right, yy), fill='#E2E8F0', width=2)
    txt(left - 15, yy, f'{10**exponent:,}', 21, '#64748B', 'rm')
dr.line((left, top, left, bottom, right, bottom), fill='#94A3B8', width=2)
for name, flag, x, color in [('图文', 0, 300, '#2563EB'), ('视频', 1, 570, '#D97706')]:
    values = sorted([p for p in data if p['video'] == flag], key=lambda p: p['观看量'])
    for i, p in enumerate(values):
        xx = x + (((i * 7) % 11) - 5) * 7
        yy = yview(p['观看量'])
        dr.ellipse((xx - 7, yy - 7, xx + 7, yy + 7), fill=color, outline='white', width=2)
        if p['source_row'] == 9:
            txt(xx, yy - 18, '77,822', 21, '#92400E', 'mb')
    median = r['genre']['观看量']['median_a' if flag else 'median_b']
    yy = yview(median)
    dr.line((x - 55, yy, x + 55, yy), fill=color, width=5)
    txt(x, bottom + 23, f'{name} n={len(values)}', 24, '#334155', 'mt')
    txt(x, bottom + 63, f'中位数 {median:,.0f}', 22, color, 'mt')
txt(95, 860, '双侧精确秩和检验；3项指标Holm校正', 22)
txt(95, 900, '观看量 p=0.00113；Cliff’s δ=0.866', 24, '#0F172A')

txt(820, 185, '控制发布月份后的观看量比', 29, '#0F172A')
txt(820, 230, '点为几何均值比，横线为HC3稳健95%置信区间', 21, '#64748B')
fleft, fright = 1035, 1570


def xratio(v):
    return fleft + math.log2(v) / 8 * (fright - fleft)


for tick in [1, 4, 16, 64, 256]:
    xx = xratio(tick)
    dr.line((xx, 300, xx, 720), fill='#CBD5E1' if tick == 1 else '#E2E8F0', width=2)
    txt(xx, 742, f'{tick}倍', 21, '#64748B', 'mt')
for i, (key, label) in enumerate([
    ('all_观看量_format_cohort', '全部24篇'),
    ('without_top_观看量_format_cohort', '去掉最高流量'),
    ('september_观看量_format_cohort', '仅9月18篇')
]):
    yy = 355 + i * 145
    term = r['regressions'][key]['terms']['video']
    ratio = math.exp(term['coef'])
    lo, hi = [math.exp(v) for v in term['hc3_ci95_t_approx']]
    txt(820, yy - 14, label, 22)
    dr.line((xratio(lo), yy, xratio(hi), yy), fill='#0F766E', width=4)
    for end in [lo, hi]:
        dr.line((xratio(end), yy - 8, xratio(end), yy + 8), fill='#0F766E', width=3)
    xx = xratio(ratio)
    dr.ellipse((xx - 8, yy - 8, xx + 8, yy + 8), fill='#0F766E')
    txt(xx, yy + 23, f'{ratio:.1f} [{lo:.1f}, {hi:.1f}]', 22, '#0F766E', 'mt')
txt(820, 860, '题材与形式高度重叠，不能稳定分离各自作用。', 22)
txt(820, 900, '视频内仅1篇非展示内容，题材比较 p=0.286。', 22, '#0F172A')
txt(65, 962, '数据来源：笔记列表明细表.xlsx；观看量轴、倍数轴均为对数刻度。', 20, '#64748B')
im.save(ROOT / '内容统计检验.png')

lines = [
    '# 内容与流量的统计分析', '',
    '源文件：/Users/zzyang/Downloads/笔记列表明细表.xlsx，Sheet1。原始26篇，排除第3、4行曝光为0的最新笔记，分析24篇。', '',
    '主要结果变量为观看量，曝光与封面点击率为辅助结果。所有结果均为单账号观察性、探索性分析。题材标签为根据标题进行的事后人工归类，未读取正文或封面。', '',
    '## 方法', '',
    '- 形式比较：双侧精确置换秩和检验，以秩和偏离零假设期望的绝对值为统计量，枚举346,104种分组；观看、曝光、点击率3项p值作Holm校正。检验反映分布/秩差异，并非在任意分布下只检验中位数。',
    '- 效应量：中位数之比、Cliff’s delta、随机一篇视频观看高于随机一篇图文的样本配对比例。中位数比及delta的区间为组内重抽样50,000次的百分位Bootstrap区间。',
    '- 题材总体比较：含并列秩修正的Kruskal–Wallis H统计量，200,000次标签置换估计p值，p=(极端次数+1)/(置换次数+1)。3项结果作Holm校正；观看量的10组两两比较另作Holm校正。',
    '- 回归：ln(观看量)=截距+视频指示变量+10月发布指示变量。9月为参照。HC3异方差稳健标准误，使用残差自由度t分布近似计算p值及95%区间。系数指数化报告条件几何均值比。回归p值未额外作多重校正。',
    '- 稳健性：去掉第9行最高流量笔记、仅分析9月18篇、在9月与10月内分别固定视频篇数作分层精确置换（61,200种分组）。',
    '- 预测评估：以帽子矩阵计算留一交叉验证PRESS及log尺度R²，区别于样本内R²。',
    '- 随机种子20261017。检验将笔记近似视为独立、可交换的观察；同账号连续发布可能不满足严格独立性。', '',
    '## 图文和视频', '',
    '|结果|视频中位数|图文中位数|U|原始p|Holm校正p|Cliff delta|',
    '|---|---:|---:|---:|---:|---:|---:|'
]
for outcome in ['观看量', '曝光', '封面点击率']:
    s = r['genre'][outcome]
    a = f'{s["median_a"]:.1%}' if outcome == '封面点击率' else f'{s["median_a"]:,.0f}'
    b = f'{s["median_b"]:.1%}' if outcome == '封面点击率' else f'{s["median_b"]:,.0f}'
    lines.append(f'|{outcome}|{a}|{b}|{s["U"]:.1f}|{s["p_exact"]:.6f}|{s["p_holm_3"]:.6f}|{s["cliff_delta"]:.3f}|')
s = r['genre']['观看量']
lines += ['', f'观看量中位数比={s["median_ratio"]:.2f}，95% Bootstrap区间={s["median_ratio_ci95_percentile"]}。Cliff delta区间={s["cliff_delta_ci95_percentile"]}。视频高于图文的样本配对比例={s["superiority_probability"]:.2%}，该比例不是未来成功概率。', '',
          '## 内容题材', '', '|题材|篇数|视频篇数|观看量中位数|曝光中位数|点击率中位数|', '|---|---:|---:|---:|---:|---:|']
for s in r['topic_summary']:
    lines.append(f'|{s["group"]}|{s["n"]}|{s["video_n"]}|{s["median_views"]:,.0f}|{s["median_impressions"]:,.1f}|{s["median_ctr"]:.1%}|')
lines += ['', '|结果|H|置换p|Holm校正p|epsilon²|', '|---|---:|---:|---:|---:|']
for outcome in ['观看量', '曝光', '封面点击率']:
    s = r['topic_omnibus'][outcome]
    lines.append(f'|{outcome}|{s["H"]:.3f}|{s["p_permutation"]:.6f}|{s["p_holm_3"]:.6f}|{s["epsilon_squared"]:.3f}|')
lines += ['', 'epsilon²为秩统计量效应量，不能解读为解释了原始观看量的相同比例方差。', '',
          '### 观看量的全部10组两两比较', '', '|A|B|原始p|Holm校正p|Cliff delta|去掉最高流量后校正p|', '|---|---|---:|---:|---:|---:|']
for s, t in zip(r['topic_pairwise_views'], r['topic_pairwise_views_without_top']):
    lines.append(f'|{s["a"]}|{s["b"]}|{s["p_exact"]:.6f}|{s["p_holm_10"]:.6f}|{s["cliff_delta"]:.3f}|{t["p_holm_10"]:.6f}|')
lines += ['', '## 回归及稳健性', '', '|样本|结果|视频系数|几何均值比|HC3 95%区间|稳健p|样本内R²|留一R²(log尺度)|', '|---|---|---:|---:|---|---:|---:|---:|']
for key, m in r['regressions'].items():
    s = m['terms']['video']
    lo, hi = [math.exp(a) for a in s['hc3_ci95_t_approx']]
    lines.append(f'|{key}|log结果|{s["coef"]:.4f}|{math.exp(s["coef"]):.2f}|[{lo:.2f}, {hi:.2f}]|{s["hc3_p_t_approx"]:.6f}|{m["r2"]:.3f}|{m["loo_r2_log_target"]:.3f}|')
lines += ['', '全样本观看量方程：ln(观看量)=4.2505+3.0640×视频−0.1126×10月。月份控制只粗略区分两个发布批次，不能替代统一发布后观察窗口。', '',
          f'按月份分层的精确置换检验，观看量p={r["genre"]["views_month_stratified_permutation"]["p_exact"]:.6f}。', '',
          '## 能下的结论与限制', '',
          '1. 视频与图文的观看分布存在较大差异，去掉最高流量笔记及仅保留9月笔记后仍存在。CTR差异未检出，不能据此证明两种形式CTR等同。',
          '2. 题材总体存在差异。展示类对模型动态类的校正p从0.039变为去掉最高流量后的0.071，具体题材结论比形式比较更脆弱。',
          '3. 展示类6篇全部为视频，7篇视频中仅1篇非展示类。视频内的题材秩和检验p=0.286，只有1篇参照，统计功效不足。不能稳定分离题材和形式各自的作用。',
          '4. 24篇不足以支持包含大量内容标签及交互项的可靠选题模型。观看量发布前信息模型的留一R²为0.230，预测精度有限，且留一法不能充分检验未来时间上的泛化。',
          '5. 所有题材分组及后续敏感性分析为探索性分析，未覆盖所有可能尝试的分组方式。Holm校正仅对应明确列出的比较族，不会消除事后选标签带来的偏差。',
          '6. 未知统一统计窗口、投流、账号粉丝基数、受众变化及平台分发状态。回归倍数为条件关联，不是把图文换成视频的因果增幅。', '',
          '## 题材标签与源行号', '']
for summary in r['topic_summary']:
    matching = [x for x in data if (x['source_row'] in {'技术作品/工具效果展示':[7,9,20,24,25,27], 'Vibe Coding指南':[11,15,16], '模型动态/观点':[5,6,8,17,18], '安全/隐私事件':[14,19,21], '其他':[10,12,13,22,23,26,28]}[summary['group']])]
    lines.append(f'### {summary["group"]}')
    lines += ['', *[f'- 第{x["source_row"]}行：{x["笔记标题"]}' for x in matching], '']
(ROOT / '内容统计分析报告.md').write_text('\n'.join(lines))