Error-budgeted local log-signature tokens / experiment.py

Failed on benchmark

Raw ⬇ ZIP
  1import math, json, time
  2import numpy as np
  3from sklearn.linear_model import LogisticRegression
  4from sklearn.model_selection import train_test_split
  5from sklearn.metrics import accuracy_score
  6
  7
  8def mobius(n):
  9    x, distinct = n, 0
 10    p = 2
 11    while p * p <= x:
 12        if x % p == 0:
 13            x //= p; distinct += 1
 14            if x % p == 0: return 0
 15            while x % p == 0: x //= p
 16        p += 1
 17    if x > 1: distinct += 1
 18    return -1 if distinct % 2 else 1
 19
 20
 21def lie_dim(k, d):
 22    return sum(mobius(q) * d ** (k // q) for q in range(1, k + 1) if k % q == 0) // k
 23
 24
 25def tail(z, N):
 26    term, partial = 1.0, 1.0
 27    for k in range(1, N + 1):
 28        term *= z / k
 29        partial += term
 30    return math.exp(z) - partial
 31
 32
 33def proxy(ells, r, N):
 34    return float(sum(tail(r * e, N) for e in ells))
 35
 36
 37def cost(m, N, d):
 38    return m * sum(lie_dim(k, d) for k in range(1, N + 1))
 39
 40
 41def select(ells, d, r, eps, Nmax=5):
 42    m, = [len(ells)]
 43    feasible = [(cost(m, N, d), N, proxy(ells, r, N))
 44                for N in range(1, Nmax + 1) if proxy(ells, r, N) <= eps]
 45    return min(feasible) if feasible else None
 46
 47
 48def variation_quantile_segments(x, m):
 49    ds = np.linalg.norm(np.diff(x, axis=0), axis=1)
 50    v = np.r_[0.0, np.cumsum(ds)]
 51    cuts = [0] + [int(np.searchsorted(v, q * v[-1] / m, side="left")) for q in range(1, m)] + [len(x)-1]
 52    cuts = np.maximum.accumulate(np.asarray(cuts)); cuts[-1] = len(x)-1
 53    ells = np.array([v[cuts[j+1]] - v[cuts[j]] for j in range(m)])
 54    return cuts, ells
 55
 56
 57def local_features(x, cuts, degree=2):
 58    # First-level increments and antisymmetric degree-2 log-signature area.
 59    feats = []
 60    for a, b in zip(cuts[:-1], cuts[1:]):
 61        p = x[a:b+1]; dx = np.diff(p, axis=0)
 62        inc = p[-1] - p[0]
 63        f = list(inc) + [float(b-a) / max(1, len(x)-1)]
 64        if degree >= 2:
 65            area = np.zeros((x.shape[1], x.shape[1]))
 66            prefix = np.zeros(x.shape[1])
 67            for u in dx:
 68                area += 0.5 * (np.outer(prefix, u) - np.outer(u, prefix))
 69                prefix += u
 70            f += [area[0,1]] if x.shape[1] >= 2 else []
 71        feats.extend(f)
 72    return np.asarray(feats, dtype=np.float32)
 73
 74
 75def make_data(n=700, length=64, seed=7):
 76    rng = np.random.default_rng(seed); X=[]; y=[]
 77    t=np.linspace(0,1,length)
 78    for i in range(n):
 79        cls=i%2; freq=3 if cls==0 else 7
 80        phase=rng.uniform(0,2*np.pi); amp=1+rng.normal(0,.04)
 81        s=amp*np.sin(2*np.pi*freq*t+phase)+rng.normal(0,.10,length)
 82        # second channel makes signed area informative for local signatures
 83        q=amp*np.cos(2*np.pi*freq*t+phase)+rng.normal(0,.10,length)
 84        X.append(np.stack([s,q],1)); y.append(cls)
 85    return np.asarray(X),np.asarray(y)
 86
 87
 88def features_fixed(X, m, degree):
 89    out=[]
 90    for x in X:
 91        cuts=np.linspace(0,len(x)-1,m+1).astype(int)
 92        out.append(local_features(x,cuts,degree))
 93    return np.asarray(out)
 94
 95
 96def run_classifier(X,y,F):
 97    a,b=train_test_split(np.arange(len(y)),test_size=.3,random_state=19,stratify=y)
 98    clf=LogisticRegression(max_iter=500,C=1.0).fit(F[a],y[a])
 99    return accuracy_score(y[b],clf.predict(F[b]))
100
101
102def main():
103    np.random.seed(7)
104    # Core math verification: exact equal-variation predictions.
105    print('LIE_DIMS', [lie_dim(k,2) for k in range(1,6)])
106    print('TAIL_SWEEP')
107    rows=[]
108    for z in [0.25,0.5,1.0,2.0,3.0]:
109        vals=[tail(z,N) for N in range(1,6)]
110        rows.append((z,vals))
111        print(z, ['%.8g'%v for v in vals])
112    # Prediction 1: for fixed total rL, equal m gives m R_N(rL/m), decreasing with m.
113    print('PREDICTION_EQUAL_M')
114    for ztot in [1.,3.,6.]:
115        vals=[proxy([ztot/m]*m,1.,2) for m in [1,2,4,8,16]]
116        print(ztot, ['%.8g'%v for v in vals])
117    # Prediction 2: raising N strictly decreases tail, with factorial-like ratios.
118    print('PREDICTION_DEGREE')
119    for z in [.5,1.,2.]:
120        vals=[tail(z,N) for N in range(1,6)]
121        print(z, ['%.8g'%v for v in vals], 'ratios', ['%.5g'%(vals[i+1]/vals[i]) for i in range(4)])
122    # Prediction 3: adaptive choice transitions toward larger m for oscillatory/high variation paths.
123    print('SELECTION_SWEEP')
124    for L in [1.,2.,4.,8.]:
125        eps=.05; picks=[]
126        for m in [1,2,4,8,16,32]:
127            # equal-variation proxy, retaining lowest cost across degrees
128            cand=select([L/m]*m,2,1.,eps)
129            if cand: picks.append((m,cand[1],cand[0],cand[2]))
130        print('L',L,'feasible_lowest_per_m',picks)
131    X,y=make_data()
132    results={}
133    t=time.time(); results['raw_8patch']=run_classifier(X,y,features_fixed(X,8,1))
134    results['fixed_sig_4_degree2']=run_classifier(X,y,features_fixed(X,4,2))
135    # adaptive: choose degree using observed proxy, and use 8 local tokens; representation is compact.
136    adaptive=[]; choices=[]
137    for x in X:
138        cuts,ells=variation_quantile_segments(x,8)
139        c=select(ells,2,1.,eps=.08)
140        N=2 if c is None else c[1]
141        choices.append(N); adaptive.append(local_features(x,cuts,N))
142    results['adaptive_8_local']=run_classifier(X,y,np.asarray(adaptive))
143    results['adaptive_degree_counts']={str(k):int(choices.count(k)) for k in sorted(set(choices))}
144    results['elapsed_sec']=time.time()-t
145    print('RESULTS',json.dumps(results,sort_keys=True))
146
147if __name__=='__main__': main()