# 隐藏维度"建造"一半: 制造商分组的残差方差分解
# 制造商为公开登记信息(RCDB/Wikipedia), 按"当前轨道的建造者"归类:
#   RMC 改造车归 RMC; GhostRider 2016 年 GCI 整体换轨归 GCI(原建 CCI)
#   不确定者已线上核对: Jetline=Schwarzkopf(设计,Zierer制造), Cu Chulainn=Gravity Group
# 运行: 在 data/ 目录 python3 class05-bt-mfr.py
import csv, math, random, unicodedata, difflib

random.seed(7)
BASE = '.'
src_code = open(f'{BASE}/class05-bt-rank.py').read()
i0 = src_code.index('STEEL = [')
i1 = src_code.index("(50,'Twister','Knoebels')]") + len("(50,'Twister','Knoebels')]")
exec(src_code[i0:i1])

MFR = {
 # ---- 钢质榜 ----
 ('S','Fury 325'):'B&M', ('S','Millennium Force'):'Intamin', ('S','Steel Vengeance'):'RMC',
 ('S','Expedition GeForce'):'Intamin', ('S','Superman: The Ride'):'Intamin',
 ('S',"Apollo's Chariot"):'B&M', ('S','Iron Rattler'):'RMC', ('S','Leviathan'):'B&M',
 ('S','Maverick'):'Intamin', ('S','Diamondback'):'B&M', ('S','Nitro'):'B&M',
 ('S','Intimidator 305'):'Intamin', ('S',"Phantom's Revenge"):'Morgan',
 ('S','Magnum XL-200'):'Arrow', ('S','Taron'):'Intamin', ('S','Top Thrill Dragster'):'Intamin',
 ('S','Mako'):'B&M', ('S','Time Traveler'):'Mack', ('S','Blue Fire'):'Mack',
 ('S','Nemesis'):'B&M', ('S','Helix'):'Mack', ('S','Intimidator'):'B&M',
 ('S','New Texas Giant'):'RMC', ('S','Twisted Colossus'):'RMC', ('S','Mind Bender'):'Schwarzkopf',
 ('S','Goliath'):'B&M',  # SFOG 与 La Ronde 两台 Goliath 均为 B&M
 ('S','Behemoth'):'B&M', ('S','Montu'):'B&M', ('S','Banshee'):'B&M', ('S','Skyrush'):'Intamin',
 ('S','X2'):'Arrow', ('S','Alpengeist'):'B&M', ('S','Wicked Cyclone'):'RMC',
 ('S','Black Mamba'):'B&M', ('S','Cheetah Hunt'):'Intamin', ('S','Kumba'):'B&M',
 ('S','Twisted Timbers'):'RMC', ('S','Jetline'):'Schwarzkopf', ('S','Superman Ride of Steel'):'Intamin',
 ('S','Lisebergbanan'):'Schwarzkopf', ('S','Griffon'):'B&M', ('S','Cannibal'):'InHouse',
 ('S','Shambhala'):'B&M', ('S','Expedition Everest'):'Vekoma', ('S','Storm Chaser'):'RMC',
 ('S','Raging Bull'):'B&M', ('S','Thunderbird'):'B&M', ('S','Whizzer'):'Schwarzkopf',
 # ---- 木质榜 ----
 ('W','Phoenix'):'PTC', ('W','El Toro'):'Intamin', ('W','Voyage'):'GravityGroup',
 ('W','Beast'):'InHouse', ('W','Outlaw Run'):'RMC', ('W','Ravine Flyer II'):'GravityGroup',
 ('W','Gold Striker'):'GCI', ('W','Thunderhead'):'GCI', ('W','Mystic Timbers'):'GCI',
 ('W','Lightning Racer'):'GCI', ('W','GhostRider'):'GCI', ('W','Balder'):'Intamin',
 ('W','Thunderbolt'):'InHouse', ('W','Wodan'):'GCI', ('W','Wildfire'):'RMC',
 ('W','Raven'):'CCI', ('W','Goliath'):'RMC', ('W','Jack Rabbit'):'Miller',
 ('W','Giant Dipper'):'InHouse', ('W','Shivering Timbers'):'CCI', ('W','Legend'):'CCI',
 ('W','White Lightning'):'GCI', ('W','Troy'):'GCI', ('W','Renegade'):'GCI',
 ('W','Cu Chulainn'):'GravityGroup', ('W','Colossos'):'Intamin', ('W','Cyclone'):'InHouse',
 ('W','Prowler'):'GCI', ('W','Rampage'):'CCI', ('W','Comet'):'PTC',
 ('W','Boardwalk Bullet'):'GravityGroup', ('W','Flying Turns'):'InHouse',
 ('W','Rutschebanan'):'InHouse', ('W','Switchback'):'GravityGroup',
 ('W','American Thunder'):'GCI', ('W',"Screamin' Eagle"):'PTC', ('W','Grizzly'):'InHouse',
 ('W','Playland Wooden Coaster'):'InHouse', ('W','Boss'):'CCI', ('W','Wild One'):'Miller',
 ('W','T Express'):'Intamin', ('W','Megafobia'):'CCI', ('W','Hades 360'):'GravityGroup',
 ('W','Mine Blower'):'GravityGroup', ('W','Wooden Warrior'):'GravityGroup',
 ('W','Twister'):'InHouse',
}
# 同名不同厂: 两台木质 Blue Streak 按乐园区分
MFR_PARK = { ('W','Blue Streak','Conneaut Lake'):'InHouse', ('W','Blue Streak','Cedar Point'):'PTC' }

def parse_dur(s):
    if not s or ':' not in s: return None
    try:
        m, sec = s.split(':'); return int(m)*60 + int(sec)
    except ValueError: return None

def norm(s):
    s = unicodedata.normalize('NFKD', s)
    s = ''.join(c for c in s if not unicodedata.combining(c))
    return ''.join(c for c in s.lower() if c.isalnum())

rows = []
with open(f'{BASE}/rcdb_strict.csv') as f:
    for r in csv.DictReader(f):
        def num(k):
            try: return float(r[k])
            except (ValueError, TypeError): return None
        rows.append({'name': r['CoasterName'].split(' / ')[0], 'park': r['Park'], 'city': r['City'],
                     'status': r['Status'], 'h': num('Height'), 's': num('Speed'), 'len': num('Length'),
                     'inv': num('Inversions'), 'dur': parse_dur(r.get('Duration', ''))})
for r in rows:
    r['nn'] = norm(r['name']); r['np'] = norm(r['park'] + r['city'])

def find(name, park):
    n, p = norm(name), norm(park)
    exact = [r for r in rows if r['nn'] == n]
    if len(exact) == 1: return exact[0]
    if len(exact) > 1:
        pk = [r for r in exact if p[:8] in r['np'] or r['np'][:8] in p]
        return pk[0] if pk else None
    sub = [r for r in rows if n in r['nn'] or (r['nn'] in n and len(r['nn']) >= 5)]
    if sub:
        pk = [r for r in sub if p[:6] in r['np']]
        if len(pk) == 1: return pk[0]
        if len(sub) == 1: return sub[0]
    best, score = None, 0
    for r in rows:
        sc = difflib.SequenceMatcher(None, n, r['nn']).ratio()
        if sc > score: best, score = r, sc
    if score >= 0.85 and (norm(park)[:5] in best['np'] or score >= 0.93):
        return best
    return None

usable = [r for r in rows if r['h'] and r['s'] and r['len']]
bothd = [(r['len']/r['s'], r['dur']) for r in usable if r['dur']]
n_ = len(bothd); mx = sum(x for x, _ in bothd)/n_; my = sum(y for _, y in bothd)/n_
tb = sum((x-mx)*(y-my) for x, y in bothd)/sum((x-mx)**2 for x, _ in bothd)
ta = my - tb*mx
for r in usable:
    r['dur_i'] = r['dur'] if r['dur'] else max(20.0, ta + tb*r['len']/r['s'])
    r['inv_i'] = r['inv'] if r['inv'] is not None else 0.0

FEATS = ['s', 'h', 'len', 'dur_i', 'inv_i']
stats = {}
for f_ in FEATS:
    v = [r[f_] for r in usable]
    m = sum(v)/len(v); sd = math.sqrt(sum((x-m)**2 for x in v)/len(v))
    stats[f_] = (m, sd)
def x_of(r): return [(r[f_]-stats[f_][0])/stats[f_][1] for f_ in FEATS]

def build(list_, tag):
    out = []
    for rank, name, park in list_:
        r = find(name, park)
        if r and r['h'] and r['s'] and r['len']:
            m = MFR_PARK.get((tag, name, park.split(' ')[0] + (' ' + park.split(' ')[1] if len(park.split(' ')) > 1 else '')))
            if m is None:
                m = MFR_PARK.get((tag, name, park)) or MFR.get((tag, name))
            out.append({'rank': rank, 'row': r, 'gname': name, 'mfr': m})
    return out

steel = build(STEEL, 'S'); wood = build(WOOD, 'W')
# Blue Streak 特判
for e in wood:
    if e['gname'] == 'Blue Streak':
        e['mfr'] = 'InHouse' if 'Conneaut' in e['row']['park'] + e['row']['city'] else 'PTC'
unk = [e['gname'] for e in steel + wood if not e['mfr']]
if unk: print('未标注厂商:', unk)

def pairs_of(lst):
    ps = []
    for i in range(len(lst)):
        for j in range(len(lst)):
            if lst[i]['rank'] < lst[j]['rank']:
                ps.append((x_of(lst[i]['row']), x_of(lst[j]['row'])))
    return ps

def fit(pairs, lam=0.02, iters=900, lr=0.1):
    k = len(FEATS); w = [0.0]*k
    for _ in range(iters):
        g = [0.0]*k
        for xa, xb in pairs:
            d = sum(w[t]*(xa[t]-xb[t]) for t in range(k))
            pr = 1/(1+math.exp(-max(-30, min(30, d))))
            for t in range(k):
                g[t] += (1-pr)*(xa[t]-xb[t])
        for t in range(k):
            w[t] += lr*(g[t]/len(pairs) - lam*w[t])
    return w

w = fit(pairs_of(steel) + pairs_of(wood))
for lst in (steel, wood):
    scored = sorted(lst, key=lambda e: -sum(w[t]*x_of(e['row'])[t] for t in range(len(w))))
    pr = {id(e['row']): i+1 for i, e in enumerate(scored)}
    for e in lst:
        e['over'] = pr[id(e['row'])] - e['rank']
allc = [e for e in steel + wood if e['mfr']]

groups = {}
for e in allc:
    groups.setdefault(e['mfr'], []).append(e['over'])
multi = {k: v for k, v in groups.items() if len(v) >= 3}
vals = [x for v in multi.values() for x in v]
gm = sum(vals)/len(vals)
sst = sum((x-gm)**2 for x in vals)
def r2_of(assign):
    return sum(len(v)*(sum(v)/len(v)-gm)**2 for v in assign)/sst
r2 = r2_of(list(multi.values()))
sizes = [len(v) for v in multi.values()]
cnt = 0; pool = vals[:]
B = 5000
for _ in range(B):
    random.shuffle(pool)
    idx = 0; assign = []
    for s_ in sizes:
        assign.append(pool[idx:idx+s_]); idx += s_
    if r2_of(assign) >= r2: cnt += 1
print(f'\n按制造商 (n>=3 的 {len(multi)} 组, 覆盖 {len(vals)}/{len(allc)} 台): '
      f'组间方差占比 R²={r2:.2f}, 置换检验 p={(cnt+1)/(B+1):.4f}')
for k, v in sorted(multi.items(), key=lambda kv: -sum(kv[1])/len(kv[1])):
    print(f'  {k:13s} n={len(v):2d}  组均值 {sum(v)/len(v):+6.1f}')
small = {k: v for k, v in groups.items() if len(v) < 3}
if small:
    print('  (n<3 未入检验:', ', '.join(f'{k}({len(v)}台,{sum(v)/len(v):+.0f})' for k, v in small.items()), ')')