Removing points from a triangular array without losing triangles

Question

I have a combinatorics problem that I'd like to put on the OEIS—the problem is that I don't have enough terms. This code challenge is to help me compute more terms, and the winner will be the user with the submission containing the greatest number of terms.

---

The Problem

Suppose I give you a triangular array of light bulbs with side length \$n\$:

     o
    o o
   o o o
  o o o o
 o o o o o
o o o o o o
1 2  ...  n

I'm going to turn on three lightbulbs that form an "upright" equilateral triangle as in the following example:

     o
    o x
   o o o
  o o o o
 o x o o x
o o o o o o

Before I turn on the lights, your job is to remove as many lightbulbs as possible from the array—without losing the ability to deduce the triangle of bulbs that has been turned on. To be clear, if a lightbulb has been removed, it is not lit up when its position is turned on.

For example, if you removed the following bulbs (marked by .) you would only see the following two lights turn on (marked by x), which is enough uniquely deduce the third (unlit) position:

     .              .
    . o            . x
   . . o          . . o
  o o o .   =>   o o o .
 o o o o .      o x o o . <- the third unlit position
o . . . o o    o . . . o o

Let a(n) be the maximum number of bulbs that can be removed without introducing any ambiguities.

---

Example

With a naive algorithm, I have checked values up to a triangle with side length 7, as seen below:

                                                                      .
                                                      .              . o
                                        .            . o            o . o
                           .           . .          . . o          . o o .
              .           . .         . o o        o o o .        o o . o .
 .           . .         . o o       o o . o      o o o o .      o . o . o o
. .         . o o       . o o o     o . . o o    o . . . o o    o . o . o o o

a(2) = 3    a(3) = 4    a(4) = 5    a(5) = 7     a(6) = 9       a(7) = 11

---

Scoring

The submission that computes the sequence [a(2), a(3), ..., a(n)] for the largest n wins. If two submissions have identical sequences, then the one that was posted earlier wins.

Although not necessary for the submission, it would be instructive to me if you post a construction of the resulting triangluar arrays, as in the example above.

Answer 1

Python 3, n=8

import itertools
from ortools.sat.python import cp_model


def solve(n):
    model = cp_model.CpModel()
    solver = cp_model.CpSolver()
    cells = {
        (y, x): model.NewBoolVar(str((y, x)))
        for y in range(n) for x in range(y + 1)}
    triangles = [
            {cells[v] for v in ((y1, x1), (y2, x1), (y2, x1 + y2 - y1))}
            for (y1, x1) in cells.keys() for y2 in range(y1 + 1, n)]
    for t1, t2 in itertools.combinations(triangles, 2):
        model.AddBoolOr(t1.symmetric_difference(t2))
    model.Minimize(sum(cells.values()))
    solver.Solve(model)
    return len(cells) - round(solver.ObjectiveValue())


for n in itertools.count(2):
    print('a(%d) = %d' % (n, solve(n)))

Uses Google OR-Tools' CP-SAT solver.

After running for ~30 seconds, it outputs the following:

a(2) = 3
a(3) = 4
a(4) = 5
a(5) = 7
a(6) = 9
a(7) = 11
a(8) = 13

I didn't event try to wait for n=9 as it would probably take hours (the number of constraints grows like \$n^6\$). After less than 30 minutes of computation I found out that a(9)=15. I'm leaving my score as n=8 because at the moment the time constraints are unclear, but half an hour is probably too long.

How it works

Take two distinct equilateral triangles \$T_1\$ and \$T_2\$. To avoid ambiguity, there should be at least one bulb on a vertex belonging to exactly one of \$T_1\$ and \$T_2\$.

Thus the question may be rephrased as a SAT problem, with one constraint for every pair of triangles.

PS: I would very much like to include an example for n=8, but I'm having issues with the SAT solver which apparently wants to keep the solutions all for itself.

Answer 2

Getting the solutions from @Delfad0r's program

I extended @Delfad0r's program to output solutions. It also gives intermediate results, so you get output like this:

Solving n = 8:
a(8) >= 9
a(8) >= 10
a(8) >= 11
a(8) >= 12
a(8) >= 13
       o
      . o
     . o o
    . o o .
   o o . o o
  o o o o . .
 o . . o o o .
o . . o . o o o
a(8) = 13

Solving n = 9:
a(9) >= 10
a(9) >= 13
a(9) >= 14
a(9) >= 15
        o
       o o
      o . .
     o . o o
    . o . o o
   . o o o o o
  o o o . o . .
 o o o . . . o o
. o o o o o o . .
a(9) = 15

This computation took several hours.

If you get impatient and press Ctrl-C after some possibly non-optimal solution was found, the program will show that solution. So it doesn't take long to get this:

                   .
                  o o
                 . o o
                . o o o
               o o . o o
              o . o o o .
             o . o . o o o
            . o o o o o . o
           o . . o o o o o o
          o o o o o o o o o .
         o o . o o o o . o o o
        o o o o o o . o . o o o
       o . o o . o o o o o o o o
      o o o . o o o o o . o o o o
     o o o . o o o o o o o o . . o
    o o o o o o o o o o o . o . . o
   o o o o . o o o o . o o o o o . o
  o o o o o o o o . o o . . o o o o .
 o o o o . o o . o . o o o o o o . o o
o o . o o . o o o o . o o o . o o o o o
a(20) >= 42

Here is the extended program:

import itertools
from ortools.sat.python import cp_model

class ReportPrinter(cp_model.CpSolverSolutionCallback):

    def __init__(self, n, total):
        cp_model.CpSolverSolutionCallback.__init__(self)
        self.__n = n
        self.__total = total

    def on_solution_callback(self):
        print('a(%d) >= %d' %
              (self.__n, self.__total-self.ObjectiveValue()) )

def solve(n):
    model = cp_model.CpModel()
    solver = cp_model.CpSolver()
    cells = {
        (y, x): model.NewBoolVar(str((y, x)))
        for y in range(n) for x in range(y + 1)}
    triangles = [
            {cells[v] for v in ((y1, x1), (y2, x1), (y2, x1 + y2 - y1))}
            for (y1, x1) in cells.keys() for y2 in range(y1 + 1, n)]
    for t1, t2 in itertools.combinations(triangles, 2):
        model.AddBoolOr(t1.symmetric_difference(t2))
    model.Minimize(sum(cells.values()))
    print('Solving n = %d:' % n)
    status = solver.SolveWithSolutionCallback(model, ReportPrinter(n,len(cells)))
    if status == cp_model.OPTIMAL:
        rel = '='
    elif status == cp_model.FEASIBLE:
        rel = '>='
    else:
        print('No result for a(%d)\n' % n)
        return
    for y in range(n):
        print(' '*(n-y-1), end='')
        for x in range(y+1):
            print('.o'[solver.Value(cells[(y,x)])],end=' ')
        print()
    print('a(%d) %s %d' % (n, rel, len(cells) - solver.ObjectiveValue()))
    print()

for n in itertools.count(2):
    solve(n)

Answer 3

Python 3

Based strongly on Delfad0r's answer, mostly follows the same logic progression by checking triangle pairs and validating the configuration if it contains no triangle pairs that fail this validation. Since I didn't use any libraries besides itertools and copy, I have full control over saving the examples encountered throughout the program.

examples = dict() # stores examples by key pair of n to a tuple with the triangle and number of lights turned off

for n in range(3, 8):
    tri = [] # model of the triangle, to be filled with booleans representing lights
    tri_points = [] # list of tuples representing points of the triangle
    for i in range(n):
        tri.append([True]*(i + 1))
        for j in range(i+1):
            tri_points.append((i, j))

    t_set = [] # list of all possible triangles from tri, represented by lists of points
    for i in range(n):
        for j in range(len(tri[i])):
            for k in range(1, n - i):
                t_set.append([(i, j), (i + k, j), (i + k, j + k)])

    from itertools import combinations
    import copy

    # validates whether or not a triangle of n lights can have i lights turned off, and saves an example to examples if validated
    def tri_validate(x):
        candidate_list = list(combinations(tri_points, x))
        tri_pairs = list(combinations(t_set, 2))
        for candidate in candidate_list:
            temp_tri = copy.deepcopy(tri)
            valid = False
            for point in candidate:
                (row, col) = point
                temp_tri[row][col] = False
            for pair in tri_pairs:
                valid = False
                (tri1, tri2) = pair
                for point in tri1:
                    if not valid:
                        if point not in tri2:
                            (row, col) = point
                            if temp_tri[row][col]:
                                valid = True
                for point in tri2:
                    if not valid:
                        if point not in tri1:
                            (row, col) = point
                            if temp_tri[row][col]:
                                valid = True
                if not valid:
                    break
            if valid:
                examples[n] = (temp_tri, x)
                return True
        return False

    # iterates up to the point that validation fails, then moves on to the next n
    for i in range(len(tri_points)):
        if tri_validate(i + 1):
            continue
        break

Problem is, it's not very efficient. It runs really fast up to n=5, but starts to slow down significantly past that point. At n=6, it takes around a minute to run, and it's much slower at n=7. I imagine there are a lot of efficiency fixes that can be made with this program, but it's a quickly made rough draft of a good solution with a lot more flexibility to check out the inner workings of this method. I'll be incrementally working on this over time.