1import itertools, random, json
  2import numpy as np
  3
  4# Exact small-strand positive Garside arithmetic. Simple braids are indexed by
  5# permutations; all operations below are exact for n=4 (24 permutations).
  6def compose(p, q):
  7    return tuple(p[q[i]] for i in range(len(p)))
  8
  9def inverse(p):
 10    z = [0] * len(p)
 11    for i, v in enumerate(p): z[v] = i
 12    return tuple(z)
 13
 14def swap_perm(n, i):
 15    p = list(range(n)); p[i], p[i+1] = p[i+1], p[i]
 16    return tuple(p)
 17
 18def invcount(p):
 19    return sum(p[i] > p[j] for i in range(len(p)) for j in range(i+1, len(p)))
 20
 21def delta_word(n):
 22    return [i for row in range(n-1, 0, -1) for i in range(1, row+1)]
 23
 24def word_perm(word, n):
 25    p = tuple(range(n))
 26    for x in word: p = compose(p, swap_perm(n, x-1))
 27    return p
 28
 29def all_simples(n):
 30    return list(itertools.permutations(range(n)))
 31
 32def left_divides(a, b):
 33    # a divides b iff a^{-1}b is positive and lengths add; for simples this
 34    # is the left weak-order test.
 35    return invcount(a) + invcount(compose(inverse(a), b)) == invcount(b)
 36
 37def left_gcd(a, b, simples):
 38    common = [c for c in simples if left_divides(c, a) and left_divides(c, b)]
 39    return max(common, key=invcount)
 40
 41def normal_form_positive(word, n):
 42    """Canonical left-greedy factors for a positive generator word.
 43
 44    Append generators as singleton factors and repeatedly left-weight adjacent
 45    factors. For simples x,y, transfer gcd(partial(x), y) from y to x.
 46    """
 47    simples = all_simples(n); identity = tuple(range(n)); delta = tuple(reversed(range(n)))
 48    fs = []
 49    for g in word:
 50        fs.append(swap_perm(n, g-1))
 51        i = len(fs)-2
 52        while i >= 0:
 53            x, y = fs[i], fs[i+1]
 54            complement = compose(inverse(x), delta)
 55            a = left_gcd(complement, y, simples)
 56            if a == identity: break
 57            fs[i] = compose(x, a)
 58            fs[i+1] = compose(inverse(a), y)
 59            if fs[i+1] == identity: fs.pop(i+1)
 60            i -= 1
 61    return fs
 62
 63def factor_id(p):
 64    b = len(p)+1
 65    return sum(v*b**i for i,v in enumerate(p))
 66
 67def strip_left_deltas(word, n):
 68    """Remove literal removable Delta powers, then normalize the suffix.
 69
 70    This is exact for the augmentation used below (literal Delta^k inserted
 71    on the left). A full arbitrary-presentation Garside normalizer needs the
 72    usual complete cycling/transport implementation and is not claimed here.
 73    """
 74    d = delta_word(n); p = 0; pos = 0
 75    while word[pos:pos + len(d)] == d:
 76        p += 1; pos += len(d)
 77    return p, normal_form_positive(word[pos:], n)
 78
 79def check_identities(n=4):
 80    d = delta_word(n)
 81    relation = [word_perm(d+[i],n) == word_perm([n-i]+d,n) for i in range(1,n)]
 82    # Key local rewrite checks, including 121=Delta for B3.
 83    local = normal_form_positive([1,2,1],3)
 84    return {'delta_word_length':len(d),
 85            'delta_permutation_is_longest':word_perm(d,n)==tuple(reversed(range(n))),
 86            'delta_sigma_relation_all':all(relation),
 87            'delta_normal_form_factor_count':len(normal_form_positive(d,n)),
 88            'b3_121_normalizes_to_delta':local==[tuple(reversed(range(3)))]}
 89
 90def make_data(seed=7,n=4,ntrain=700,ntest=300):
 91    rng=random.Random(seed); d=delta_word(n); train=[]; test=[]
 92    for _ in range(ntrain):
 93        base=[rng.randint(1,n-1) for _ in range(rng.randint(8,22))]
 94        y=sum(x==1 for x in base)%2; train.append((base,y,0,base))
 95    for _ in range(ntest):
 96        base=[rng.randint(1,n-1) for _ in range(rng.randint(8,22))]
 97        y=sum(x==1 for x in base)%2; k=rng.randint(1,3)
 98        test.append((d*k+base,y,k,base))
 99    return train,test
100
101def featurize(examples,n,mode):
102    out=[]
103    for word,_,_,_ in examples:
104        if mode=='raw': z,size=word,n
105        else:
106            _,fs=strip_left_deltas(word,n); z=[factor_id(p) for p in fs]; size=(n+1)**n
107        v=np.bincount(z,minlength=size).astype(float); v/=max(1,len(z)); out.append(v)
108    return np.asarray(out)
109
110def nearest_centroid(x,y,tx):
111    cents=np.vstack([x[y==c].mean(0) for c in [0,1]])
112    return np.argmin(((tx[:,None,:]-cents[None,:,:])**2).sum(2),axis=1)
113
114def run(seed=7):
115    n=4; train,test=make_data(seed,n); metrics={}
116    for mode in ['raw','factor']:
117        x=featurize(train,n,mode); tx=featurize(test,n,mode)
118        y=np.array([e[1] for e in train]); ty=np.array([e[1] for e in test])
119        metrics[mode+'_accuracy']=float(np.mean(nearest_centroid(x,y,tx)==ty))
120    ratios=[]; exact=[]; raw_factor=[]
121    for word,_,k,base in test:
122        p,fs=strip_left_deltas(word,n); ratios.append(len(word)/max(1,len(fs)))
123        exact.append([factor_id(x) for x in fs]==[factor_id(x) for x in normal_form_positive(base,n)])
124        raw_factor.append((len(word),len(fs)))
125    metrics.update({'mean_presented_generator_length':float(np.mean([a for a,b in raw_factor])),
126      'mean_delta_free_factor_length':float(np.mean([b for a,b in raw_factor])),
127      'mean_generator_to_factor_ratio':float(np.mean(ratios)),
128      'delta_removal_factor_invariance':float(np.mean(exact)),
129      'expected_delta_added_generators':len(delta_word(n))})
130    return {'identity_check':check_identities(n),'metrics':metrics}
131
132if __name__=='__main__': print(json.dumps(run(),sort_keys=True,indent=2))