import json
import math
import warnings
from pathlib import Path

import numpy as np
import openpyxl
import pandas as pd

warnings.filterwarnings('ignore', message='Workbook contains no default style')
OUT = Path(__file__).parent
SOURCE = '/Users/zzyang/Downloads/笔记列表明细表.xlsx'


def beta_cf(a, b, x):
    qab, qap, qam = a + b, a + 1, a - 1
    c, d = 1.0, 1 - qab * x / qap
    d = 1 / (d if abs(d) > 1e-300 else 1e-300)
    h = d
    for m in range(1, 301):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1 + aa * d
        c = 1 + aa / c
        d = 1 / (d if abs(d) > 1e-300 else 1e-300)
        c = c if abs(c) > 1e-300 else 1e-300
        h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1 + aa * d
        c = 1 + aa / c
        d = 1 / (d if abs(d) > 1e-300 else 1e-300)
        c = c if abs(c) > 1e-300 else 1e-300
        change = d * c
        h *= change
        if abs(change - 1) < 3e-14:
            return h
    raise ArithmeticError('Incomplete beta did not converge')


def ibeta(a, b, x):
    if x <= 0:
        return 0.0
    if x >= 1:
        return 1.0
    bt = math.exp(math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
                  + a * math.log(x) + b * math.log1p(-x))
    if x < (a + 1) / (a + b + 2):
        return bt * beta_cf(a, b, x) / a
    return 1 - bt * beta_cf(b, a, 1 - x) / b


def t_pvalue(t, df):
    return ibeta(df / 2, 0.5, df / (df + t * t))


def t_critical(df):
    lo, hi = 0.0, 100.0
    for _ in range(80):
        mid = (lo + hi) / 2
        if t_pvalue(mid, df) > 0.05:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def ols(y, x, names):
    y, x = np.asarray(y, dtype=float), np.asarray(x, dtype=float)
    if x.ndim == 1:
        x = x[:, None]
    X = np.column_stack([np.ones(len(y)), x])
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    residual = y - X @ beta
    n, k = X.shape
    df = n - k
    inv = np.linalg.inv(X.T @ X)
    sse = residual @ residual
    s2 = sse / df
    se = np.sqrt(np.diag(inv) * s2)
    leverage = np.sum((X @ inv) * X, axis=1)
    scaled = residual**2 / (1 - leverage)**2
    hc3_cov = inv @ (X.T @ (scaled[:, None] * X)) @ inv
    hc3_se = np.sqrt(np.diag(hc3_cov))
    critical = t_critical(df)
    r2 = 1 - sse / np.sum((y - y.mean())**2)
    result = {'n': n, 'df': df, 'r2': float(r2),
              'adj_r2': float(1 - (1 - r2) * (n - 1) / df),
              'residual_se': float(np.sqrt(s2)), 'terms': {}}
    for i, name in enumerate(['intercept'] + names):
        result['terms'][name] = {
            'coef': float(beta[i]), 'se': float(se[i]),
            'p': t_pvalue(beta[i] / se[i], df),
            'ci95': [float(beta[i] - critical * se[i]), float(beta[i] + critical * se[i])],
            'hc3_se': float(hc3_se[i]),
            'hc3_p_t_approx': t_pvalue(beta[i] / hc3_se[i], df),
            'hc3_ci95_t_approx': [float(beta[i] - critical * hc3_se[i]), float(beta[i] + critical * hc3_se[i])]
        }
    return result


assert abs(t_pvalue(1, 1) - .5) < 1e-10
assert abs(t_critical(22) - 2.0738730679) < 1e-8
w = openpyxl.load_workbook(SOURCE, data_only=True)
rows = list(w.active.values)
d = pd.DataFrame(rows[2:], columns=rows[1])
d['source_row'] = np.arange(3, 3 + len(d))
v = d[d['曝光'] > 0].copy()
v['ctr_pp'] = v['封面点击率'] * 100
v['video'] = (v['体裁'] == '视频').astype(int)
v['published'] = pd.to_datetime(v['首次发布时间'], format='%Y年%m月%d日%H时%M分%S秒')
v['publish_day'] = (v['published'] - v['published'].min()).dt.total_seconds() / 86400
v['ln_imp'] = np.log(v['曝光'])
top = v['观看量'].idxmax()
models = {}
for subset_name, sub in [('all', v), ('without_top', v.drop(index=top))]:
    for outcome, label in [('曝光', 'impressions'), ('观看量', 'views')]:
        for transform in ['raw', 'ln']:
            y = sub[outcome] if transform == 'raw' else np.log(sub[outcome])
            models[f'{subset_name}_{label}_{transform}'] = ols(y, sub['ctr_pp'], ['ctr_pp'])
        models[f'{subset_name}_{label}_ln_controls'] = ols(
            np.log(sub[outcome]), sub[['ctr_pp', 'video', 'publish_day']],
            ['ctr_pp', 'video', 'publish_day'])
    models[f'{subset_name}_views_ln_imp_controls'] = ols(
        np.log(sub['观看量']), sub[['ctr_pp', 'ln_imp', 'video']],
        ['ctr_pp', 'ln_imp', 'video'])

result = {'source': SOURCE, 'total_n': len(d), 'excluded_rows': d[d['曝光'] == 0]['source_row'].tolist(),
          'top_title': v.loc[top, '笔记标题'], 'top_source_row': int(v.loc[top, 'source_row']),
          'models': models}
(OUT / 'regression_results.json').write_text(json.dumps(result, ensure_ascii=False, indent=2))
for key, model in models.items():
    t = model['terms']['ctr_pp']
    print(key, json.dumps({'n': model['n'], 'intercept': model['terms']['intercept']['coef'],
          'r2': model['r2'], **t}, ensure_ascii=False))
print('INDEPENDENT_MODEL', json.dumps(models['all_views_ln_imp_controls'], ensure_ascii=False))
v.to_json(OUT / 'analysis_data.json', orient='records', force_ascii=False, date_format='iso')