import itertools, random, json
import numpy as np

# Exact small-strand positive Garside arithmetic. Simple braids are indexed by
# permutations; all operations below are exact for n=4 (24 permutations).
def compose(p, q):
    return tuple(p[q[i]] for i in range(len(p)))

def inverse(p):
    z = [0] * len(p)
    for i, v in enumerate(p): z[v] = i
    return tuple(z)

def swap_perm(n, i):
    p = list(range(n)); p[i], p[i+1] = p[i+1], p[i]
    return tuple(p)

def invcount(p):
    return sum(p[i] > p[j] for i in range(len(p)) for j in range(i+1, len(p)))

def delta_word(n):
    return [i for row in range(n-1, 0, -1) for i in range(1, row+1)]

def word_perm(word, n):
    p = tuple(range(n))
    for x in word: p = compose(p, swap_perm(n, x-1))
    return p

def all_simples(n):
    return list(itertools.permutations(range(n)))

def left_divides(a, b):
    # a divides b iff a^{-1}b is positive and lengths add; for simples this
    # is the left weak-order test.
    return invcount(a) + invcount(compose(inverse(a), b)) == invcount(b)

def left_gcd(a, b, simples):
    common = [c for c in simples if left_divides(c, a) and left_divides(c, b)]
    return max(common, key=invcount)

def normal_form_positive(word, n):
    """Canonical left-greedy factors for a positive generator word.

    Append generators as singleton factors and repeatedly left-weight adjacent
    factors. For simples x,y, transfer gcd(partial(x), y) from y to x.
    """
    simples = all_simples(n); identity = tuple(range(n)); delta = tuple(reversed(range(n)))
    fs = []
    for g in word:
        fs.append(swap_perm(n, g-1))
        i = len(fs)-2
        while i >= 0:
            x, y = fs[i], fs[i+1]
            complement = compose(inverse(x), delta)
            a = left_gcd(complement, y, simples)
            if a == identity: break
            fs[i] = compose(x, a)
            fs[i+1] = compose(inverse(a), y)
            if fs[i+1] == identity: fs.pop(i+1)
            i -= 1
    return fs

def factor_id(p):
    b = len(p)+1
    return sum(v*b**i for i,v in enumerate(p))

def strip_left_deltas(word, n):
    """Remove literal removable Delta powers, then normalize the suffix.

    This is exact for the augmentation used below (literal Delta^k inserted
    on the left). A full arbitrary-presentation Garside normalizer needs the
    usual complete cycling/transport implementation and is not claimed here.
    """
    d = delta_word(n); p = 0; pos = 0
    while word[pos:pos + len(d)] == d:
        p += 1; pos += len(d)
    return p, normal_form_positive(word[pos:], n)

def check_identities(n=4):
    d = delta_word(n)
    relation = [word_perm(d+[i],n) == word_perm([n-i]+d,n) for i in range(1,n)]
    # Key local rewrite checks, including 121=Delta for B3.
    local = normal_form_positive([1,2,1],3)
    return {'delta_word_length':len(d),
            'delta_permutation_is_longest':word_perm(d,n)==tuple(reversed(range(n))),
            'delta_sigma_relation_all':all(relation),
            'delta_normal_form_factor_count':len(normal_form_positive(d,n)),
            'b3_121_normalizes_to_delta':local==[tuple(reversed(range(3)))]}

def make_data(seed=7,n=4,ntrain=700,ntest=300):
    rng=random.Random(seed); d=delta_word(n); train=[]; test=[]
    for _ in range(ntrain):
        base=[rng.randint(1,n-1) for _ in range(rng.randint(8,22))]
        y=sum(x==1 for x in base)%2; train.append((base,y,0,base))
    for _ in range(ntest):
        base=[rng.randint(1,n-1) for _ in range(rng.randint(8,22))]
        y=sum(x==1 for x in base)%2; k=rng.randint(1,3)
        test.append((d*k+base,y,k,base))
    return train,test

def featurize(examples,n,mode):
    out=[]
    for word,_,_,_ in examples:
        if mode=='raw': z,size=word,n
        else:
            _,fs=strip_left_deltas(word,n); z=[factor_id(p) for p in fs]; size=(n+1)**n
        v=np.bincount(z,minlength=size).astype(float); v/=max(1,len(z)); out.append(v)
    return np.asarray(out)

def nearest_centroid(x,y,tx):
    cents=np.vstack([x[y==c].mean(0) for c in [0,1]])
    return np.argmin(((tx[:,None,:]-cents[None,:,:])**2).sum(2),axis=1)

def run(seed=7):
    n=4; train,test=make_data(seed,n); metrics={}
    for mode in ['raw','factor']:
        x=featurize(train,n,mode); tx=featurize(test,n,mode)
        y=np.array([e[1] for e in train]); ty=np.array([e[1] for e in test])
        metrics[mode+'_accuracy']=float(np.mean(nearest_centroid(x,y,tx)==ty))
    ratios=[]; exact=[]; raw_factor=[]
    for word,_,k,base in test:
        p,fs=strip_left_deltas(word,n); ratios.append(len(word)/max(1,len(fs)))
        exact.append([factor_id(x) for x in fs]==[factor_id(x) for x in normal_form_positive(base,n)])
        raw_factor.append((len(word),len(fs)))
    metrics.update({'mean_presented_generator_length':float(np.mean([a for a,b in raw_factor])),
      'mean_delta_free_factor_length':float(np.mean([b for a,b in raw_factor])),
      'mean_generator_to_factor_ratio':float(np.mean(ratios)),
      'delta_removal_factor_invariance':float(np.mean(exact)),
      'expected_delta_added_generators':len(delta_word(n))})
    return {'identity_check':check_identities(n),'metrics':metrics}

if __name__=='__main__': print(json.dumps(run(),sort_keys=True,indent=2))