#!/usr/bin/env python
# coding: utf-8

# In[39]:


n = 10
r = 4
meja = floor(sqrt((n * (n + 1) * (2*n + 1)) / (r * 6)))
p = MixedIntegerLinearProgram(maximization=True)

w = p.new_variable(binary=True)
y = p.new_variable(binary=True)
z = p.new_variable(binary=True)
p.set_objective(sum(z[l] for l in range(1, meja + 1))) ##ciljna funkcija
p.add_constraint(z[0] == 1) ## robni pogoj, kvadrat velikosti 0 lahko vedno pokrijemo
for i in range(1, n + 1):
    p.add_constraint(sum(w[i,u,v] for u in range(1, meja + 1) for v in range(1, meja + 1)) == 1) ##pogoj, da imamo največ en kvadrat dolžine i
for j in range(1, meja + 1):
    for k in range(1, meja + 1):
        p.add_constraint(r*y[j,k] <= sum(sum(sum
                           (w[i,u,v] for v in range(max(1, k - i + 1), k + 1)) ## pogoj, da je kvadratek (j,k) pokrit vsaj r-krat
                               for u in range(max(1, j - i + 1), j + 1))
                                  for i in range(1, n + 1))
                        )
for l in range(1, meja + 1):
    p.add_constraint(2*l*z[l] <= (z[l-1] + y[l,l] + sum((y[m,l] +y[l,m]) for m in range(1, l)))) ## pogoj, da je kvadrat lxl pokrit



# In[40]:


p.solve()  ##vrne velikost največjega kvadrata, ki ga dobimo s pokrivanjem.


# In[41]:


resw = p.get_values(w)
resy = p.get_values(y)
resz = p.get_values(z)
kvadratki = []
for (x,y) in resy.items():
    if y == 1:
        kvadratki.append(x)
KVADRATI = []
for (x,y) in resw.items():
    if y == 1:
        KVADRATI.append(x)


# In[42]:


KVADRATI  ## vrne izhodišča kvadratov v optimalni rešitvi