#!/usr/bin/env python3
"""
Exact minimum-size covering designs invariant under the cyclic group Z_v
acting regularly on {0..v-1}, for a given cell (v,k,t).

Method (Kramer-Mesner style): enumerate Z_v-orbits of k-subsets; a
Z_v-invariant design is a selection of whole orbits. Solve the ILP
  min sum_j o_j * y_j      s.t. for every t-subset T: sum_j c_{j,T} y_j >= 1,
where o_j = |orbit j| and c_{j,T} in {0,1} marks whether the orbit covers T.
An OPTIMAL solution is the exact minimum over all cyclically invariant
coverings; feasibility/optimal status comes from the solver.

Usage: python cyclic_exact.py v k t seconds workers
"""
import sys, json
from itertools import combinations
from ortools.sat.python import cp_model

def schoenheim(v,k,t):
    val=1
    for i in range(t-1,-1,-1):
        val=-((-(v-i)*val)//(k-i))
    return val

def orbit_of(block, v):
    rots=set()
    cur=tuple(sorted(block))
    for _ in range(v):
        rots.add(cur)
        cur=tuple(sorted((x+1)%v for x in cur))
    return frozenset(rots)

def main():
    v,k,t=int(sys.argv[1]),int(sys.argv[2]),int(sys.argv[3])
    secs=float(sys.argv[4]) if len(sys.argv)>4 else 300.0
    workers=int(sys.argv[5]) if len(sys.argv)>5 else 6
    L=schoenheim(v,k,t)

    blocks=list(combinations(range(v),k))
    seen={}
    orbits=[]  # each: list of blocks (as tuples), representative first
    for b in blocks:
        o=orbit_of(b,v)
        key=min(o)
        if key not in seen:
            seen[key]=len(orbits)
            orbits.append(o)
    print(f"C({v},{k},{t}): {len(blocks)} blocks -> {len(orbits)} cyclic orbits; Schoenheim L={L}",flush=True)

    tid={}; tsubs=[]
    def rk(s):
        key=tuple(sorted(s))
        if key not in tid:
            tid[key]=len(tsubs); tsubs.append(key)
        return tid[key]
    orb_lists=[]
    orb_sizes=[]
    covers=[]
    for o in orbits:
        st=set()
        for rep in o:
            s=set(rep)
            for c in combinations(rep,t):
                st.add(rk(c))
        orb_sizes.append(len(o)); covers.append(st)
    nt=len(tsubs)
    cover_by=[[] for _ in range(nt)]
    for oi,st in enumerate(covers):
        for j in st: cover_by[j].append(oi)

    m=cp_model.CpModel()
    y=[m.NewBoolVar(f"o{i}") for i in range(len(orbits))]
    for j in range(nt):
        m.Add(sum(y[i] for i in cover_by[j])>=1)
    m.Minimize(sum(orb_sizes[i]*y[i] for i in range(len(orbits))))
    solver=cp_model.CpSolver()
    solver.parameters.max_time_in_seconds=secs
    solver.parameters.num_search_workers=workers
    status=solver.Solve(m)
    name=solver.StatusName(status)
    print("status:",name,flush=True)
    print("cyclic bound:",solver.ObjectiveValue(),"proven optimal:" , status==cp_model.OPTIMAL,"| best-bound:",solver.BestObjectiveBound(),flush=True)
    result={"cell":[v,k,t],"status":name,"optimal":status==cp_model.OPTIMAL,
      "objective":solver.ObjectiveValue(),"best_bound":solver.BestObjectiveBound(),
      "wall_time":solver.WallTime(),"branches":solver.NumBranches(),"conflicts":solver.NumConflicts(),
      "block_orbits":len(orbits),"t_subsets":nt}
    if status in (cp_model.OPTIMAL,cp_model.FEASIBLE):
        chosen=[i for i in range(len(orbits)) if solver.Value(y[i])]
        result["chosen_orbit_indices"]=chosen
        result["base_blocks"]=[list(min(orbits[i])) for i in chosen]
        print("base blocks:",result["base_blocks"],flush=True)
        blocks_out=[]
        for i in chosen:
            for rep in sorted(orbits[i]):
                blocks_out.append(tuple(sorted(rep)))
        blocks_out.sort()
        fn=f"cyclic_{v}_{k}_{t}_b{len(blocks_out)}.txt"
        with open(fn,"w") as f:
            f.write(f"{v} {k} {t}\n")
            for blk in blocks_out:
                f.write(" ".join(map(str,blk))+"\n")
        print("saved",fn,"with",len(blocks_out),"blocks from",len(chosen),"orbits",flush=True)
    result_fn=f"cyclic_exact_{v}_{k}_{t}.json"
    json.dump(result,open(result_fn,"w"),indent=1)
    print("solver stats:\n"+solver.ResponseStats(),flush=True)
    print("result metadata saved",result_fn,flush=True)

if __name__=="__main__":
    main()