Subsets of the square grid without isoscales right triangles

isoscales.py

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt
import matplotlib.patches as patches


def enumerate_isoscales(n):
    triangles = []
    for x1 in range(n):
        for y1 in range(n):
            p1 = np.array([x1, y1])
            for x2 in range(n):
                for y2 in range(n):
                    if x1 == x2 and y1 == y2:
                        continue
                    p2 = np.array([x2, y2])
                    delta = p2 - p1
                    d_o = delta[::-1]
                    d_o[1] = - d_o[1]
                    # print(p1.shape, p2.shape, d_o.shape) ; exit()
                    for p3 in (p2 + d_o, p2 - d_o):
                        if 0 <= p3[0] < n and 0 <= p3[1] < n:
                            triangles.append(np.array([p1, p2, p3]))
    triangles = np.array(triangles)
    return triangles


n = 8

triangles = enumerate_isoscales(n)
print(triangles.shape)


x = cp.Variable((n, n), boolean=True)


objective = cp.Maximize(cp.sum(x))

constraints = []
for triangle in triangles:
    constraints.append(cp.sum([x[i, j] for (i, j) in triangle]) <= 2)


problem = cp.Problem(objective, constraints)
problem.solve(verbose=True)
x = x.value.round().astype(int)
print(x)
print(n, x.sum(), x.sum() / n / n)

exit()


for triangle in triangles[::99][:10]:
    print(triangle)
    fig, ax = plt.subplots()
    polygon = patches.Polygon(triangle, closed=True, edgecolor='r', fill=None)
    ax.add_patch(polygon)

    # Set the aspect of the plot to be equal
    ax.set_aspect('equal')

    # Set limits based on the data
    ax.set_xlim(0, n - 1)
    ax.set_ylim(0, n - 1)

    plt.show()