import facile

circuit = [facile.variable(0, Nv) for i in range(Nvmax + 1)]

facile.constraint(circuit[0] == 1)

for j in range(1, Nvmax + 1):
    for i in range(j):
        facile.constraint((circuit[j] == 0) | (circuit[i] != 0))

va = facile.array(Va)
ov = facile.array(Ov)
dv = facile.array(Dv)

c = facile.array(circuit)
count_0 = sum([x == 0 for x in c])
csts = circuit[1] == 0
for city in set(Vh):
    origin = (va[ov[c[1]]] == city)
    destination = (va[dv[c[Nvmax - count_0]]] == city)
    # et hop, la contrainte!
    csts = csts | (origin & destination)

facile.constraint(csts)

td = facile.array(Td)
ta = facile.array(Ta)
dt = facile.array([item for sub in Dt for item in sub])

for k in range(1, Nvmax):
    same_city = va[dv[circuit[k]]] == va[ov[circuit[k+1]]]
    # Attention à la linéarisation des indices
    transfert = (td[circuit[k+1]] - ta[circuit[k]] >=
                 dt[dv[circuit[k]] * (Na + 1) + ov[circuit[k+1]]])
    facile.constraint((same_city & transfert) | (circuit[k+1] == 0))

duration = Dda + ta[c[Nvmax - count_0]] - td[circuit[1]] <= Dmax
facile.constraint(duration | (circuit[1] == 0))

circuits = facile.solve_all(circuit[1:])
print (circuits)