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
min(stable_triangles(7), key=flattend_sum)
would yield the correct result in a day's time. \$\endgroup\$