#!/usr/bin/env python3
"""Independent exact branch-and-bound for cyclic covering designs.

Reduces a Z_v-invariant C(v,k,t) covering to set cover on the orbits of
t-subsets.  For the three target cells gcd(v,k)=1, all block orbits have
size v, so minimizing blocks is exactly minimizing selected block orbits.

Decision search proves that q orbit-unions cannot cover the universe; a
witness with q+1 orbits plus the failed q search proves exactness.

Usage: python cyclic_bnb.py v k t incumbent_orbits [seconds]
"""
import sys,time,json
from itertools import combinations

def orbit(s,v):
    s=tuple(sorted(s)); out=set()
    cur=s
    for _ in range(v):
        out.add(cur)
        cur=tuple(sorted((x+1)%v for x in cur))
    return frozenset(out)

def orbit_partition(v,r):
    reps={}
    for s in combinations(range(v),r):
        o=orbit(s,v); key=min(o)
        if key not in reps: reps[key]=o
    return [(key,reps[key]) for key in sorted(reps)]

def main():
    v,k,t=int(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3])
    incumbent=int(sys.argv[4])
    limit=float(sys.argv[5]) if len(sys.argv)>5 else 600.0
    torbs=orbit_partition(v,t); korbs=orbit_partition(v,k)
    tmap={s:i for i,(_,o) in enumerate(torbs) for s in o}
    weights=[len(o) for _,o in korbs]
    if len(set(weights))!=1:
        raise SystemExit(f"nonuniform block-orbit weights {sorted(set(weights))}; this decision solver assumes uniform")
    w=weights[0]
    covers=[]; reps=[]
    for rep,o in korbs:
        mask=0
        for sub in combinations(rep,t): mask |= 1<<tmap[tuple(sorted(sub))]
        covers.append(mask); reps.append(rep)
    # Remove duplicate and dominated masks: with equal cost, a strict subset is never useful.
    uniq={m:i for i,m in enumerate(covers)}
    masks=list(uniq)
    masks.sort(key=int.bit_count,reverse=True)
    kept=[]
    for m in masks:
        if any((m|q)==q for q in kept): continue
        kept.append(m)
    mask_to_rep={covers[i]:reps[i] for i in range(len(covers))}
    covers=kept; reps=[mask_to_rep[m] for m in covers]
    U=len(torbs); ALL=(1<<U)-1
    byelem=[[] for _ in range(U)]
    for i,m in enumerate(covers):
        for e in range(U):
            if (m>>e)&1: byelem[e].append(i)
    # Pairwise-incompatibility packing bound: selected universe elements that no
    # single block-orbit can cover together each require distinct choices.
    cooccur=[0]*U
    for m in covers:
        z=m
        while z:
            bit=z & -z; e=bit.bit_length()-1
            cooccur[e] |= m
            z ^= bit
    pack_order=sorted(range(U),key=lambda e:cooccur[e].bit_count())
    def packing_lb(um):
        n=0
        while um:
            e=next(e for e in pack_order if (um>>e)&1)
            n+=1; um &= ~cooccur[e]
        return n

    deadline=time.time()+limit; nodes=0; memo={}; witness=None; timed=False

    def dfs(uncovered,slots,chosen):
        nonlocal nodes,witness,timed
        nodes+=1
        if nodes%10000==0 and time.time()>deadline:
            timed=True; return False
        if uncovered==0:
            witness=list(chosen); return True
        if slots==0: return False
        prev=memo.get(uncovered,-1)
        if prev>=slots: return False
        memo[uncovered]=slots
        maxgain=0
        um=uncovered
        for m in covers:
            g=(m&um).bit_count()
            if g>maxgain: maxgain=g
        if maxgain==0 or (um.bit_count()+maxgain-1)//maxgain>slots: return False
        if packing_lb(um)>slots: return False
        # MRV uncovered orbit, then greatest gain first.
        elems=[e for e in range(U) if (um>>e)&1]
        e=min(elems,key=lambda z:sum(1 for i in byelem[z] if covers[i]&um))
        cand=[i for i in byelem[e] if covers[i]&um]
        cand.sort(key=lambda i:(covers[i]&um).bit_count(),reverse=True)
        for i in cand:
            if dfs(um & ~covers[i],slots-1,chosen+[i]): return True
            if timed: return False
        return False

    q=incumbent-1
    start=time.time(); sat=dfs(ALL,q,[]); elapsed=time.time()-start
    result={"cell":[v,k,t],"t_orbits":len(torbs),"block_orbits_raw":len(korbs),
      "block_orbits_after_dominance":len(covers),"orbit_weight":w,
      "decision_orbits":q,"decision_blocks":q*w,"status":("TIMEOUT" if timed else ("SAT" if sat else "UNSAT")),"sat":(None if timed else sat),"timed_out":timed,
      "nodes":nodes,"elapsed_seconds":elapsed,"incumbent_orbits":incumbent,
      "incumbent_blocks":incumbent*w}
    if witness is not None: result["witness_reps"]=[list(reps[i]) for i in witness]
    print(json.dumps(result,indent=1),flush=True)
    fn=f"cyclic_bnb_{v}_{k}_{t}.json"
    json.dump(result,open(fn,"w"),indent=1)
    if timed: sys.exit(2)
    if sat:
        print(f"DECISION SAT: cyclic cover exists with <= {q} orbits; incumbent not exact",flush=True); sys.exit(1)
    print(f"DECISION UNSAT by exhaustive B&B: no cyclic cover with {q} orbits; with incumbent {incumbent}, exact value={incumbent*w}",flush=True)

if __name__=="__main__": main()