Stable reversed Tartaglia's triangles

Question

Inspired by Minimising the triangle

---

I am writing a fully tested program to solve the following problem:

A triangle needs a good foundation. Every row in the triangle is derived from the sum of the two values below it. However, there can be no repeated values, if a value shows up more than once the triangle crumbles. Find the base which minimises the value in the top of the triangle satisfying the condition of no duplicates.

Example:

      20
    8   12
 [3]   5   7
1   2  [3]  4

Here 3 occurs twice, so the triangle is invalid.

My main concern is that my code has time complexity \$O(n!)\$ in the function stable_triangles, so it is not scalable at all.

"""
Finds `Reversed Tartaglia's (or Pascals') triangles` where there are
no repeated values.
"""

import doctest
from itertools import permutations

def next_level(level):
    """
    Builds the layer above of the triangle by
    summing the near pairs of numbers.

    >>> next_level([1, 4, 5])
    [5, 9]
    """
    return [curr + level[index + 1]
               for index, curr in enumerate(level[:-1])]

def make_triangle(base):
    """
    Builds a full triangle by `next_level` until the
    top has length 1.

    >>> make_triangle([1,2,3])
    [[1, 2, 3], [3, 5], [8]]
    """
    triangle = [base]
    while len(triangle[-1]) > 1:
        triangle.append(next_level(triangle[-1]))
    return triangle

def all_unique(lst):
    """
    Are the items all unique?

    >>> all_unique([4, 6, 7])
    True
    >>> all_unique([4, 4, 6, 7])
    False
    """
    return len(lst) == len(set(lst))

def flatten1(lst):
    """
    Flattens one level, shallow.

    >>> flatten1([ [1,2], [3,4] ])
    [1, 2, 3, 4]
    """
    return [i for sublst in lst for i in sublst]

def is_stable(triangle):
    """
    A triangle is stable if no number is repeated,
    that is: all the numbers are unique.

    >>> is_stable([[1,2,3], [3,5], [8]])
    False
    >>> is_stable([[1,2], [2]])
    False
    >>> is_stable([[1,3], [4]])
    True
    """
    return all_unique(flatten1(triangle))

def stable_triangles(base_length):
    """
    Yields many stable triangles with the given base length.

    >>> flattend_sum = lambda tri: sum(flatten1(tri))
    >>> min(stable_triangles(3), key=flattend_sum)
    [[2, 1, 4], [3, 5], [8]]
    >>> min(stable_triangles(4), key=flattend_sum)
    [[2, 3, 1, 6], [5, 4, 7], [9, 11], [20]]
    >>> min(stable_triangles(5), key=flattend_sum)
    [[6, 4, 1, 2, 7], [10, 5, 3, 9], [15, 8, 12], [23, 20], [43]]
    >>> next(stable_triangles(6))
    [[1, 2, 4, 7, 5, 8], [3, 6, 11, 12, 13], [9, 17, 23, 25], [26, 40, 48], [66, 88], [154]]

    # >>> next(stable_triangles(7)) # This takes about a minute. Comment out for faster testing.
    # [[1, 2, 4, 7, 5, 8, 10], [3, 6, 11, 12, 13, 18], [9, 17, 23, 25, 31], [26, 40, 48, 56], [66, 88, 104], [154, 192], [346]]
    """
    for base in permutations(range(base_length * 2), base_length):
        triangle = make_triangle(list(base))
        if is_stable(triangle):
            yield triangle

if __name__ == "__main__":
    doctest.testmod()
    for base_length in range(10):
        print(next(stable_triangles(base_length)))

Answer 1

This code is written very idiomatically, so I don't have any style points.

For speed, the first thing to do is use PyPy. That gives a 5x throughput improvement.

Then one should consider your algorithm. A quick run with cProfile gives

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
  8420076    1.231    0.000    1.231    0.000 p.py:5(<listcomp>)
  1414943    0.966    0.000    0.992    0.000 p.py:14(all_unique)
  1414943    0.557    0.000    2.462    0.000 p.py:8(make_triangle)
  1414943    0.504    0.000    0.504    0.000 p.py:18(<listcomp>)
  8420076    0.477    0.000    1.708    0.000 p.py:4(next_level)
        8    0.275    0.034    4.347    0.543 p.py:23(stable_triangles)
  8421251    0.114    0.000    0.114    0.000 {method 'append' of 'list' objects}
12667944/12667894    0.109    0.000    0.109    0.000 {built-in function len}
  1414943    0.077    0.000    1.610    0.000 p.py:20(is_stable)
  1414943    0.037    0.000    0.541    0.000 p.py:17(flatten1)
        1    0.020    0.020    4.466    4.466 p.py:1(<module>)

This means that next_level is taking a considerable amount of time. I noticed it's better written

[left + right for left, right in zip(level[:-1], level[1:])]

but this doesn't change speed much.

The first useful thing I would consider is interspersing generation and testing to prevent many calls to next_level. This also prevents the need to keep the whole triangle - one only needs the base and current row. This can be done by making triangle yield rows instead.

from itertools import permutations

def next_level(level):
    return [left + right for left, right in zip(level[:-1], level[1:])]

def make_triangle(base):
    row = base

    while row:
        yield row
        row = next_level(row)

def is_stable(triangle):
    seen = set()

    for row in triangle:
        if not seen.isdisjoint(row) or len(row) != len(set(row)):
            return False
        seen.update(row)

    return True

def stable_triangles(base_length):
    for base in permutations(range(base_length * 2), base_length):
        if is_stable(make_triangle(base)):
            yield make_triangle(base)

This is pretty fast comparatively. Not fast enough, though!

---

Here's where we get to the really interesting parts.

Consider some unique triangle:

    9  
  3   6
1   2   4

and the next

     26
    9  17
  3   6  11
1   2   4   7

It's not a coincidence that one is a superset of another. In fact, is must be so. For any non-trivial unique triangle, there are subsets you can find by dropping the left or right hand side.

So consider building them up in this manner, filtering as early as possible. One problem is that we need a value limit to prevent infinite recursion. Luckily, I can choose the same rule you already have to avoid changing behaviours:

def next_level(level):
    return [left + right for left, right in zip(level[:-1], level[1:])]

def make_triangle(base):
    row = base

    while row:
        yield row
        row = next_level(row)

def is_stable(triangle):
    seen = set()

    for row in triangle:
        if not seen.isdisjoint(row) or len(row) != len(set(row)):
            return False
        seen.update(row)

    return True

def stable_triangle_bases(base_length, limit):
    if base_length < 0:
        raise ValueError("Base length must be nonnegative")

    if not base_length:
        yield ()
        return

    for base in stable_triangle_bases(base_length - 1, limit):
        for i in range(limit):
            if is_stable(make_triangle(base + (i,))):
                yield base + (i,)

if __name__ == "__main__":
    for base_length in range(40):
        print(list(next(stable_triangle_bases(base_length, limit=base_length * 2))))

This works really well - up to about length 300 or so, in fact.

---

I suppose we could try to speed this up more, but that seems beside the point. We want to find the minimum. Well, this isn't totally pointless - we can just run a min with a good key on our result. Unfortunately, though, I haven't managed to get past

                                1001
                             523     478
                         309     214     264
                     195     114     100     164
                 125      70      44      56     108
              79      46      24      20      36      72
          48      31      15       9      11      25      47
      27      21      10       5       4       7      18      29
  14      13       8       2       3       1       6      12      17

with that code.

Further, a limit of base_length * 2 is not optimal. Consider the triangle

                              1000
                           485     515
                       273     212     303
                   171     102     110     193
               112      59      43      67     126
            73      39      20      23      44      82
        46      27      12       8      15      29      53
    28      18       9       3       5      10      19      34
17      11       7       2       1       4       6      13      21

!

Raising the limit makes things even slower - an unbounded limit would take forever!

---

One conceptually easy, but extremely effective, change is to apply a heuristic to the search. This says that if we have found a triangle of cost \$k\$ and our new path could not give us a cost lower than that, do not try the path.

I originally avoided this because I could not think of a strong enough heuristic. It turns out, though, that it's really simple. Consider this unfinished triangle:

       ??
     ??  ??
    9  ??  ??
  3   6  ??  ??
1   2   4   x  ??

We are going to place onto x next, for all values that will fit. Fill in the unknowns on the base with the minimal values that fit, pushing larger numbers to the edges.

In this case, we cannot add 1, 2, 3 or 4, so we fill in a 5. Treat the x as a 0 for now:

       38
     19  19
    9  10   9
  3   6   4   5
1   2   4   x   5

Note that the duplicates here are allowed, as removing them requires knowing the value of x.

Further, consider the contribution of x:

       4x
      x  3x
    0   x  2x
  0   0   x   x
0   0   0   x   0

We can add these two triangles together, giving a top value of

4x + 38

as our heuristic. Now when outputting results, we can ignore any branches where the heuristic brings us above (or equal to) that. (If you want all of the best options, use strict inequality. If you only need one, use non-strict inequality.)

Here is some code to try. It is not particularly pretty, but it works and is pretty fast.

from itertools import chain, count


def pascal_row(row):
    val = 1
    yield val
    for k in range(row):
        val *= row - k
        val //= k + 1
        yield val


def make_triangle(row):
    while row:
        yield tuple(row)
        row = [left + right for left, right in zip(row[:-1], row[1:])]

def heuristic_minimum(base_so_far, length):
    used = set().union(*make_triangle(base_so_far))
    unused = (x for x in count(1) if x not in used)

    factors = pascal_row(length - 1)
    raw_cost =  sum(a * b for a, b in zip(base_so_far, factors))
    x_factor = next(factors)
    raw_cost += sum(a * b for a, b in zip(unused, sorted(factors, reverse=True)))

    return x_factor, raw_cost


def is_stable(triangle):
    seen = set()

    for row in triangle:
        if not seen.isdisjoint(row) or len(row) != len(set(row)):
            return False
        seen.update(row)

    return True


def stable_triangle_bases(base_length, upper_bound, full_length):
    if base_length < 0:
        raise ValueError("Base length must be nonnegative")

    if not base_length:
        yield ()
        return

    for base in stable_triangle_bases(base_length - 1, upper_bound, full_length):
        x_factor, min_const_cost = heuristic_minimum(base, full_length)

        for i in count(1):
            cost = i * x_factor + min_const_cost

            if cost >= upper_bound.value:
                break

            if is_stable(make_triangle(base + (i,))):
                if base_length == full_length:
                    upper_bound.value = cost

                yield base + (i,)


class UpperBound:
    def __init__(self):
        self.value = float("inf")


if __name__ == "__main__":
    base_length = 10

    bases = stable_triangle_bases(
        base_length,
        UpperBound(),
        base_length
    )

    for result in bases:
        print(result)

This will print the improving bound as time progresses.

To top this off, here's a nice big triangle:

                                      4497
                                  2301    2196
                              1298    1003    1193
                           778     520     483     710
                       479     299     221     262     448
                   296     183     116     105     157     291
               181     115      68      48      57     100     191
           109      72      43      25      23      34      66     125
        65      44      28      15      10      13      21      45      80
    38      27      17      11       4       6       7      14      31      49
20      18       9       8       3       1       5       2      12      19      30