import numpy as np, itertools
from twonight import step
rng=np.random.default_rng(43)
pats=[p for l in (1,2,3) for p in itertools.product((0,1),repeat=l) if any(p)]
def feats(a):
    n=len(a); f=[]
    for p in pats:
        l=len(p); hits=np.ones(n-l+1,dtype=bool)
        for j in range(l): hits&=(a[j:n-l+1+j]==p[j])
        pos=np.arange(1,n-l+2)
        f+= [int(hits[pos%2==0].sum())%2, int(hits[pos%2==1].sum())%2]
    return f
def rank2(M):
    A=M.copy()%2; r=0
    for c in range(A.shape[1]):
        p=next((i for i in range(r,len(A)) if A[i,c]),None)
        if p is None: continue
        A[[r,p]]=A[[p,r]]
        for i in range(len(A)):
            if i!=r and A[i,c]: A[i]^=A[r]
        r+=1
    return r
def nullsp(M):
    A=M.copy()%2; n=A.shape[1]; r=0; piv=[]
    for c in range(n):
        p=next((i for i in range(r,len(A)) if A[i,c]),None)
        if p is None: continue
        A[[r,p]]=A[[p,r]]
        for i in range(len(A)):
            if i!=r and A[i,c]: A[i]^=A[r]
        piv.append(c); r+=1
    basis=[]
    for f in [c for c in range(n) if c not in piv]:
        v=np.zeros(n,dtype=np.uint8); v[f]=1
        for i,c in enumerate(piv):
            if A[i,f]: v[c]=1
        basis.append(v)
    return np.array(basis).T if basis else np.zeros((n,0),dtype=np.uint8)
K=6
for d in (1,0):
    X=[];Y=[];Lc=[]
    for _ in range(1500):
        Lw=int(rng.integers(8,26)); a=np.zeros(Lw+8,dtype=np.uint8); a[:Lw]=rng.integers(0,2,size=Lw)
        X.append(feats(a)); Y.append(feats(step(a,d))); Lc.append([int(a[i]) for i in range(K)]+[1])
    X=np.array(X,dtype=np.uint8); Y=np.array(Y,dtype=np.uint8); Lc=np.array(Lc,dtype=np.uint8)
    nf=X.shape[1]
    # remove combinations that are identically zero on the data (trivial)
    V=np.eye(nf,dtype=np.uint8)
    for it in range(30):
        m=V.shape[1]
        big=np.concatenate([(Y@V)%2,(X@V)%2,Lc],1)
        Nb=nullsp(big)
        C=Nb[:m,:]                      # coefficients c (in V basis) whose tomorrow-value is expressible
        if C.shape[1]==0: V=np.zeros((nf,0),dtype=np.uint8); break
        Vn=(V@C)%2
        r=rank2(Vn.T)
        if r==m: break
        # keep a basis of the column space
        A=Vn.T.copy(); keep=[]; R=np.zeros((0,nf),dtype=np.uint8)
        for col in A:
            T=np.vstack([R,col])
            if rank2(T)>len(R): R=T
        V=R.T
    triv=rank2(((X@V)%2).T) if V.shape[1] else 0
    print(f"{'LIT' if d else 'DARK'} night: closed subspace of global pattern-parities has dimension {V.shape[1]}, of which non-trivial (varies across crowds): {triv}")
names=[]
for p in pats:
    s="".join(map(str,p)); names+=[f"#{s}@even",f"#{s}@odd"]
# find the single varying closed combination: take V columns, pick one whose values vary
vals=(X@V)%2
for j in range(V.shape[1]):
    if len(set(vals[:,j]))>1:
        v=V[:,j]; break
print("closed parity  Pi =", " + ".join(names[i] for i in range(nf) if v[i]))
# its rule: Pi(next) as combination of Pi(today) and lanterns 1..K
for d in (1,0):
    rows=[];ys=[]
    for _ in range(3000):
        Lw=int(rng.integers(8,26)); a=np.zeros(Lw+8,dtype=np.uint8); a[:Lw]=rng.integers(0,2,size=Lw)
        pi=int(np.dot(feats(a),v)%2); pin=int(np.dot(feats(step(a,d)),v)%2)
        rows.append([pi]+[int(a[i]) for i in range(K)]+[1]); ys.append(pin)
    A=np.array(rows); y=np.array(ys)
    import itertools as it
    best=None
    for mask in it.product((0,1),repeat=A.shape[1]):
        if np.all((A@np.array(mask))%2==y): best=mask; break
    nm=["Pi"]+[f"a{i+1}" for i in range(K)]+["1"]
    print(f"   {'LIT' if d else 'DARK'} night:  Pi(next) =", " + ".join(nm[i] for i in range(len(nm)) if best and best[i]) if best else "linear rule not found")