#!/usr/bin/env python3

"""A model of how the V9938 hands VRAM access slots to the CPU, and its fit.

Usage:

    cpu-arbitration-model.py <measurement-repo>            # both sets

    cpu-arbitration-model.py <measurement-repo> --2013     # only the old set

    cpu-arbitration-model.py <measurement-repo> --2026     # only the new set

The model (all times in VDP cycles):

  1. The CPU's request register is one deep.  It is loaded by a port access to

     #98 and stays occupied until the VRAM access has taken place -- not until

     the slot has been granted.  A port access that arrives while the register

     is occupied is lost: no VRAM access happens for it (and the VDP's VRAM

     pointer does not advance, so the CPU sees the previous byte again).

  2. LEAD cycles before every access slot, the VDP decides who gets it.  A

     pending CPU request wins; otherwise the command engine gets it if it has a

     request ready; otherwise the slot goes unused.  LEAD is counted in the

     VDP's *memory cycles*: the two memory cycles per line that are stretched by

     two cycles each in the horizontal blanking region do not count towards it,

     exactly as for the command engine's delays.

  3. A slot that has another slot only 6 cycles before it cannot carry a CPU

     access.  Such slots exist only in the sprites-off table (25 of the 88), and

     when the CPU's request wins one, the VDP spends that memory cycle on a dummy

     read of 0x1FFFF and the CPU's access happens in the next slot.  The decision

     point of these slots is two cycles earlier than the uniform rule.

What the two measurement sets contribute:

  * 2013 (Philips NMS8250): the CPU clock is exactly VDP/6, so the port accesses

    sit on an exact grid -- 40 of them 72 cycles apart, then 252 cycles for the

    loop overhead.  The phase of that grid relative to the VDP is unknown and is

    fitted per capture; a constant delay between the port access and the

    arbitration is indistinguishable from it.

  * 2026 (Philips NMS8280): the two clocks are asynchronous, but /CSR and /CSW

    were recorded, so the port accesses are known individually.  What is fitted

    per capture is then only the constant offset, and LEAD becomes measurable:

    26 cycles from the rising edge of /CSx.

Both sets are scored in both directions: accesses that happened but are not

predicted, and accesses that are predicted but did not happen.

"""

from __future__ import annotations

import argparse

import bisect

import math

import collections

import glob

import importlib.util

import os

import re

import statistics

import sys

TICKS = 1368

LEAD = 16                # the arbiter's lookahead, in memory cycles

DEAD = 2                 # cycles after an access during which a new port

                         # access is still lost

PAIR_EXTRA_LEAD = 2      # extra lead of a slot 6 cycles after another slot

FREE_BEFORE = [-DEAD]

PERIOD_2013 = 3060       # 39 * 72 + 252

IN_PERIOD = 72

BURST = 40

_here = os.path.dirname(os.path.abspath(__file__))

def _load(name, mod):

    spec = importlib.util.spec_from_file_location(mod,

                                                  os.path.join(_here, name))

    m = importlib.util.module_from_spec(spec)

    spec.loader.exec_module(m)

    return m

C = _load('2026-measurement-check.py', 'check')

G = _load('vcd-slot-grid.py', 'grid')

TABLES = C.TABLES

PAIR = {m: {s for s in TABLES[m] if (s - 6) in set(TABLES[m])} for m in TABLES}

# --------------------------------------------------------------------------

# the model

# --------------------------------------------------------------------------

def V(mode, t):

    """Memory-cycle time: real time minus the blanking stretch seen so far."""

    return C.V(mode, t)

def decision(mode, s):

    """Real time at which slot s is handed out.

    The lookahead is counted in the VDP's memory cycles, exactly as the command

    engine's delays are: the cycles the VDP pads the line with in the blanking

    region do not count towards it.

    """

    lead = LEAD + (PAIR_EXTRA_LEAD if (s % TICKS) in PAIR[mode] else 0)

    want = V(mode, s) - lead

    t = s - lead

    while V(mode, t) > want:

        t -= 1

    while V(mode, t + 1) <= want:

        t += 1

    return t

def slots(mode, lo, hi):

    tab = TABLES[mode]

    out = []

    line = lo // TICKS - 1

    while line * TICKS <= hi:

        out += [line * TICKS + s for s in tab if lo <= line * TICKS + s <= hi]

        line += 1

    return sorted(out)

def simulate(mode, requests, lo, hi, drop=True):

    """(cpu accesses, dummy reads, lost requests).

    'drop' selects what happens when a port access arrives while the register is

    still occupied: True loses the new one (what the measurements show), False

    overwrites and loses the old one (what openMSX assumes).

    """

    ss = slots(mode, lo - 200, hi + 200)

    nxt = {s: n for s, n in zip(ss, ss[1:])}

    # 0 = access completes, 1 = port access, 2 = slot is handed out.  A port

    # access that lands exactly on a decision point still wins the slot, i.e.

    # the lead is a shade under 26 rather than exactly 26; five accesses in the

    # 2013 set say so, and the 2026 measurement of it is 25.98 +- 0.11.

    ev = ([(s - FREE_BEFORE[0], 0, s) for s in ss]

          + [(r, 1, r) for r in requests]

          + [(decision(mode, s), 2, s) for s in ss])

    ev.sort()

    busy = False            # register holds a request that has not been served

    pending = False         # ... and it has not been given a slot yet

    at = None               # slot in which its access will happen

    cpu, dummy, lost = [], [], 0

    for t, kind, val in ev:

        if kind == 1:

            if busy:

                lost += 1

                if not drop:

                    pending, at = True, None

            else:

                busy, pending, at = True, True, None

        elif kind == 2:

            if pending and at is None:

                if (val % TICKS) in PAIR[mode]:

                    dummy.append(val)

                    at = nxt.get(val)

                else:

                    at = val

                pending = False

        elif at == val:

            cpu.append(val)

            busy, at = False, None

        elif at is None and not pending:

            pass

    return cpu, dummy, lost

# --------------------------------------------------------------------------

# the 2013 set: exact request grid, unknown phase

# --------------------------------------------------------------------------

CMDS = ('HMMM', 'HMMV', 'LMMM', 'LMMV', 'YMMM', 'HYMM', 'LINE')

def meta13(path):

    b = os.path.basename(path)[:-4]

    if not b.startswith('screen') or 'wrong' in b:

        return None

    scr = re.match(r'screen(\d+)', b)

    if scr is None or scr[1] not in ('5', '6', '7', '8'):

        return None

    mode = ('dispOff' if 'screenoff' in b else

            'sprOff' if 'nosprites' in b else

            'sprOn' if 'sprite' in b else None)

    cpu = ('rd' if 'cpuread' in b else 'wr' if 'cpuwrite' in b else None)

    if mode is None or cpu is None:

        return None

    return {'mode': mode, 'name': b,

            'cmd': next((c for c in CMDS if c in b), None)}

def parse13(path, mode):

    rel = C.RELABEL[mode]

    out = []

    with open(path) as f:

        for ln in f:

            if ':' not in ln:

                continue

            row, rest = ln.rstrip('\n').split(':', 1)

            t, col = int(row), 0

            while rest:

                cell = rest[2:13]

                if cell.strip():

                    line, c = divmod(TICKS * col + t - 1, TICKS)

                    out.append((line * TICKS + rel.get(c, c), cell[0:3],

                                int(cell[6:11], 16)))

                rest = rest[13:]

                col += 1

    return sorted(out)

def requests13(phase, lo, hi):

    out = []

    n = (lo - phase) // PERIOD_2013 - 1

    while phase + n * PERIOD_2013 <= hi:

        base = phase + n * PERIOD_2013

        out += [base + j * IN_PERIOD for j in range(BURST)

                if lo <= base + j * IN_PERIOD <= hi]

        n += 1

    return sorted(out)

def run2013(repo, drop=True):

    print('=== 2013, Philips NMS8250: CPU clock exactly VDP/6 ===')

    print(f'{"capture":42} {"mode":8} {"n":>4} {"missed":>7} {"spurious":>9} '

          f'{"lost":>5} {"phases":>7}')

    tot = collections.Counter()

    for path in sorted(glob.glob(os.path.join(repo, '8.final-analysis',

                                              'screen*.txt'))):

        m = meta13(path)

        if m is None:

            continue

        acc = parse13(path, m['mode'])

        obs = sorted(t for t, code, a in acc if code in ('R.c', 'W.c'))

        if len(obs) < 10:

            continue

        lo, hi = obs[0] - 100, obs[-1] + 100

        best = None

        for phase in range(PERIOD_2013):

            cpu, dummy, lost = simulate(m['mode'], requests13(phase, lo, hi),

                                        lo, hi, drop)

            p = {x for x in cpu if obs[0] <= x <= obs[-1]}

            miss, extra = len(set(obs) - p), len(p - set(obs))

            score = (miss + extra, )

            if best is None or score < best[0]:

                best = (score, miss, extra, lost, [phase])

            elif score == best[0]:

                best[4].append(phase)

        _, miss, extra, lost, phases = best

        tot['n'] += len(obs)

        tot['miss'] += miss

        tot['extra'] += extra

        print(f'{m["name"]:42} {m["mode"]:8} {len(obs):4} {miss:7} {extra:9} '

              f'{lost:5} {len(phases):7}')

    print(f'\n{tot["n"]} CPU accesses: {tot["miss"]} not predicted '

          f'({100.0 * tot["miss"] / tot["n"]:.2f}%), {tot["extra"]} predicted '

          f'but absent\n')

# --------------------------------------------------------------------------

# the 2026 set: /CSx recorded, asynchronous clock

# --------------------------------------------------------------------------

def timebase(ref, cyc, anchor):

    """Anchors from the refresh chain, which is exactly 128 cycles apart."""

    out = []

    for t in sorted(ref):

        c = (t - anchor) * cyc + G.REFRESH_CYCLE

        line, off = divmod(c, TICKS)

        k = round((off - G.REFRESH_CYCLE) / 128)

        if 0 <= k <= 7:

            exact = line * TICKS + G.REFRESH_CYCLE + 128 * k

            if abs(exact - c) <= 3:

                out.append((t, exact))

    return out

def decode26(vcd, rng):

    trace = G.parse_vcd(vcd)

    edges = G.ras_edges(trace)

    g = G.geometry(edges)

    if g is None:

        return None

    period, _ = g

    cyc = TICKS / period

    acc = G.accesses(trace)

    anchor, ref = G.refresh_chain_start(acc, cyc)

    if anchor is None:

        return None

    fixed = timebase(ref, cyc, anchor)

    if len(fixed) < 4:

        return None

    xs = [t for t, _ in fixed]

    ys = [c for _, c in fixed]

    def cycle(t):

        i = bisect.bisect_left(xs, t)

        if i == 0:

            return ys[0] + (t - xs[0]) * cyc

        if i >= len(xs):

            return ys[-1] + (t - xs[-1]) * cyc

        f = (t - xs[i - 1]) / (xs[i] - xs[i - 1])

        return ys[i - 1] + f * (ys[i] - ys[i - 1])

    req, prev = [], {'CSR': None, 'CSW': None}

    for t, s in trace:

        for sig in ('CSR', 'CSW'):

            if prev[sig] == 0 and s[sig] == 1:

                req.append(cycle(t))

            prev[sig] = s[sig]

    by = collections.defaultdict(list)

    for cas, ras, addr, rd in acc:

        by[ras].append(addr)

    cpu, dummy = [], []

    for ras, addrs in sorted(by.items()):

        if len(addrs) != 1 or ras in ref:

            continue

        if addrs[0] == 0x1FFFF:

            dummy.append(cycle(ras))

        elif rng[0] <= addrs[0] <= rng[1]:

            cpu.append(cycle(ras))

    return sorted(req), sorted(cpu), sorted(dummy)

def snap(mode, t):

    line, off = divmod(t, TICKS)

    return min(((line + k) * TICKS + s for k in (-1, 0, 1)

                for s in TABLES[mode]), key=lambda v: abs(v - t))

def fit_pace(req):

    """(cycles per Z80 T-state, phase, residual rms, T-state indices, pulses).

    The port accesses are 12 T-states apart inside the burst of 40 IN/OUT

    instructions and 42 from the last of one burst to the first of the next.

    Fitting one straight line through the whole capture recovers the sequence

    far more accurately than the individual /CSx edges do -- the residuals come

    out at ~0.15 cycles, well under the analyzer's 0.27-cycle sampling step.

    A capture occasionally holds a stray pulse, so outliers are rejected.

    """

    # two /CSx edges a few cycles apart cannot both be port accesses

    req = [r for i, r in enumerate(req) if i == 0 or r - req[i - 1] > 40]

    g = [b - a for a, b in zip(req, req[1:])]

    # The test loops, in Z80 T-states between port accesses and from the last

    # access of a burst to the first of the next:

    #   IN/OUT x40                  12, 42   ('rdCpu', 'wrCpu', 'rdwrCpu')

    #   IN ; NOP ; NOP              22, 52   ('rdCpu132', 'wrCpu132')

    #   IN x20, slow                37, 67   ('rdCpu222')

    inner = statistics.median(x for x in g if x < 300)

    step, tail = min(((12, 42), (22, 52), (37, 67)),

                     key=lambda st: abs(inner / st[0] - 5.96))

    s0 = inner / step

    x = [0]

    for gap in g:

        n = round(gap / (step * s0))

        x.append(x[-1] + (tail if abs(gap - tail * s0) < 20 * s0 / step

                          else step * max(1, n)))

    keep = list(range(len(req)))

    s, phi, rms = s0, 0.0, 0.0

    for _ in range(4):

        n = len(keep)

        mx = sum(x[i] for i in keep) / n

        my = sum(req[i] for i in keep) / n

        sxx = sum((x[i] - mx) ** 2 for i in keep)

        sxy = sum((x[i] - mx) * (req[i] - my) for i in keep)

        s = sxy / sxx

        phi = my - s * mx

        res = [req[i] - (phi + s * x[i]) for i in keep]

        rms = (sum(e * e for e in res) / n) ** 0.5

        nk = [i for i, e in zip(keep, res) if abs(e) <= max(0.4, 3 * rms)]

        if len(nk) == len(keep) or len(nk) < 20:

            break

        keep = nk

    return s, phi, rms, x, req

def run2026(repo, drop=True, quantum=1):

    print('=== 2026, Philips NMS8280: CPU clock asynchronous to the VDP ===')

    vcd_dir = os.path.join(repo, 'part2', '1.vcd')

    txt_dir = os.path.join(repo, 'part2', '5.slots')

    agg = collections.defaultdict(collections.Counter)

    offs = collections.defaultdict(list)

    rates, rmss = [], []

    for vcd in sorted(glob.glob(os.path.join(vcd_dir, 'scr5-*.vcd'))):

        name = os.path.basename(vcd)[:-4]

        if 'noCpu' in name:

            continue

        txt = os.path.join(txt_dir, name + '.txt')

        if not os.path.exists(txt) or C.meta(txt) is None:

            continue

        mode = C.meta(txt)['mode']

        v = [a for t, code, a in C.parse(txt) if code in C.CPU_CODES]

        if not v:

            continue

        d = decode26(vcd, (min(v), max(v)))

        if d is None:

            continue

        req, cpu, dummy = d

        obs = sorted({snap(mode, t) for t in cpu})

        if len(obs) < 10:

            continue

        obsd = {snap(mode, t) for t in dummy

                if (snap(mode, t) % TICKS) in PAIR[mode]}

        s, phi, rms, xs, _ = fit_pace(req)

        rates.append(s)

        rmss.append(rms)

        base = [phi + s * a for a in xs]

        lo, hi = obs[0] - 100, obs[-1] + 100

        # the constant delay between the port access and the arbiter, fitted

        # per capture: it also absorbs the residual alignment of the capture's

        # own time base

        best = None

        o = -8.0

        while o < 24.0:

            # the VDP can only see a request at one of its own clock edges

            rq = [(quantum * math.ceil((b + o) / quantum) if quantum

                   else b + o) for b in base if lo - 150 < b + o < hi]

            p, pd, lost = simulate(mode, rq, lo, hi, drop)

            ps = {x for x in p if obs[0] <= x <= obs[-1]}

            pds = {x for x in pd if obs[0] <= x <= obs[-1]}

            score = (len(set(obs) - ps) + len(ps - set(obs))

                     + len(obsd - pds) + len(pds - obsd))

            if best is None or score < best[0]:

                best = (score, o, ps, pds, lost, len(rq))

            o += 0.05

        _, off, ps, pds, lost, nreq = best

        pace = ('222' if 'Cpu222' in name else

                '132' if 'Cpu132' in name else

                'r/w' if 'rdwrCpu' in name else ' 72')

        k = agg[(pace, mode)]

        k['n'] += len(obs)

        k['miss'] += len(set(obs) - ps)

        k['extra'] += len(ps - set(obs))

        k['dn'] += len(obsd)

        k['dmiss'] += len(obsd - pds)

        k['dextra'] += len(pds - obsd)

        k['lost'] += lost

        k['req'] += nreq

        k['cap'] += 1

        offs[(pace, mode)].append(off)

    print(f'pace: {statistics.fmean(rates):.5f} +- {statistics.stdev(rates):.5f}'

          f' VDP cycles per Z80 T-state ({12 * statistics.fmean(rates):.4f} '

          f'cycles between port accesses; 6 and 72 on a single-clock machine), '

          f'fit residual {statistics.median(rmss):.3f} cycles')

    print(f'requests quantised to {quantum} VDP cycle(s)'